Repository: scikit-learn/scikit-learn
Pull request: [MRG+2] OPTICS: add extract_xi method
Related: #12036, #12053, #11677
Fixes #12376

Adding the Xi extraction method to OPTICS.

@espg, I'll need the `predecessor_` to add that final touch to it. I'm still testing the implementation, but it'll be nice if I could get some feedback. (done)

Differences with the [paper](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.407.5572&rep=rep1&type=pdf):
- Definition 11 section 4: comparisons are done with `r(e_U + 1)`, but section 4.3.2 compares values with `r(e_U)`, this implementation takes `r(e_U + 1)` which seems more consistent.
- The article assimes `min_samples` as the minimum size of clusters, which is the default behavior of this implementation, but we also have a `min_cluster_size` to give more freedom to the user.

Problem statement:
OPTICS Return all the clusters found by the extraction algorithm
In OPTICS, the automatic extraction algorithm actually returns all possible clusters, but currently we only provide users with the clusters on the leaves. I guess we can provide users with an option to get all these clusters. (R provides clusters_xi attribute).
ping @espg @adrinjalali (I think it's even more important for extract_Xi)

Hints from issue discussion:
Both `extract_sqlnk` and `extract_xi` provide a hierarchy of clusters. I've already put `return_clusters` as a parameter to `extract_xi` in #12077. We can also extract the cluster hierarchy from the tree constructed by sqlnk.

I'd say we can add `store_cluster_hierarchy` (or something) to the `OPTICS`'s init, and store the hierarchy in `fit` accordingly.
Thanks. I've gone through your PR and this issue is for so-called sqlnk
We've correct our RD and I'll submit a PR to correct CD in a couple of days. At this point, an important thing in my mind is to decide the API design of OPTICS. If we use extract_method, then if users want to extract clusters with different methods, they'll need to refit the model (i.e., compute CD and RD again).
And another thing is the way we extract labels for each sample from the hierarchy structure. In sqlnk, we only consider the leaves, and I doubt whether it's correct. E.g., I guess there might be some unexpected behavior if one leaf of a node is pruned. AFAIK, R dbscan also consider the parent node. I'll need some time to think about it carefully and your insights on these questions are deeply appreciated.
> . If we use extract_method, then if users want to extract clusters with
different methods, they'll need to refit the model (i.e., compute CD and RD
again).

See #12044

The R OPTICS extract_xi has an option where you can tell if it should use the smallest clusters or also use the parent clusters when assigning the labels. The `_xi_cluster` in #12077 puts the clusters from smallest to largest in the reported cluster list, that's why `_extract_xi_labels` works fine. It'll be easy to also include the larger clusters in labeling if we want to. But that can also be seen as adding `include_larger_clusters` or something with default to `False` later as an enhancement without breaking backward compatibility.
Return the leaves but store the hierarchy on the estimator??? Return at a specified max depth according to some parameter?
Or prune to a specified `n_clusters`
It's more like calculate the `labels_` according to the leaves rather than returning the leaves. And yes, to store the hierarchy in the estimator. Does "???" mean you really don't like it? I guess I've seen it in other estimators, like store_covariances or something.

We can start by calculating the `labels_` according to the leaves (or smallest found clusters), then add storing the hierarchy, and helper functions and parameters to fine tune extraction of `labels_` from the hierarchy.
??? is an expression of uncertainty.

We should absolutely provide a way to encapsulate multiple cluster hierarchies. I don't know how to do it, especially within the current sklearn API, but I certainly support enabling it as a feature...

Commit message: fix min_cluster_size as 1

Patch to review:
diff --git a/doc/modules/classes.rst b/doc/modules/classes.rst
index 7464729a52e0d..aa4fa44e5bb12 100644
--- a/doc/modules/classes.rst
+++ b/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
diff --git a/doc/modules/clustering.rst b/doc/modules/clustering.rst
index 9a77f657e2257..81d75fb6515b0 100644
--- a/doc/modules/clustering.rst
+++ b/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
@@ -796,6 +802,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
@@ -814,6 +825,92 @@ by black points below.
    In Proceedings of the 2nd International Conference on Knowledge Discovery
    and Data Mining, Portland, OR, AAAI Press, pp. 226–231. 1996
 
+.. _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
diff --git a/examples/cluster/plot_cluster_comparison.py b/examples/cluster/plot_cluster_comparison.py
index 39d8bca458cc2..69bdaf98d5565 100644
--- a/examples/cluster/plot_cluster_comparison.py
+++ b/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': 30, 'xi': 0.001, '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)
     )
diff --git a/examples/cluster/plot_optics.py b/examples/cluster/plot_optics.py
new file mode 100644
index 0000000000000..1618cbd18530f
--- /dev/null
+++ b/examples/cluster/plot_optics.py
@@ -0,0 +1,101 @@
+"""
+===================================
+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 clustering is first used in its automatic settings, which is the
+:class:`sklearn.cluster.OPTICS` algorithm, and then setting specific
+thresholds on the reachability, which corresponds to DBSCAN.
+We can see that the different clusters of OPTICS 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 numpy as np
+
+import matplotlib.pyplot as plt
+
+# 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()
diff --git a/sklearn/cluster/__init__.py b/sklearn/cluster/__init__.py
index 3f6f7d50d3a4c..da4cfdb6f0734 100644
--- a/sklearn/cluster/__init__.py
+++ b/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',
diff --git a/sklearn/cluster/optics_.py b/sklearn/cluster/optics_.py
index c947d12fb8f8e..38c149533f746 100755
--- a/sklearn/cluster/optics_.py
+++ b/sklearn/cluster/optics_.py
@@ -47,8 +47,8 @@ class OPTICS(BaseEstimator, ClusterMixin):
 
     max_eps : float, optional (default=np.inf)
         The maximum distance between two samples for them to be considered
-        as in the same neighborhood. Default value of "np.inf" will identify
-        clusters across all scales; reducing `max_eps` will result in
+        as in the same neighborhood. Default value of ``np.inf`` will identify
+        clusters across all scales; reducing ``max_eps`` will result in
         shorter run times.
 
     metric : string or callable, optional (default='minkowski')
@@ -86,19 +86,33 @@ class OPTICS(BaseEstimator, ClusterMixin):
     metric_params : dict, optional (default=None)
         Additional keyword arguments for the metric function.
 
-    cluster_method : string, optional (default='dbscan')
+    cluster_method : string, optional (default='xi')
         The extraction method used to extract clusters using the calculated
-        reachability and ordering. Possible values are "dbscan".
+        reachability and ordering. Possible values are "xi" and "dbscan".
 
-    eps : float, optional (default=0.5)
+    eps : float, optional (default=None)
         The maximum distance between two samples for them to be considered
-        as in the same neighborhood. Used ony when `cluster_method='dbscan'`.
-
-    min_cluster_size : int > 1 or float between 0 and 1 (default=0.005)
+        as in the same neighborhood. By default it assumes the same value as
+        ``max_eps``.
+        Used only when ``cluster_method='dbscan'``.
+
+    xi : float, between 0 and 1, optional (default=0.05)
+        Determines the minimum steepness on the reachability plot that
+        constitutes a cluster boundary. For example, an upwards point in the
+        reachability plot is defined by the ratio from one point to its
+        successor being at most 1-xi.
+        Used only when ``cluster_method='xi'``.
+
+    predecessor_correction : bool, optional (default=True)
+        Correct clusters according to the predecessors calculated by OPTICS
+        [2]_. This parameter has minimal effect on most datasets.
+        Used only when ``cluster_method='xi'``.
+
+    min_cluster_size : int > 1 or float between 0 and 1 (default=None)
         Minimum number of samples in an OPTICS cluster, expressed as an
-        absolute number or a fraction of the number of samples (rounded
-        to be at least 2).
-        Used only when `extract_method='xi'`.
+        absolute number or a fraction of the number of samples (rounded to be
+        at least 2). If ``None``, the value of ``min_samples`` is used instead.
+        Used only when ``cluster_method='xi'``.
 
     algorithm : {'auto', 'ball_tree', 'kd_tree', 'brute'}, optional
         Algorithm used to compute the nearest neighbors:
@@ -128,7 +142,8 @@ class OPTICS(BaseEstimator, ClusterMixin):
     ----------
     labels_ : array, shape (n_samples,)
         Cluster labels for each point in the dataset given to fit().
-        Noisy samples are given the label -1.
+        Noisy samples and points which are not included in a leaf cluster
+        of ``cluster_hierarchy_`` are labeled as -1.
 
     reachability_ : array, shape (n_samples,)
         Reachability distances per sample, indexed by object order. Use
@@ -146,6 +161,17 @@ class OPTICS(BaseEstimator, ClusterMixin):
         Point that a sample was reached from, indexed by object order.
         Seed points have a predecessor of -1.
 
+    cluster_hierarchy_ : array, shape (n_clusters, 2)
+        The list of clusters in the form of ``[start, end]`` in each row, with
+        all indices inclusive. The clusters are ordered according to
+        ``(end, -start)`` (ascending) so that larger clusters encompassing
+        smaller clusters come after those smaller ones. Since ``labels_`` does
+        not reflect the hierarchy, usually
+        ``len(cluster_hierarchy_) > np.unique(optics.labels_)``. Please also
+        note that these indices are of the ``ordering_``, i.e.
+        ``X[ordering_][start:end + 1]`` form a cluster.
+        Only available when ``cluster_method='xi'``.
+
     See also
     --------
 
@@ -166,9 +192,9 @@ class OPTICS(BaseEstimator, ClusterMixin):
     """
 
     def __init__(self, min_samples=5, max_eps=np.inf, metric='minkowski', p=2,
-                 metric_params=None, cluster_method='dbscan', eps=0.5,
-                 min_cluster_size=.005, algorithm='auto', leaf_size=30,
-                 n_jobs=None):
+                 metric_params=None, cluster_method='xi', eps=None, xi=0.05,
+                 predecessor_correction=True, min_cluster_size=None,
+                 algorithm='auto', leaf_size=30, n_jobs=None):
 
         self.max_eps = max_eps
         self.min_samples = min_samples
@@ -180,13 +206,15 @@ def __init__(self, min_samples=5, max_eps=np.inf, metric='minkowski', p=2,
         self.leaf_size = leaf_size
         self.cluster_method = cluster_method
         self.eps = eps
+        self.xi = xi
+        self.predecessor_correction = predecessor_correction
         self.n_jobs = n_jobs
 
     def fit(self, X, y=None):
         """Perform OPTICS clustering
 
         Extracts an ordered list of points and reachability distances, and
-        performs initial clustering using `max_eps` distance specified at
+        performs initial clustering using ``max_eps`` distance specified at
         OPTICS object instantiation.
 
         Parameters
@@ -205,16 +233,26 @@ def fit(self, X, y=None):
 
         n_samples = len(X)
 
-        if self.min_cluster_size <= 0 or (self.min_cluster_size !=
-                                          int(self.min_cluster_size)
-                                          and self.min_cluster_size > 1):
+        if self.eps is None:
+            eps = self.max_eps
+        else:
+            eps = self.eps
+
+        if self.min_cluster_size is None:
+            min_cluster_size = self.min_samples
+        else:
+            min_cluster_size = self.min_cluster_size
+
+        if min_cluster_size <= 0 or (min_cluster_size !=
+                                     int(min_cluster_size)
+                                     and min_cluster_size > 1):
             raise ValueError('min_cluster_size must be a positive integer or '
                              'a float between 0 and 1. Got %r' %
-                             self.min_cluster_size)
-        elif self.min_cluster_size > n_samples:
+                             min_cluster_size)
+        elif min_cluster_size > n_samples:
             raise ValueError('min_cluster_size must be no greater than the '
                              'number of samples (%d). Got %d' %
-                             (n_samples, self.min_cluster_size))
+                             (n_samples, min_cluster_size))
 
         if self.min_samples > n_samples:
             raise ValueError("Number of training samples (n_samples=%d) must "
@@ -222,23 +260,15 @@ def fit(self, X, y=None):
                              "used for clustering." %
                              (n_samples, self.min_samples))
 
-        if self.cluster_method not in ['dbscan']:
+        if self.cluster_method not in ['dbscan', 'xi']:
             raise ValueError("cluster_method should be one of"
-                             " 'dbscan', but is %s" %
+                             " 'dbscan' or 'xi' but is %s" %
                              self.cluster_method)
 
         if self.cluster_method == 'dbscan':
-            if self.eps > self.max_eps:
+            if eps > self.max_eps:
                 raise ValueError('Specify an epsilon smaller than %s. Got %s.'
-                                 % (self.max_eps, self.eps))
-
-            if self.eps * 5.0 > self.max_eps * 1.05:
-                warnings.warn(
-                    "Warning, max_eps (%s) is close to eps (%s): "
-                    "Output may be unstable." % (self.max_eps, self.eps),
-                    RuntimeWarning, stacklevel=2)
-                # Stability warning is documented in cluster_optics_dbscan
-                # method...
+                                 % (self.max_eps, eps))
 
         (self.ordering_, self.core_distances_, self.reachability_,
          self.predecessor_) = compute_optics_graph(
@@ -246,13 +276,25 @@ def fit(self, X, y=None):
              leaf_size=self.leaf_size, metric=self.metric,
              metric_params=self.metric_params, p=self.p, n_jobs=self.n_jobs,
              max_eps=self.max_eps)
+        # Start all points as noise ##
+        self.labels_ = np.full(n_samples, -1, dtype=int)
 
         # Extract clusters from the calculated orders and reachability
-        if self.cluster_method == 'dbscan':
+        if self.cluster_method == 'xi':
+            labels_, clusters_ = cluster_optics_xi(
+                self.reachability_,
+                self.predecessor_,
+                self.ordering_,
+                self.min_samples,
+                min_cluster_size,
+                self.xi,
+                self.predecessor_correction)
+            self.cluster_hierarchy_ = clusters_
+        elif self.cluster_method == 'dbscan':
             labels_ = cluster_optics_dbscan(self.reachability_,
                                             self.core_distances_,
                                             self.ordering_,
-                                            self.eps)
+                                            eps)
 
         self.labels_ = labels_
         return self
@@ -489,29 +531,28 @@ def _set_reach_dist(core_distances_, reachability_, predecessor_,
     predecessor_[unproc[improved]] = point_index
 
 
-def cluster_optics_dbscan(reachability, core_distances, ordering, eps=0.5):
+def cluster_optics_dbscan(reachability, core_distances, ordering, eps):
     """Performs DBSCAN extraction for an arbitrary epsilon.
 
-    Extracting the clusters runs in linear time. Note that if the `max_eps`
-    OPTICS parameter was set to < inf for extracting reachability and ordering
-    arrays, DBSCAN extractions will be unstable for `eps` values close to
-    `max_eps`. Setting `eps` < (`max_eps` / 5.0) will guarantee extraction
-    parity with DBSCAN.
+    Extracting the clusters runs in linear time. Note that this results in
+    ``labels_`` which are close to a `DBSCAN` with similar settings and
+    ``eps``, only if ``eps`` is close to ``max_eps``.
 
     Parameters
     ----------
     reachability : array, shape (n_samples,)
-        Reachability distances calculated by OPTICS (`reachability_`)
+        Reachability distances calculated by OPTICS (``reachability_``)
 
     core_distances : array, shape (n_samples,)
-        Distances at which points become core (`core_distances_`)
+        Distances at which points become core (``core_distances_``)
 
     ordering : array, shape (n_samples,)
-        OPTICS ordered point indices (`ordering_`)
+        OPTICS ordered point indices (``ordering_``)
 
-    eps : float, optional (default=0.5)
-        DBSCAN `eps` parameter. Must be set to < `max_eps`. Results
-        will be close to DBSCAN algorithm if `eps` is < (`max_eps` / 5)
+    eps : float
+        DBSCAN ``eps`` parameter. Must be set to < ``max_eps``. Results
+        will be close to DBSCAN algorithm if ``eps`` and ``max_eps`` are close
+        to one another.
 
     Returns
     -------
@@ -527,3 +568,382 @@ def cluster_optics_dbscan(reachability, core_distances, ordering, eps=0.5):
     labels[ordering] = np.cumsum(far_reach[ordering] & near_core[ordering]) - 1
     labels[far_reach & ~near_core] = -1
     return labels
+
+
+def cluster_optics_xi(reachability, predecessor, ordering, min_samples,
+                      min_cluster_size=None, xi=0.05,
+                      predecessor_correction=False):
+    """Automatically extract clusters according to the Xi-steep method.
+
+    Parameters
+    ----------
+    reachability : array, shape (n_samples,)
+        Reachability distances calculated by OPTICS (`reachability_`)
+
+    predecessor : array, shape (n_samples,)
+        Predecessors calculated by OPTICS.
+
+    ordering : array, shape (n_samples,)
+        OPTICS ordered point indices (`ordering_`)
+
+    min_samples : integer
+       The same as the min_samples given to OPTICS. Up and down steep
+       regions can't have more then `min_samples` consecutive non-steep
+       points.
+
+    min_cluster_size : int > 1 or float between 0 and 1 (default=None)
+        Minimum number of samples in an OPTICS cluster, expressed as an
+        absolute number or a fraction of the number of samples (rounded to be
+        at least 2). If ``None``, the value of ``min_samples`` is used instead.
+
+    xi : float, between 0 and 1, optional (default=0.05)
+        Determines the minimum steepness on the reachability plot that
+        constitutes a cluster boundary. For example, an upwards point in the
+        reachability plot is defined by the ratio from one point to its
+        successor being at most 1-xi.
+
+    predecessor_correction : bool, optional (default=True)
+        Correct clusters based on the calculated predecessors.
+
+    Returns
+    -------
+    labels : array, shape (n_samples)
+        The labels assigned to samples. Points which are not included
+        in any cluster are labeled as -1.
+
+    clusters : array, shape (n_clusters, 2)
+        The list of clusters in the form of ``[start, end]`` in each row, with
+        all indices inclusive. The clusters are ordered according to ``(end,
+        -start)`` (ascending) so that larger clusters encompassing smaller
+        clusters come after such nested smaller clusters. Since ``labels`` does
+        not reflect the hierarchy, usually ``len(clusters) >
+        np.unique(labels)``.
+    """
+    if min_cluster_size is None:
+        min_cluster_size = min_samples
+
+    clusters = _xi_cluster(reachability[ordering], predecessor[ordering], xi,
+                           min_samples, min_cluster_size,
+                           predecessor_correction)
+    labels = _extract_xi_labels(ordering, clusters)
+    return labels, clusters
+
+
+class _Area:
+    """An (upward or downward) area.
+
+    Attributes
+    ----------
+    start : integer
+        The start of the region.
+
+    end : integer
+        The end of the region.
+
+    maximum : float
+        The maximum reachability in this region, which is the
+        start of the region for a downward area and the end of
+        the region for an upward area.
+
+    mib : float
+        Maximum in between value, i.e. the maximum value between
+        the end of a steep down area and the current index. It is
+        irrelevant for steep up areas.
+    """
+    def __init__(self, start, end, maximum, mib):
+        self.start = start
+        self.end = end
+        # maximum is the beginning of a downward area
+        # and the end of an upward area
+        self.maximum = maximum
+        # maximum in between (ref Sec. 4.3.2)
+        self.mib = mib
+
+    def __repr__(self):
+        return ("start: %s, end: %s, max: %4g, mib: %4g" %
+                (self.start, self.end, self.maximum, self.mib))
+
+
+def _extend_region(steep_point, xward_point, start, min_samples):
+    """Extend the area until it's maximal.
+
+    It's the same function for both upward and downward reagions, depending on
+    the given input parameters. To extend an upward reagion,
+    ``steep_point=steep_upward`` and ``xward_point=downward`` are expected, and
+    to extend a downward region, ``steep_point=steep_downward`` and
+    ``xward_point=upward``.
+
+    Parameters
+    ----------
+    steep_point : bool array, shape (n_samples)
+        True if the point is steep downward (upward).
+
+    xward_point : bool array, shape (n_samples)
+        True if the point is an upward (respectively downward) point.
+
+    start : integer
+        The start of the xward region.
+
+    min_samples : integer
+       The same as the min_samples given to OPTICS. Up and down steep
+       regions can't have more then ``min_samples`` consecutive non-steep
+       points.
+
+    Returns
+    -------
+    index : integer
+        The current index iterating over all the samples, i.e. where we are up
+        to in our search.
+
+    end : integer
+        The end of the region, which can be behind the index. The region
+        includes the ``end`` index.
+    """
+    n_samples = len(steep_point)
+    non_xward_points = 0
+    index = start
+    end = start
+    # find a maximal area
+    while index < n_samples:
+        if steep_point[index]:
+            non_xward_points = 0
+            end = index
+        elif not xward_point[index]:
+            # it's not a steep point, but still goes up.
+            non_xward_points += 1
+            # region should include no more than min_samples consecutive
+            # non steep xward points.
+            if non_xward_points > min_samples:
+                break
+        else:
+            return min(index, n_samples - 1), end
+
+        index += 1
+    return min(index + 1, n_samples - 1), end
+
+
+def _extend_upward_region(steep_upward_mask, downward_mask, start,
+                          min_samples):
+    return _extend_region(steep_upward_mask, downward_mask,
+                          start, min_samples)
+
+
+def _extend_downward_region(steep_downward_mask, upward_mask, start,
+                            min_samples):
+    return _extend_region(steep_downward_mask, upward_mask,
+                          start, min_samples)
+
+
+def _update_filter_sdas(sdas, mib, xi_complement):
+    """Update steep down areas (SDAs) using the new maximum in between (mib)
+    value, and the given complement of xi, i.e. ``1 - xi``.
+    """
+    if np.isinf(mib):
+        return []
+    res = [sda for sda in sdas if mib <= sda.maximum * xi_complement]
+    for sda in res:
+        sda.mib = max(sda.mib, mib)
+    return res
+
+
+def _correct_predecessor(reachability_plot, predecessor, s, e):
+    """Correct for predecessors.
+
+    Applies Algorithm 2 of [1]_.
+
+    Input parameters are ordered by the computer OPTICS ordering.
+
+    .. [1] Schubert, Erich, Michael Gertz.
+       "Improving the Cluster Structure Extracted from OPTICS Plots." Proc. of
+       the Conference "Lernen, Wissen, Daten, Analysen" (LWDA) (2018): 318-329.
+    """
+
+    while s < e:
+        if reachability_plot[s] > reachability_plot[e]:
+            return s, e
+        for i in range(s, e):
+            if predecessor[e] == i:
+                return s, e
+        e -= 1
+    return None, None
+
+
+def _xi_cluster(reachability_plot, predecessor, xi, min_samples,
+                min_cluster_size, predecessor_correction):
+    """Automatically extract clusters according to the Xi-steep method.
+
+    This is rouphly an implementation of Figure 19 of the OPTICS paper.
+
+    Parameters
+    ----------
+    reachability_plot : array, shape (n_samples)
+        The reachability plot, i.e. reachability ordered according to
+        the calculated ordering, all computed by OPTICS.
+
+    predecessor : array, shape (n_samples)
+        Predecessors ordered according to the calculated ordering.
+
+    xi : float, between 0 and 1
+        Determines the minimum steepness on the reachability plot that
+        constitutes a cluster boundary. For example, an upwards point in the
+        reachability plot is defined by the ratio from one point to its
+        successor being at most 1-xi.
+
+    min_samples : integer
+       The same as the min_samples given to OPTICS. Up and down steep
+       regions can't have more then `min_samples` consecutive non-steep
+       points.
+
+    min_cluster_size : int > 1 or float between 0 and 1
+        Minimum number of samples in an OPTICS cluster, expressed as an
+        absolute number or a fraction of the number of samples (rounded
+        to be at least 2).
+
+    predecessor_correction : bool
+        Correct clusters based on the calculated predecessors.
+
+    Returns
+    -------
+    clusters : array, shape (n_clusters, 2)
+        The list of clusters in the form of [start, end] in each row, with all
+        indices inclusive. The clusters are ordered in a way that larger
+        clusters encompassing smaller clusters come after those smaller
+        clusters.
+    """
+
+    # all indices are inclusive (specially at the end)
+    n_samples = len(reachability_plot)
+    # add an inf to the end of reachability plot
+    # but process the data only until the last actual point
+    # this is fine since the last point is considered upward anyway
+    reachability_plot = np.hstack((reachability_plot, np.inf))
+
+    if min_cluster_size <= 1:
+        min_cluster_size = max(2, min_cluster_size * n_samples)
+
+    xi_complement = 1 - xi
+    sdas = []  # steep down areas, introduced in section 4.3.2 of the paper
+    clusters = []
+    index = 0
+    mib = 0.  # maximum in between
+
+    ratio = reachability_plot[:-1] / reachability_plot[1:]
+    steep_upward = ratio <= xi_complement
+    steep_downward = ratio >= 1 / xi_complement
+    downward = ratio > 1
+    upward = ratio < 1
+
+    # the following loop is is almost exactly as Figure 19 of the paper.
+    # it jumps over the areas which are not either steep down or up areas
+    for steep_index in iter(np.flatnonzero(steep_upward | steep_downward)):
+        # just continue if steep_index has been a part of a discovered xward
+        # area.
+        if steep_index < index:
+            continue
+
+        mib = max(mib, np.max(reachability_plot[index:steep_index + 1]))
+        index = steep_index
+
+        # check if a steep downward area starts
+        if steep_downward[index]:
+            sdas = _update_filter_sdas(sdas, mib, xi_complement)
+            D_start = index
+            index, end = _extend_downward_region(steep_downward, upward,
+                                                 D_start, min_samples)
+            D = _Area(start=D_start, end=end,
+                      maximum=reachability_plot[D_start], mib=0.)
+            sdas.append(D)
+            mib = reachability_plot[end + 1]
+
+        elif steep_upward[index]:
+            sdas = _update_filter_sdas(sdas, mib, xi_complement)
+            U_start = index
+            index, end = _extend_upward_region(steep_upward, downward, U_start,
+                                               min_samples)
+            U = _Area(start=U_start, end=end, maximum=reachability_plot[end],
+                      mib=-1)
+            mib = reachability_plot[end + 1]
+
+            U_clusters = []
+            for D in sdas:
+                c_start = D.start
+                c_end = min(U.end, n_samples - 1)
+
+                # line (**), sc2*
+                if reachability_plot[c_end + 1] * xi_complement < D.mib:
+                    continue
+
+                # Definition 11: criterion 4
+                if D.maximum * xi_complement >= reachability_plot[c_end + 1]:
+                    # Find the first index from the left side which is almost
+                    # at the same level as the end of the detected cluster.
+                    while (reachability_plot[c_start + 1] >
+                           reachability_plot[c_end + 1]
+                           and c_start < c_end):
+                        c_start += 1
+                elif reachability_plot[c_end] * xi_complement >= D.maximum:
+                    # Find the first index from the right side which is almost
+                    # at the same level as the beginning of the detected
+                    # cluster.
+                    while (reachability_plot[c_end] > D.maximum
+                           and c_end > c_start):
+                        c_end -= 1
+
+                # predecessor correction
+                c_start, c_end = _correct_predecessor(reachability_plot,
+                                                      predecessor, c_start,
+                                                      c_end)
+                if c_start is None:
+                    continue
+
+                # Definition 11: criterion 3.a
+                if c_end - c_start + 1 < min_cluster_size:
+                    continue
+
+                # Definition 11: criterion 1
+                if c_start > D.end:
+                    continue
+
+                # Definition 11: criterion 2
+                if c_end < U.start:
+                    continue
+
+                U_clusters.append((c_start, c_end))
+
+            # add smaller clusters first.
+            U_clusters.reverse()
+            clusters.extend(U_clusters)
+
+            if U.end == n_samples - 1:
+                break
+
+    return np.array(clusters)
+
+
+def _extract_xi_labels(ordering, clusters):
+    """Extracts the labels from the clusters returned by `_xi_cluster`.
+    We rely on the fact that clusters are stored
+    with the smaller clusters coming before the larger ones.
+
+    Parameters
+    ----------
+    ordering : array, shape (n_samples)
+        The ordering of points calculated by OPTICS
+
+    clusters : array, shape (n_clusters, 2)
+        List of clusters i.e. (start, end) tuples,
+        as returned by `_xi_cluster`.
+
+    Returns
+    -------
+    labels : array, shape (n_samples)
+    """
+
+    labels = np.full(len(ordering), -1, dtype=int)
+    label = 0
+    for c in clusters:
+        if not np.any(labels[c[0]:(c[1] + 1)] != -1):
+            labels[c[0]:(c[1] + 1)] = label
+            label += 1
+    labels[ordering] = labels.copy()
+    return labels
diff --git a/sklearn/cluster/tests/test_optics.py b/sklearn/cluster/tests/test_optics.py
old mode 100755
new mode 100644
index 37166e7f8dd20..b421a9f9b1a0f
--- a/sklearn/cluster/tests/test_optics.py
+++ b/sklearn/cluster/tests/test_optics.py
@@ -6,10 +6,13 @@
 import pytest
 
 from sklearn.datasets.samples_generator import make_blobs
-from sklearn.cluster.optics_ import OPTICS
+from sklearn.cluster.optics_ import (OPTICS,
+                                     _extend_region,
+                                     _extract_xi_labels)
 from sklearn.metrics.cluster import contingency_matrix
 from sklearn.metrics.pairwise import pairwise_distances
 from sklearn.cluster.dbscan_ import DBSCAN
+from sklearn.utils import shuffle
 from sklearn.utils.testing import assert_equal, assert_warns
 from sklearn.utils.testing import assert_array_equal
 from sklearn.utils.testing import assert_raise_message
@@ -29,6 +32,114 @@
 X = np.vstack((C1, C2, C3, C4, C5, C6))
 
 
+@pytest.mark.parametrize(
+    ('r_plot', 'end', 'index'),
+    [[[10, 8.9, 8.8, 8.7, 7, 10], 3, 4],
+     [[10, 8.9, 8.8, 8.7, 8.6, 7, 10], 0, 4],
+     [[10, 8.9, 8.8, 8.7, 7, 6, np.inf], 4, 5],
+     [[10, 8.9, 8.8, 8.7, 7, 6, np.inf], 4, 5],
+     ])
+def test_extend_downward(r_plot, end, index):
+    r_plot = np.array(r_plot)
+    ratio = r_plot[:-1] / r_plot[1:]
+    steep_downward = ratio >= 1 / .9
+    upward = ratio < 1
+
+    i, e = _extend_region(steep_downward, upward, 0, 2)
+    assert e == end
+    assert i == index
+
+
+@pytest.mark.parametrize(
+    ('r_plot', 'end', 'index'),
+    [[[1, 2, 2.1, 2.2, 4, 8, 8, np.inf], 6, 6],
+     [[1, 2, 2.1, 2.2, 2.3, 4, 8, 8, np.inf], 0, 4],
+     [[1, 2, 2.1, 2, np.inf], 0, 2],
+     [[1, 2, 2.1, np.inf], 2, 2],
+     ])
+def test_extend_upward(r_plot, end, index):
+    r_plot = np.array(r_plot)
+    ratio = r_plot[:-1] / r_plot[1:]
+    steep_upward = ratio <= .9
+    downward = ratio > 1
+
+    i, e = _extend_region(steep_upward, downward, 0, 2)
+    assert e == end
+    assert i == index
+
+
+@pytest.mark.parametrize(
+    ('ordering', 'clusters', 'expected'),
+    [[[0, 1, 2, 3], [[0, 1], [2, 3]], [0, 0, 1, 1]],
+     [[0, 1, 2, 3], [[0, 1], [3, 3]], [0, 0, -1, 1]],
+     [[0, 1, 2, 3], [[0, 1], [3, 3], [0, 3]], [0, 0, -1, 1]],
+     [[3, 1, 2, 0], [[0, 1], [3, 3], [0, 3]], [1, 0, -1, 0]],
+     ])
+def test_the_extract_xi_labels(ordering, clusters, expected):
+    labels = _extract_xi_labels(ordering, clusters)
+
+    assert_array_equal(labels, expected)
+
+
+def test_extract_xi():
+    # small and easy test (no clusters around other clusters)
+    # but with a clear noise data.
+    rng = np.random.RandomState(0)
+    n_points_per_cluster = 5
+
+    C1 = [-5, -2] + .8 * rng.randn(n_points_per_cluster, 2)
+    C2 = [4, -1] + .1 * rng.randn(n_points_per_cluster, 2)
+    C3 = [1, -2] + .2 * rng.randn(n_points_per_cluster, 2)
+    C4 = [-2, 3] + .3 * rng.randn(n_points_per_cluster, 2)
+    C5 = [3, -2] + .6 * rng.randn(n_points_per_cluster, 2)
+    C6 = [5, 6] + .2 * rng.randn(n_points_per_cluster, 2)
+
+    X = np.vstack((C1, C2, C3, C4, C5, np.array([[100, 100]]), C6))
+    expected_labels = np.r_[[2] * 5, [0] * 5, [1] * 5, [3] * 5, [1] * 5,
+                            -1, [4] * 5]
+    X, expected_labels = shuffle(X, expected_labels, random_state=rng)
+
+    clust = OPTICS(min_samples=3, min_cluster_size=2,
+                   max_eps=np.inf, cluster_method='xi',
+                   xi=0.1).fit(X)
+    assert_array_equal(clust.labels_, expected_labels)
+
+    X = np.vstack((C1, C2, C3, C4, C5, np.array([[100, 100]] * 2), C6))
+    expected_labels = np.r_[[1] * 5, [3] * 5, [2] * 5, [0] * 5, [2] * 5,
+                            -1, -1, [4] * 5]
+    X, expected_labels = shuffle(X, expected_labels, random_state=rng)
+
+    clust = OPTICS(min_samples=3, min_cluster_size=3,
+                   max_eps=np.inf, cluster_method='xi',
+                   xi=0.1).fit(X)
+    assert_array_equal(clust.labels_, expected_labels)
+
+    C1 = [[0, 0], [0, 0.1], [0, -.1], [0.1, 0]]
+    C2 = [[10, 10], [10, 9], [10, 11], [9, 10]]
+    C3 = [[100, 100], [100, 90], [100, 110], [90, 100]]
+    X = np.vstack((C1, C2, C3))
+    expected_labels = np.r_[[0] * 4, [1] * 4, [2] * 4]
+    X, expected_labels = shuffle(X, expected_labels, random_state=rng)
+
+    clust = OPTICS(min_samples=2, min_cluster_size=2,
+                   max_eps=np.inf, cluster_method='xi',
+                   xi=0.04).fit(X)
+    assert_array_equal(clust.labels_, expected_labels)
+
+
+def test_cluster_hierarchy_():
+    n_points_per_cluster = 100
+    C1 = [0, 0] + 2 * rng.randn(n_points_per_cluster, 2)
+    C2 = [0, 0] + 10 * rng.randn(n_points_per_cluster, 2)
+    X = np.vstack((C1, C2))
+    X = shuffle(X, random_state=0)
+
+    clusters = OPTICS(min_samples=20, xi=.1).fit(X).cluster_hierarchy_
+    assert clusters.shape == (2, 2)
+    diff = np.sum(clusters - np.array([[0, 99], [0, 199]]))
+    assert diff / len(X) < 0.05
+
+
 def test_correct_number_of_clusters():
     # in 'auto' mode
 
@@ -36,7 +147,7 @@ def test_correct_number_of_clusters():
     X = generate_clustered_data(n_clusters=n_clusters)
     # Parameters chosen specifically for this task.
     # Compute OPTICS
-    clust = OPTICS(max_eps=5.0 * 6.0, min_samples=4)
+    clust = OPTICS(max_eps=5.0 * 6.0, min_samples=4, xi=.1)
     clust.fit(X)
     # number of clusters, ignoring noise if present
     n_clusters_1 = len(set(clust.labels_)) - int(-1 in clust.labels_)
@@ -64,7 +175,7 @@ def test_minimum_number_of_sample_check():
 
     # Compute OPTICS
     X = [[1, 1]]
-    clust = OPTICS(max_eps=5.0 * 0.3, min_samples=10)
+    clust = OPTICS(max_eps=5.0 * 0.3, min_samples=10, min_cluster_size=1)
 
     # Run the fit
     assert_raise_message(ValueError, msg, clust.fit, X)
@@ -96,7 +207,7 @@ def test_bad_reachability():
 
 
 def test_close_extract():
-    # Test extract where extraction eps is close to scaled epsPrime
+    # Test extract where extraction eps is close to scaled max_eps
 
     centers = [[1, 1], [-1, -1], [1, -1]]
     X, labels_true = make_blobs(n_samples=750, centers=centers,
@@ -104,9 +215,7 @@ def test_close_extract():
 
     # Compute OPTICS
     clust = OPTICS(max_eps=1.0, cluster_method='dbscan',
-                   eps=0.3, min_samples=10)
-    # check warning when centers are passed
-    assert_warns(RuntimeWarning, clust.fit, X)
+                   eps=0.3, min_samples=10).fit(X)
     # Cluster ordering starts at 0; max cluster label = 2 is 3 clusters
     assert_equal(max(clust.labels_), 2)
 
@@ -138,6 +247,31 @@ def test_dbscan_optics_parity(eps, min_samples):
     assert percent_mismatch <= 0.05
 
 
+def test_min_samples_edge_case():
+    C1 = [[0, 0], [0, 0.1], [0, -.1]]
+    C2 = [[10, 10], [10, 9], [10, 11]]
+    C3 = [[100, 100], [100, 96], [100, 106]]
+    X = np.vstack((C1, C2, C3))
+
+    expected_labels = np.r_[[0] * 3, [1] * 3, [2] * 3]
+    clust = OPTICS(min_samples=3,
+                   max_eps=7, cluster_method='xi',
+                   xi=0.04).fit(X)
+    assert_array_equal(clust.labels_, expected_labels)
+
+    expected_labels = np.r_[[0] * 3, [1] * 3, [-1] * 3]
+    clust = OPTICS(min_samples=3,
+                   max_eps=3, cluster_method='xi',
+                   xi=0.04).fit(X)
+    assert_array_equal(clust.labels_, expected_labels)
+
+    expected_labels = np.r_[[-1] * 9]
+    clust = OPTICS(min_samples=4,
+                   max_eps=3, cluster_method='xi',
+                   xi=0.04).fit(X)
+    assert_array_equal(clust.labels_, expected_labels)
+
+
 # try arbitrary minimum sizes
 @pytest.mark.parametrize('min_cluster_size', range(2, X.shape[0] // 10, 23))
 def test_min_cluster_size(min_cluster_size):
diff --git a/sklearn/utils/estimator_checks.py b/sklearn/utils/estimator_checks.py
index 73c98ea4685be..2ee43657fa02a 100644
--- a/sklearn/utils/estimator_checks.py
+++ b/sklearn/utils/estimator_checks.py
@@ -867,6 +867,10 @@ def check_fit2d_1sample(name, estimator_orig):
 
     set_random_state(estimator, 1)
 
+    # min_cluster_size cannot be less than the data size for OPTICS.
+    if name == 'OPTICS':
+        estimator.set_params(min_cluster_size=1)
+
     msgs = ["1 sample", "n_samples = 1", "n_samples=1", "one sample",
             "1 class", "one class"]
 
@@ -1316,6 +1320,9 @@ def check_clustering(name, clusterer_orig, readonly_memmap=False):
         clusterer.set_params(preference=-100)
         clusterer.set_params(max_iter=100)
 
+    if name == 'OPTICS':
+        clusterer.set_params(min_samples=5)
+
     # fit
     clusterer.fit(X)
     # with lists


Code review comments:
diff_hunk: @@ -527,3 +568,382 @@ def cluster_optics_dbscan(reachability, core_distances, ordering, eps=0.5):
     labels[ordering] = np.cumsum(far_reach[ordering] & near_core[ordering]) - 1
     labels[far_reach & ~near_core] = -1
     return labels
+
+
+def cluster_optics_xi(reachability, predecessor, ordering, min_samples,
+                      min_cluster_size=None, xi=0.05,
+                      predecessor_correction=False):
path: sklearn/cluster/optics_.py
text: @user1:
predecessor_correction=True?

diff_hunk: @@ -205,54 +233,68 @@ def fit(self, X, y=None):
 
         n_samples = len(X)
 
-        if self.min_cluster_size <= 0 or (self.min_cluster_size !=
-                                          int(self.min_cluster_size)
-                                          and self.min_cluster_size > 1):
+        if self.eps is None:
+            eps = self.max_eps
+        else:
+            eps = self.eps
+
+        if self.min_cluster_size is None:
path: sklearn/cluster/optics_.py
text: @user1:
move these things to ``if self.cluster_method == 'xi':``?

diff_hunk: @@ -205,54 +233,68 @@ def fit(self, X, y=None):
 
         n_samples = len(X)
 
-        if self.min_cluster_size <= 0 or (self.min_cluster_size !=
-                                          int(self.min_cluster_size)
-                                          and self.min_cluster_size > 1):
+        if self.eps is None:
+            eps = self.max_eps
+        else:
+            eps = self.eps
+
+        if self.min_cluster_size is None:
+            min_cluster_size = self.min_samples
+        else:
+            min_cluster_size = self.min_cluster_size
+
+        if min_cluster_size <= 0 or (min_cluster_size !=
+                                     int(min_cluster_size)
+                                     and min_cluster_size > 1):
             raise ValueError('min_cluster_size must be a positive integer or '
                              'a float between 0 and 1. Got %r' %
-                             self.min_cluster_size)
-        elif self.min_cluster_size > n_samples:
+                             min_cluster_size)
+        elif min_cluster_size > n_samples:
             raise ValueError('min_cluster_size must be no greater than the '
                              'number of samples (%d). Got %d' %
-                             (n_samples, self.min_cluster_size))
+                             (n_samples, min_cluster_size))
 
         if self.min_samples > n_samples:
             raise ValueError("Number of training samples (n_samples=%d) must "
                              "be greater than min_samples (min_samples=%d) "
                              "used for clustering." %
                              (n_samples, self.min_samples))
 
-        if self.cluster_method not in ['dbscan']:
+        if self.cluster_method not in ['dbscan', 'xi']:
             raise ValueError("cluster_method should be one of"
-                             " 'dbscan', but is %s" %
+                             " 'dbscan' or 'xi' but is %s" %
                              self.cluster_method)
 
         if self.cluster_method == 'dbscan':
path: sklearn/cluster/optics_.py
text: @user1:
move these things to ``elif self.cluster_method == 'dbscan':``?

diff_hunk: @@ -527,3 +568,382 @@ def cluster_optics_dbscan(reachability, core_distances, ordering, eps=0.5):
     labels[ordering] = np.cumsum(far_reach[ordering] & near_core[ordering]) - 1
     labels[far_reach & ~near_core] = -1
     return labels
+
+
+def cluster_optics_xi(reachability, predecessor, ordering, min_samples,
+                      min_cluster_size=None, xi=0.05,
+                      predecessor_correction=False):
+    """Automatically extract clusters according to the Xi-steep method.
+
+    Parameters
+    ----------
+    reachability : array, shape (n_samples,)
+        Reachability distances calculated by OPTICS (`reachability_`)
+
+    predecessor : array, shape (n_samples,)
+        Predecessors calculated by OPTICS.
+
+    ordering : array, shape (n_samples,)
+        OPTICS ordered point indices (`ordering_`)
+
+    min_samples : integer
+       The same as the min_samples given to OPTICS. Up and down steep
+       regions can't have more then `min_samples` consecutive non-steep
+       points.
+
+    min_cluster_size : int > 1 or float between 0 and 1 (default=None)
+        Minimum number of samples in an OPTICS cluster, expressed as an
+        absolute number or a fraction of the number of samples (rounded to be
+        at least 2). If ``None``, the value of ``min_samples`` is used instead.
+
+    xi : float, between 0 and 1, optional (default=0.05)
+        Determines the minimum steepness on the reachability plot that
+        constitutes a cluster boundary. For example, an upwards point in the
+        reachability plot is defined by the ratio from one point to its
+        successor being at most 1-xi.
+
+    predecessor_correction : bool, optional (default=True)
+        Correct clusters based on the calculated predecessors.
+
+    Returns
+    -------
+    labels : array, shape (n_samples)
+        The labels assigned to samples. Points which are not included
+        in any cluster are labeled as -1.
+
+    clusters : array, shape (n_clusters, 2)
+        The list of clusters in the form of ``[start, end]`` in each row, with
+        all indices inclusive. The clusters are ordered according to ``(end,
+        -start)`` (ascending) so that larger clusters encompassing smaller
+        clusters come after such nested smaller clusters. Since ``labels`` does
+        not reflect the hierarchy, usually ``len(clusters) >
+        np.unique(labels)``.
+    """
+    if min_cluster_size is None:
+        min_cluster_size = min_samples
+
+    clusters = _xi_cluster(reachability[ordering], predecessor[ordering], xi,
+                           min_samples, min_cluster_size,
+                           predecessor_correction)
+    labels = _extract_xi_labels(ordering, clusters)
+    return labels, clusters
+
+
+class _Area:
+    """An (upward or downward) area.
+
+    Attributes
+    ----------
+    start : integer
+        The start of the region.
+
+    end : integer
+        The end of the region.
+
+    maximum : float
+        The maximum reachability in this region, which is the
+        start of the region for a downward area and the end of
+        the region for an upward area.
+
+    mib : float
+        Maximum in between value, i.e. the maximum value between
+        the end of a steep down area and the current index. It is
+        irrelevant for steep up areas.
+    """
+    def __init__(self, start, end, maximum, mib):
+        self.start = start
+        self.end = end
+        # maximum is the beginning of a downward area
+        # and the end of an upward area
+        self.maximum = maximum
+        # maximum in between (ref Sec. 4.3.2)
+        self.mib = mib
+
+    def __repr__(self):
+        return ("start: %s, end: %s, max: %4g, mib: %4g" %
+                (self.start, self.end, self.maximum, self.mib))
+
+
+def _extend_region(steep_point, xward_point, start, min_samples):
+    """Extend the area until it's maximal.
+
+    It's the same function for both upward and downward reagions, depending on
+    the given input parameters. To extend an upward reagion,
+    ``steep_point=steep_upward`` and ``xward_point=downward`` are expected, and
path: sklearn/cluster/optics_.py
text: @user1:
Be more specific, it's difficult for others to understand the meaning of ``steep_point=steep_upward`` without reading through the code carefully.

@author:
I'm not sure how it could be improved, please provide alternative text.

@user1:
``steep_upward`` -> ``whether a point is a steep upward point``. The main point to avoid using the name of a variable.
Trivial issue, feel free to ignore.

@author:
added some clarification.

diff_hunk: @@ -205,54 +233,68 @@ def fit(self, X, y=None):
 
         n_samples = len(X)
 
-        if self.min_cluster_size <= 0 or (self.min_cluster_size !=
-                                          int(self.min_cluster_size)
-                                          and self.min_cluster_size > 1):
+        if self.eps is None:
path: sklearn/cluster/optics_.py
text: @user1:
move these things to ``elif self.cluster_method == 'dbscan':``?

diff_hunk: @@ -527,3 +568,382 @@ def cluster_optics_dbscan(reachability, core_distances, ordering, eps=0.5):
     labels[ordering] = np.cumsum(far_reach[ordering] & near_core[ordering]) - 1
     labels[far_reach & ~near_core] = -1
     return labels
+
+
+def cluster_optics_xi(reachability, predecessor, ordering, min_samples,
+                      min_cluster_size=None, xi=0.05,
+                      predecessor_correction=False):
+    """Automatically extract clusters according to the Xi-steep method.
+
+    Parameters
+    ----------
+    reachability : array, shape (n_samples,)
+        Reachability distances calculated by OPTICS (`reachability_`)
+
+    predecessor : array, shape (n_samples,)
+        Predecessors calculated by OPTICS.
+
+    ordering : array, shape (n_samples,)
+        OPTICS ordered point indices (`ordering_`)
+
+    min_samples : integer
+       The same as the min_samples given to OPTICS. Up and down steep
+       regions can't have more then `min_samples` consecutive non-steep
+       points.
+
+    min_cluster_size : int > 1 or float between 0 and 1 (default=None)
+        Minimum number of samples in an OPTICS cluster, expressed as an
+        absolute number or a fraction of the number of samples (rounded to be
+        at least 2). If ``None``, the value of ``min_samples`` is used instead.
+
+    xi : float, between 0 and 1, optional (default=0.05)
+        Determines the minimum steepness on the reachability plot that
+        constitutes a cluster boundary. For example, an upwards point in the
+        reachability plot is defined by the ratio from one point to its
+        successor being at most 1-xi.
+
+    predecessor_correction : bool, optional (default=True)
+        Correct clusters based on the calculated predecessors.
+
+    Returns
+    -------
+    labels : array, shape (n_samples)
+        The labels assigned to samples. Points which are not included
+        in any cluster are labeled as -1.
+
+    clusters : array, shape (n_clusters, 2)
+        The list of clusters in the form of ``[start, end]`` in each row, with
+        all indices inclusive. The clusters are ordered according to ``(end,
+        -start)`` (ascending) so that larger clusters encompassing smaller
+        clusters come after such nested smaller clusters. Since ``labels`` does
+        not reflect the hierarchy, usually ``len(clusters) >
+        np.unique(labels)``.
+    """
+    if min_cluster_size is None:
+        min_cluster_size = min_samples
+
+    clusters = _xi_cluster(reachability[ordering], predecessor[ordering], xi,
+                           min_samples, min_cluster_size,
+                           predecessor_correction)
+    labels = _extract_xi_labels(ordering, clusters)
+    return labels, clusters
+
+
+class _Area:
+    """An (upward or downward) area.
+
+    Attributes
+    ----------
+    start : integer
+        The start of the region.
+
+    end : integer
+        The end of the region.
+
+    maximum : float
+        The maximum reachability in this region, which is the
+        start of the region for a downward area and the end of
+        the region for an upward area.
+
+    mib : float
+        Maximum in between value, i.e. the maximum value between
+        the end of a steep down area and the current index. It is
+        irrelevant for steep up areas.
+    """
+    def __init__(self, start, end, maximum, mib):
+        self.start = start
+        self.end = end
+        # maximum is the beginning of a downward area
+        # and the end of an upward area
+        self.maximum = maximum
+        # maximum in between (ref Sec. 4.3.2)
+        self.mib = mib
+
+    def __repr__(self):
+        return ("start: %s, end: %s, max: %4g, mib: %4g" %
+                (self.start, self.end, self.maximum, self.mib))
+
+
+def _extend_region(steep_point, xward_point, start, min_samples):
+    """Extend the area until it's maximal.
+
+    It's the same function for both upward and downward reagions, depending on
+    the given input parameters. To extend an upward reagion,
+    ``steep_point=steep_upward`` and ``xward_point=downward`` are expected, and
+    to extend a downward region, ``steep_point=steep_downward`` and
+    ``xward_point=upward``.
+
+    Parameters
+    ----------
+    steep_point : bool array, shape (n_samples)
+        True if the point is steep downward (upward).
+
+    xward_point : bool array, shape (n_samples)
+        True if the point is an upward (respectively downward) point.
+
+    start : integer
+        The start of the xward region.
+
+    min_samples : integer
+       The same as the min_samples given to OPTICS. Up and down steep
+       regions can't have more then ``min_samples`` consecutive non-steep
+       points.
+
+    Returns
+    -------
+    index : integer
+        The current index iterating over all the samples, i.e. where we are up
+        to in our search.
+
+    end : integer
+        The end of the region, which can be behind the index. The region
+        includes the ``end`` index.
+    """
+    n_samples = len(steep_point)
+    non_xward_points = 0
+    index = start
+    end = start
+    # find a maximal area
+    while index < n_samples:
+        if steep_point[index]:
+            non_xward_points = 0
+            end = index
+        elif not xward_point[index]:
+            # it's not a steep point, but still goes up.
+            non_xward_points += 1
+            # region should include no more than min_samples consecutive
+            # non steep xward points.
+            if non_xward_points > min_samples:
+                break
+        else:
+            return min(index, n_samples - 1), end
+
+        index += 1
+    return min(index + 1, n_samples - 1), end
path: sklearn/cluster/optics_.py
text: @user1:
``end`` seems redundant, the extra points you go through can't be steep_upward points or steep_downward points?

@author:
`end` is not redundant since it can be on the left side of `index` when the function returns. The loop continues for `downward`/`upward` points, even if they're not _steep_.

@user1:
I can't understand your point, a steep area has to been maximum?

@author:
could you please rephrase the question?

@user2:
As far as I'm concerned, `end` is the main thing this function is trying to calculate. `index` is returned merely to update the index at which we are searching, and may only be an efficiency bonus to avoid duplicated work. The one cannot be derived directly from the other.

@user2:
The next index to explore is only for efficiency, I think. It is certainly not explicit from the algorithm in the paper. It is just a convenience to, as in the description, indicate "where we are up to in our search."

diff_hunk: @@ -527,3 +568,382 @@ def cluster_optics_dbscan(reachability, core_distances, ordering, eps=0.5):
     labels[ordering] = np.cumsum(far_reach[ordering] & near_core[ordering]) - 1
     labels[far_reach & ~near_core] = -1
     return labels
+
+
+def cluster_optics_xi(reachability, predecessor, ordering, min_samples,
+                      min_cluster_size=None, xi=0.05,
+                      predecessor_correction=False):
+    """Automatically extract clusters according to the Xi-steep method.
+
+    Parameters
+    ----------
+    reachability : array, shape (n_samples,)
+        Reachability distances calculated by OPTICS (`reachability_`)
+
+    predecessor : array, shape (n_samples,)
+        Predecessors calculated by OPTICS.
+
+    ordering : array, shape (n_samples,)
+        OPTICS ordered point indices (`ordering_`)
+
+    min_samples : integer
+       The same as the min_samples given to OPTICS. Up and down steep
+       regions can't have more then `min_samples` consecutive non-steep
+       points.
+
+    min_cluster_size : int > 1 or float between 0 and 1 (default=None)
+        Minimum number of samples in an OPTICS cluster, expressed as an
+        absolute number or a fraction of the number of samples (rounded to be
+        at least 2). If ``None``, the value of ``min_samples`` is used instead.
+
+    xi : float, between 0 and 1, optional (default=0.05)
+        Determines the minimum steepness on the reachability plot that
+        constitutes a cluster boundary. For example, an upwards point in the
+        reachability plot is defined by the ratio from one point to its
+        successor being at most 1-xi.
+
+    predecessor_correction : bool, optional (default=True)
+        Correct clusters based on the calculated predecessors.
+
+    Returns
+    -------
+    labels : array, shape (n_samples)
+        The labels assigned to samples. Points which are not included
+        in any cluster are labeled as -1.
+
+    clusters : array, shape (n_clusters, 2)
+        The list of clusters in the form of ``[start, end]`` in each row, with
+        all indices inclusive. The clusters are ordered according to ``(end,
+        -start)`` (ascending) so that larger clusters encompassing smaller
+        clusters come after such nested smaller clusters. Since ``labels`` does
+        not reflect the hierarchy, usually ``len(clusters) >
+        np.unique(labels)``.
+    """
+    if min_cluster_size is None:
+        min_cluster_size = min_samples
+
+    clusters = _xi_cluster(reachability[ordering], predecessor[ordering], xi,
+                           min_samples, min_cluster_size,
+                           predecessor_correction)
+    labels = _extract_xi_labels(ordering, clusters)
+    return labels, clusters
+
+
+class _Area:
+    """An (upward or downward) area.
+
+    Attributes
+    ----------
+    start : integer
+        The start of the region.
+
+    end : integer
+        The end of the region.
+
+    maximum : float
+        The maximum reachability in this region, which is the
+        start of the region for a downward area and the end of
+        the region for an upward area.
+
+    mib : float
+        Maximum in between value, i.e. the maximum value between
+        the end of a steep down area and the current index. It is
+        irrelevant for steep up areas.
+    """
+    def __init__(self, start, end, maximum, mib):
+        self.start = start
+        self.end = end
+        # maximum is the beginning of a downward area
+        # and the end of an upward area
+        self.maximum = maximum
+        # maximum in between (ref Sec. 4.3.2)
+        self.mib = mib
+
+    def __repr__(self):
+        return ("start: %s, end: %s, max: %4g, mib: %4g" %
+                (self.start, self.end, self.maximum, self.mib))
+
+
+def _extend_region(steep_point, xward_point, start, min_samples):
+    """Extend the area until it's maximal.
+
+    It's the same function for both upward and downward reagions, depending on
+    the given input parameters. To extend an upward reagion,
+    ``steep_point=steep_upward`` and ``xward_point=downward`` are expected, and
+    to extend a downward region, ``steep_point=steep_downward`` and
+    ``xward_point=upward``.
+
+    Parameters
+    ----------
+    steep_point : bool array, shape (n_samples)
+        True if the point is steep downward (upward).
+
+    xward_point : bool array, shape (n_samples)
+        True if the point is an upward (respectively downward) point.
+
+    start : integer
+        The start of the xward region.
+
+    min_samples : integer
+       The same as the min_samples given to OPTICS. Up and down steep
+       regions can't have more then ``min_samples`` consecutive non-steep
+       points.
+
+    Returns
+    -------
+    index : integer
+        The current index iterating over all the samples, i.e. where we are up
+        to in our search.
+
+    end : integer
+        The end of the region, which can be behind the index. The region
+        includes the ``end`` index.
+    """
+    n_samples = len(steep_point)
+    non_xward_points = 0
+    index = start
+    end = start
+    # find a maximal area
+    while index < n_samples:
+        if steep_point[index]:
+            non_xward_points = 0
+            end = index
+        elif not xward_point[index]:
+            # it's not a steep point, but still goes up.
+            non_xward_points += 1
+            # region should include no more than min_samples consecutive
+            # non steep xward points.
+            if non_xward_points > min_samples:
+                break
+        else:
+            return min(index, n_samples - 1), end
+
+        index += 1
+    return min(index + 1, n_samples - 1), end
+
+
+def _extend_upward_region(steep_upward_mask, downward_mask, start,
+                          min_samples):
+    return _extend_region(steep_upward_mask, downward_mask,
+                          start, min_samples)
+
+
+def _extend_downward_region(steep_downward_mask, upward_mask, start,
+                            min_samples):
+    return _extend_region(steep_downward_mask, upward_mask,
+                          start, min_samples)
+
+
+def _update_filter_sdas(sdas, mib, xi_complement):
+    """Update steep down areas (SDAs) using the new maximum in between (mib)
+    value, and the given complement of xi, i.e. ``1 - xi``.
+    """
+    if np.isinf(mib):
+        return []
+    res = [sda for sda in sdas if mib <= sda.maximum * xi_complement]
+    for sda in res:
+        sda.mib = max(sda.mib, mib)
+    return res
+
+
+def _correct_predecessor(reachability_plot, predecessor, s, e):
+    """Correct for predecessors.
+
+    Applies Algorithm 2 of [1]_.
+
+    Input parameters are ordered by the computer OPTICS ordering.
+
+    .. [1] Schubert, Erich, Michael Gertz.
+       "Improving the Cluster Structure Extracted from OPTICS Plots." Proc. of
+       the Conference "Lernen, Wissen, Daten, Analysen" (LWDA) (2018): 318-329.
+    """
+
+    while s < e:
+        if reachability_plot[s] > reachability_plot[e]:
+            return s, e
+        for i in range(s, e):
+            if predecessor[e] == i:
+                return s, e
+        e -= 1
+    return None, None
+
+
+def _xi_cluster(reachability_plot, predecessor, xi, min_samples,
+                min_cluster_size, predecessor_correction):
+    """Automatically extract clusters according to the Xi-steep method.
+
+    This is rouphly an implementation of Figure 19 of the OPTICS paper.
+
+    Parameters
+    ----------
+    reachability_plot : array, shape (n_samples)
+        The reachability plot, i.e. reachability ordered according to
+        the calculated ordering, all computed by OPTICS.
+
+    predecessor : array, shape (n_samples)
+        Predecessors ordered according to the calculated ordering.
+
+    xi : float, between 0 and 1
+        Determines the minimum steepness on the reachability plot that
+        constitutes a cluster boundary. For example, an upwards point in the
+        reachability plot is defined by the ratio from one point to its
+        successor being at most 1-xi.
+
+    min_samples : integer
+       The same as the min_samples given to OPTICS. Up and down steep
+       regions can't have more then `min_samples` consecutive non-steep
+       points.
+
+    min_cluster_size : int > 1 or float between 0 and 1
+        Minimum number of samples in an OPTICS cluster, expressed as an
+        absolute number or a fraction of the number of samples (rounded
+        to be at least 2).
+
+    predecessor_correction : bool
+        Correct clusters based on the calculated predecessors.
+
+    Returns
+    -------
+    clusters : array, shape (n_clusters, 2)
+        The list of clusters in the form of [start, end] in each row, with all
+        indices inclusive. The clusters are ordered in a way that larger
+        clusters encompassing smaller clusters come after those smaller
+        clusters.
+    """
+
+    # all indices are inclusive (specially at the end)
+    n_samples = len(reachability_plot)
+    # add an inf to the end of reachability plot
+    # but process the data only until the last actual point
+    # this is fine since the last point is considered upward anyway
+    reachability_plot = np.hstack((reachability_plot, np.inf))
+
+    if min_cluster_size <= 1:
+        min_cluster_size = max(2, min_cluster_size * n_samples)
+
+    xi_complement = 1 - xi
+    sdas = []  # steep down areas, introduced in section 4.3.2 of the paper
+    clusters = []
+    index = 0
+    mib = 0.  # maximum in between
+
+    ratio = reachability_plot[:-1] / reachability_plot[1:]
+    steep_upward = ratio <= xi_complement
+    steep_downward = ratio >= 1 / xi_complement
+    downward = ratio > 1
+    upward = ratio < 1
+
+    # the following loop is is almost exactly as Figure 19 of the paper.
+    # it jumps over the areas which are not either steep down or up areas
+    for steep_index in iter(np.flatnonzero(steep_upward | steep_downward)):
+        # just continue if steep_index has been a part of a discovered xward
+        # area.
+        if steep_index < index:
+            continue
+
+        mib = max(mib, np.max(reachability_plot[index:steep_index + 1]))
+        index = steep_index
+
+        # check if a steep downward area starts
+        if steep_downward[index]:
+            sdas = _update_filter_sdas(sdas, mib, xi_complement)
+            D_start = index
+            index, end = _extend_downward_region(steep_downward, upward,
+                                                 D_start, min_samples)
+            D = _Area(start=D_start, end=end,
+                      maximum=reachability_plot[D_start], mib=0.)
+            sdas.append(D)
+            mib = reachability_plot[end + 1]
+
+        elif steep_upward[index]:
+            sdas = _update_filter_sdas(sdas, mib, xi_complement)
+            U_start = index
+            index, end = _extend_upward_region(steep_upward, downward, U_start,
+                                               min_samples)
+            U = _Area(start=U_start, end=end, maximum=reachability_plot[end],
+                      mib=-1)
+            mib = reachability_plot[end + 1]
+
+            U_clusters = []
+            for D in sdas:
+                c_start = D.start
+                c_end = min(U.end, n_samples - 1)
+
+                # line (**), sc2*
+                if reachability_plot[c_end + 1] * xi_complement < D.mib:
+                    continue
+
+                # Definition 11: criterion 4
+                if D.maximum * xi_complement >= reachability_plot[c_end + 1]:
+                    # Find the first index from the left side which is almost
+                    # at the same level as the end of the detected cluster.
+                    while (reachability_plot[c_start + 1] >
+                           reachability_plot[c_end + 1]
+                           and c_start < c_end):
+                        c_start += 1
+                elif reachability_plot[c_end] * xi_complement >= D.maximum:
+                    # Find the first index from the right side which is almost
+                    # at the same level as the beginning of the detected
+                    # cluster.
+                    while (reachability_plot[c_end] > D.maximum
+                           and c_end > c_start):
+                        c_end -= 1
+
+                # predecessor correction
+                c_start, c_end = _correct_predecessor(reachability_plot,
path: sklearn/cluster/optics_.py
text: @user1:
you ignore parameter ``predecessor_correction``?

diff_hunk: @@ -527,3 +568,382 @@ def cluster_optics_dbscan(reachability, core_distances, ordering, eps=0.5):
     labels[ordering] = np.cumsum(far_reach[ordering] & near_core[ordering]) - 1
     labels[far_reach & ~near_core] = -1
     return labels
+
+
+def cluster_optics_xi(reachability, predecessor, ordering, min_samples,
+                      min_cluster_size=None, xi=0.05,
+                      predecessor_correction=False):
+    """Automatically extract clusters according to the Xi-steep method.
+
+    Parameters
+    ----------
+    reachability : array, shape (n_samples,)
+        Reachability distances calculated by OPTICS (`reachability_`)
+
+    predecessor : array, shape (n_samples,)
+        Predecessors calculated by OPTICS.
+
+    ordering : array, shape (n_samples,)
+        OPTICS ordered point indices (`ordering_`)
+
+    min_samples : integer
+       The same as the min_samples given to OPTICS. Up and down steep
+       regions can't have more then `min_samples` consecutive non-steep
+       points.
+
+    min_cluster_size : int > 1 or float between 0 and 1 (default=None)
+        Minimum number of samples in an OPTICS cluster, expressed as an
+        absolute number or a fraction of the number of samples (rounded to be
+        at least 2). If ``None``, the value of ``min_samples`` is used instead.
+
+    xi : float, between 0 and 1, optional (default=0.05)
+        Determines the minimum steepness on the reachability plot that
+        constitutes a cluster boundary. For example, an upwards point in the
+        reachability plot is defined by the ratio from one point to its
+        successor being at most 1-xi.
+
+    predecessor_correction : bool, optional (default=True)
+        Correct clusters based on the calculated predecessors.
+
+    Returns
+    -------
+    labels : array, shape (n_samples)
+        The labels assigned to samples. Points which are not included
+        in any cluster are labeled as -1.
+
+    clusters : array, shape (n_clusters, 2)
+        The list of clusters in the form of ``[start, end]`` in each row, with
+        all indices inclusive. The clusters are ordered according to ``(end,
+        -start)`` (ascending) so that larger clusters encompassing smaller
+        clusters come after such nested smaller clusters. Since ``labels`` does
+        not reflect the hierarchy, usually ``len(clusters) >
+        np.unique(labels)``.
+    """
+    if min_cluster_size is None:
+        min_cluster_size = min_samples
+
+    clusters = _xi_cluster(reachability[ordering], predecessor[ordering], xi,
+                           min_samples, min_cluster_size,
+                           predecessor_correction)
+    labels = _extract_xi_labels(ordering, clusters)
+    return labels, clusters
+
+
+class _Area:
+    """An (upward or downward) area.
+
+    Attributes
+    ----------
+    start : integer
+        The start of the region.
+
+    end : integer
+        The end of the region.
+
+    maximum : float
+        The maximum reachability in this region, which is the
+        start of the region for a downward area and the end of
+        the region for an upward area.
+
+    mib : float
+        Maximum in between value, i.e. the maximum value between
+        the end of a steep down area and the current index. It is
+        irrelevant for steep up areas.
+    """
+    def __init__(self, start, end, maximum, mib):
+        self.start = start
+        self.end = end
+        # maximum is the beginning of a downward area
+        # and the end of an upward area
+        self.maximum = maximum
+        # maximum in between (ref Sec. 4.3.2)
+        self.mib = mib
+
+    def __repr__(self):
+        return ("start: %s, end: %s, max: %4g, mib: %4g" %
+                (self.start, self.end, self.maximum, self.mib))
+
+
+def _extend_region(steep_point, xward_point, start, min_samples):
+    """Extend the area until it's maximal.
+
+    It's the same function for both upward and downward reagions, depending on
+    the given input parameters. To extend an upward reagion,
+    ``steep_point=steep_upward`` and ``xward_point=downward`` are expected, and
+    to extend a downward region, ``steep_point=steep_downward`` and
+    ``xward_point=upward``.
+
+    Parameters
+    ----------
+    steep_point : bool array, shape (n_samples)
+        True if the point is steep downward (upward).
+
+    xward_point : bool array, shape (n_samples)
+        True if the point is an upward (respectively downward) point.
+
+    start : integer
+        The start of the xward region.
+
+    min_samples : integer
+       The same as the min_samples given to OPTICS. Up and down steep
+       regions can't have more then ``min_samples`` consecutive non-steep
+       points.
+
+    Returns
+    -------
+    index : integer
+        The current index iterating over all the samples, i.e. where we are up
+        to in our search.
+
+    end : integer
+        The end of the region, which can be behind the index. The region
+        includes the ``end`` index.
+    """
+    n_samples = len(steep_point)
+    non_xward_points = 0
+    index = start
+    end = start
+    # find a maximal area
+    while index < n_samples:
+        if steep_point[index]:
+            non_xward_points = 0
+            end = index
+        elif not xward_point[index]:
+            # it's not a steep point, but still goes up.
+            non_xward_points += 1
+            # region should include no more than min_samples consecutive
+            # non steep xward points.
+            if non_xward_points > min_samples:
+                break
+        else:
+            return min(index, n_samples - 1), end
+
+        index += 1
+    return min(index + 1, n_samples - 1), end
+
+
+def _extend_upward_region(steep_upward_mask, downward_mask, start,
+                          min_samples):
+    return _extend_region(steep_upward_mask, downward_mask,
+                          start, min_samples)
+
+
+def _extend_downward_region(steep_downward_mask, upward_mask, start,
+                            min_samples):
+    return _extend_region(steep_downward_mask, upward_mask,
+                          start, min_samples)
+
+
+def _update_filter_sdas(sdas, mib, xi_complement):
+    """Update steep down areas (SDAs) using the new maximum in between (mib)
+    value, and the given complement of xi, i.e. ``1 - xi``.
+    """
+    if np.isinf(mib):
+        return []
+    res = [sda for sda in sdas if mib <= sda.maximum * xi_complement]
+    for sda in res:
+        sda.mib = max(sda.mib, mib)
+    return res
+
+
+def _correct_predecessor(reachability_plot, predecessor, s, e):
+    """Correct for predecessors.
+
+    Applies Algorithm 2 of [1]_.
+
+    Input parameters are ordered by the computer OPTICS ordering.
+
+    .. [1] Schubert, Erich, Michael Gertz.
+       "Improving the Cluster Structure Extracted from OPTICS Plots." Proc. of
+       the Conference "Lernen, Wissen, Daten, Analysen" (LWDA) (2018): 318-329.
+    """
+
+    while s < e:
+        if reachability_plot[s] > reachability_plot[e]:
+            return s, e
+        for i in range(s, e):
+            if predecessor[e] == i:
+                return s, e
+        e -= 1
+    return None, None
+
+
+def _xi_cluster(reachability_plot, predecessor, xi, min_samples,
+                min_cluster_size, predecessor_correction):
+    """Automatically extract clusters according to the Xi-steep method.
+
+    This is rouphly an implementation of Figure 19 of the OPTICS paper.
+
+    Parameters
+    ----------
+    reachability_plot : array, shape (n_samples)
+        The reachability plot, i.e. reachability ordered according to
+        the calculated ordering, all computed by OPTICS.
+
+    predecessor : array, shape (n_samples)
+        Predecessors ordered according to the calculated ordering.
+
+    xi : float, between 0 and 1
+        Determines the minimum steepness on the reachability plot that
+        constitutes a cluster boundary. For example, an upwards point in the
+        reachability plot is defined by the ratio from one point to its
+        successor being at most 1-xi.
+
+    min_samples : integer
+       The same as the min_samples given to OPTICS. Up and down steep
+       regions can't have more then `min_samples` consecutive non-steep
+       points.
+
+    min_cluster_size : int > 1 or float between 0 and 1
+        Minimum number of samples in an OPTICS cluster, expressed as an
+        absolute number or a fraction of the number of samples (rounded
+        to be at least 2).
+
+    predecessor_correction : bool
+        Correct clusters based on the calculated predecessors.
+
+    Returns
+    -------
+    clusters : array, shape (n_clusters, 2)
+        The list of clusters in the form of [start, end] in each row, with all
+        indices inclusive. The clusters are ordered in a way that larger
+        clusters encompassing smaller clusters come after those smaller
+        clusters.
+    """
+
+    # all indices are inclusive (specially at the end)
+    n_samples = len(reachability_plot)
+    # add an inf to the end of reachability plot
+    # but process the data only until the last actual point
+    # this is fine since the last point is considered upward anyway
+    reachability_plot = np.hstack((reachability_plot, np.inf))
+
+    if min_cluster_size <= 1:
+        min_cluster_size = max(2, min_cluster_size * n_samples)
+
+    xi_complement = 1 - xi
+    sdas = []  # steep down areas, introduced in section 4.3.2 of the paper
+    clusters = []
+    index = 0
+    mib = 0.  # maximum in between
+
+    ratio = reachability_plot[:-1] / reachability_plot[1:]
+    steep_upward = ratio <= xi_complement
+    steep_downward = ratio >= 1 / xi_complement
+    downward = ratio > 1
+    upward = ratio < 1
+
+    # the following loop is is almost exactly as Figure 19 of the paper.
+    # it jumps over the areas which are not either steep down or up areas
+    for steep_index in iter(np.flatnonzero(steep_upward | steep_downward)):
+        # just continue if steep_index has been a part of a discovered xward
+        # area.
+        if steep_index < index:
+            continue
+
+        mib = max(mib, np.max(reachability_plot[index:steep_index + 1]))
+        index = steep_index
+
+        # check if a steep downward area starts
+        if steep_downward[index]:
+            sdas = _update_filter_sdas(sdas, mib, xi_complement)
+            D_start = index
+            index, end = _extend_downward_region(steep_downward, upward,
+                                                 D_start, min_samples)
+            D = _Area(start=D_start, end=end,
+                      maximum=reachability_plot[D_start], mib=0.)
+            sdas.append(D)
+            mib = reachability_plot[end + 1]
+
+        elif steep_upward[index]:
+            sdas = _update_filter_sdas(sdas, mib, xi_complement)
+            U_start = index
+            index, end = _extend_upward_region(steep_upward, downward, U_start,
+                                               min_samples)
+            U = _Area(start=U_start, end=end, maximum=reachability_plot[end],
path: sklearn/cluster/optics_.py
text: @user1:
We don't need to create an object here?

diff_hunk: @@ -527,3 +568,382 @@ def cluster_optics_dbscan(reachability, core_distances, ordering, eps=0.5):
     labels[ordering] = np.cumsum(far_reach[ordering] & near_core[ordering]) - 1
     labels[far_reach & ~near_core] = -1
     return labels
+
+
+def cluster_optics_xi(reachability, predecessor, ordering, min_samples,
+                      min_cluster_size=None, xi=0.05,
+                      predecessor_correction=False):
+    """Automatically extract clusters according to the Xi-steep method.
+
+    Parameters
+    ----------
+    reachability : array, shape (n_samples,)
+        Reachability distances calculated by OPTICS (`reachability_`)
+
+    predecessor : array, shape (n_samples,)
+        Predecessors calculated by OPTICS.
+
+    ordering : array, shape (n_samples,)
+        OPTICS ordered point indices (`ordering_`)
+
+    min_samples : integer
+       The same as the min_samples given to OPTICS. Up and down steep
+       regions can't have more then `min_samples` consecutive non-steep
+       points.
+
+    min_cluster_size : int > 1 or float between 0 and 1 (default=None)
+        Minimum number of samples in an OPTICS cluster, expressed as an
+        absolute number or a fraction of the number of samples (rounded to be
+        at least 2). If ``None``, the value of ``min_samples`` is used instead.
+
+    xi : float, between 0 and 1, optional (default=0.05)
+        Determines the minimum steepness on the reachability plot that
+        constitutes a cluster boundary. For example, an upwards point in the
+        reachability plot is defined by the ratio from one point to its
+        successor being at most 1-xi.
+
+    predecessor_correction : bool, optional (default=True)
+        Correct clusters based on the calculated predecessors.
+
+    Returns
+    -------
+    labels : array, shape (n_samples)
+        The labels assigned to samples. Points which are not included
+        in any cluster are labeled as -1.
+
+    clusters : array, shape (n_clusters, 2)
+        The list of clusters in the form of ``[start, end]`` in each row, with
+        all indices inclusive. The clusters are ordered according to ``(end,
+        -start)`` (ascending) so that larger clusters encompassing smaller
+        clusters come after such nested smaller clusters. Since ``labels`` does
+        not reflect the hierarchy, usually ``len(clusters) >
+        np.unique(labels)``.
+    """
+    if min_cluster_size is None:
+        min_cluster_size = min_samples
+
+    clusters = _xi_cluster(reachability[ordering], predecessor[ordering], xi,
+                           min_samples, min_cluster_size,
+                           predecessor_correction)
+    labels = _extract_xi_labels(ordering, clusters)
+    return labels, clusters
+
+
+class _Area:
+    """An (upward or downward) area.
+
+    Attributes
+    ----------
+    start : integer
+        The start of the region.
+
+    end : integer
+        The end of the region.
+
+    maximum : float
+        The maximum reachability in this region, which is the
+        start of the region for a downward area and the end of
+        the region for an upward area.
+
+    mib : float
+        Maximum in between value, i.e. the maximum value between
+        the end of a steep down area and the current index. It is
+        irrelevant for steep up areas.
+    """
+    def __init__(self, start, end, maximum, mib):
+        self.start = start
+        self.end = end
+        # maximum is the beginning of a downward area
+        # and the end of an upward area
+        self.maximum = maximum
+        # maximum in between (ref Sec. 4.3.2)
+        self.mib = mib
+
+    def __repr__(self):
path: sklearn/cluster/optics_.py
text: @user1:
remove this?

diff_hunk: @@ -527,3 +568,382 @@ def cluster_optics_dbscan(reachability, core_distances, ordering, eps=0.5):
     labels[ordering] = np.cumsum(far_reach[ordering] & near_core[ordering]) - 1
     labels[far_reach & ~near_core] = -1
     return labels
+
+
+def cluster_optics_xi(reachability, predecessor, ordering, min_samples,
+                      min_cluster_size=None, xi=0.05,
+                      predecessor_correction=False):
+    """Automatically extract clusters according to the Xi-steep method.
+
+    Parameters
+    ----------
+    reachability : array, shape (n_samples,)
+        Reachability distances calculated by OPTICS (`reachability_`)
+
+    predecessor : array, shape (n_samples,)
+        Predecessors calculated by OPTICS.
+
+    ordering : array, shape (n_samples,)
+        OPTICS ordered point indices (`ordering_`)
+
+    min_samples : integer
+       The same as the min_samples given to OPTICS. Up and down steep
+       regions can't have more then `min_samples` consecutive non-steep
+       points.
+
+    min_cluster_size : int > 1 or float between 0 and 1 (default=None)
+        Minimum number of samples in an OPTICS cluster, expressed as an
+        absolute number or a fraction of the number of samples (rounded to be
+        at least 2). If ``None``, the value of ``min_samples`` is used instead.
+
+    xi : float, between 0 and 1, optional (default=0.05)
+        Determines the minimum steepness on the reachability plot that
+        constitutes a cluster boundary. For example, an upwards point in the
+        reachability plot is defined by the ratio from one point to its
+        successor being at most 1-xi.
+
+    predecessor_correction : bool, optional (default=True)
+        Correct clusters based on the calculated predecessors.
+
+    Returns
+    -------
+    labels : array, shape (n_samples)
+        The labels assigned to samples. Points which are not included
+        in any cluster are labeled as -1.
+
+    clusters : array, shape (n_clusters, 2)
+        The list of clusters in the form of ``[start, end]`` in each row, with
+        all indices inclusive. The clusters are ordered according to ``(end,
+        -start)`` (ascending) so that larger clusters encompassing smaller
+        clusters come after such nested smaller clusters. Since ``labels`` does
+        not reflect the hierarchy, usually ``len(clusters) >
+        np.unique(labels)``.
+    """
+    if min_cluster_size is None:
+        min_cluster_size = min_samples
+
+    clusters = _xi_cluster(reachability[ordering], predecessor[ordering], xi,
+                           min_samples, min_cluster_size,
+                           predecessor_correction)
+    labels = _extract_xi_labels(ordering, clusters)
+    return labels, clusters
+
+
+class _Area:
+    """An (upward or downward) area.
+
+    Attributes
+    ----------
+    start : integer
+        The start of the region.
+
+    end : integer
+        The end of the region.
+
+    maximum : float
+        The maximum reachability in this region, which is the
+        start of the region for a downward area and the end of
+        the region for an upward area.
+
+    mib : float
+        Maximum in between value, i.e. the maximum value between
+        the end of a steep down area and the current index. It is
+        irrelevant for steep up areas.
+    """
+    def __init__(self, start, end, maximum, mib):
+        self.start = start
+        self.end = end
+        # maximum is the beginning of a downward area
+        # and the end of an upward area
+        self.maximum = maximum
+        # maximum in between (ref Sec. 4.3.2)
+        self.mib = mib
+
+    def __repr__(self):
+        return ("start: %s, end: %s, max: %4g, mib: %4g" %
+                (self.start, self.end, self.maximum, self.mib))
+
+
+def _extend_region(steep_point, xward_point, start, min_samples):
+    """Extend the area until it's maximal.
+
+    It's the same function for both upward and downward reagions, depending on
+    the given input parameters. To extend an upward reagion,
+    ``steep_point=steep_upward`` and ``xward_point=downward`` are expected, and
+    to extend a downward region, ``steep_point=steep_downward`` and
+    ``xward_point=upward``.
+
+    Parameters
+    ----------
+    steep_point : bool array, shape (n_samples)
+        True if the point is steep downward (upward).
+
+    xward_point : bool array, shape (n_samples)
+        True if the point is an upward (respectively downward) point.
+
+    start : integer
+        The start of the xward region.
+
+    min_samples : integer
+       The same as the min_samples given to OPTICS. Up and down steep
+       regions can't have more then ``min_samples`` consecutive non-steep
+       points.
+
+    Returns
+    -------
+    index : integer
+        The current index iterating over all the samples, i.e. where we are up
+        to in our search.
+
+    end : integer
+        The end of the region, which can be behind the index. The region
+        includes the ``end`` index.
+    """
+    n_samples = len(steep_point)
+    non_xward_points = 0
+    index = start
+    end = start
+    # find a maximal area
+    while index < n_samples:
+        if steep_point[index]:
+            non_xward_points = 0
+            end = index
+        elif not xward_point[index]:
+            # it's not a steep point, but still goes up.
+            non_xward_points += 1
+            # region should include no more than min_samples consecutive
+            # non steep xward points.
+            if non_xward_points > min_samples:
+                break
+        else:
+            return min(index, n_samples - 1), end
+
+        index += 1
+    return min(index + 1, n_samples - 1), end
+
+
+def _extend_upward_region(steep_upward_mask, downward_mask, start,
+                          min_samples):
+    return _extend_region(steep_upward_mask, downward_mask,
+                          start, min_samples)
+
+
+def _extend_downward_region(steep_downward_mask, upward_mask, start,
+                            min_samples):
+    return _extend_region(steep_downward_mask, upward_mask,
+                          start, min_samples)
+
+
+def _update_filter_sdas(sdas, mib, xi_complement):
+    """Update steep down areas (SDAs) using the new maximum in between (mib)
+    value, and the given complement of xi, i.e. ``1 - xi``.
+    """
+    if np.isinf(mib):
+        return []
+    res = [sda for sda in sdas if mib <= sda.maximum * xi_complement]
+    for sda in res:
+        sda.mib = max(sda.mib, mib)
+    return res
+
+
+def _correct_predecessor(reachability_plot, predecessor, s, e):
+    """Correct for predecessors.
+
+    Applies Algorithm 2 of [1]_.
+
+    Input parameters are ordered by the computer OPTICS ordering.
+
+    .. [1] Schubert, Erich, Michael Gertz.
+       "Improving the Cluster Structure Extracted from OPTICS Plots." Proc. of
+       the Conference "Lernen, Wissen, Daten, Analysen" (LWDA) (2018): 318-329.
+    """
+
+    while s < e:
+        if reachability_plot[s] > reachability_plot[e]:
+            return s, e
+        for i in range(s, e):
+            if predecessor[e] == i:
+                return s, e
+        e -= 1
+    return None, None
+
+
+def _xi_cluster(reachability_plot, predecessor, xi, min_samples,
+                min_cluster_size, predecessor_correction):
+    """Automatically extract clusters according to the Xi-steep method.
+
+    This is rouphly an implementation of Figure 19 of the OPTICS paper.
+
+    Parameters
+    ----------
+    reachability_plot : array, shape (n_samples)
+        The reachability plot, i.e. reachability ordered according to
+        the calculated ordering, all computed by OPTICS.
+
+    predecessor : array, shape (n_samples)
+        Predecessors ordered according to the calculated ordering.
+
+    xi : float, between 0 and 1
+        Determines the minimum steepness on the reachability plot that
+        constitutes a cluster boundary. For example, an upwards point in the
+        reachability plot is defined by the ratio from one point to its
+        successor being at most 1-xi.
+
+    min_samples : integer
+       The same as the min_samples given to OPTICS. Up and down steep
+       regions can't have more then `min_samples` consecutive non-steep
+       points.
+
+    min_cluster_size : int > 1 or float between 0 and 1
+        Minimum number of samples in an OPTICS cluster, expressed as an
+        absolute number or a fraction of the number of samples (rounded
+        to be at least 2).
+
+    predecessor_correction : bool
+        Correct clusters based on the calculated predecessors.
+
+    Returns
+    -------
+    clusters : array, shape (n_clusters, 2)
+        The list of clusters in the form of [start, end] in each row, with all
+        indices inclusive. The clusters are ordered in a way that larger
+        clusters encompassing smaller clusters come after those smaller
+        clusters.
+    """
+
+    # all indices are inclusive (specially at the end)
+    n_samples = len(reachability_plot)
+    # add an inf to the end of reachability plot
+    # but process the data only until the last actual point
+    # this is fine since the last point is considered upward anyway
+    reachability_plot = np.hstack((reachability_plot, np.inf))
+
+    if min_cluster_size <= 1:
+        min_cluster_size = max(2, min_cluster_size * n_samples)
+
+    xi_complement = 1 - xi
+    sdas = []  # steep down areas, introduced in section 4.3.2 of the paper
+    clusters = []
+    index = 0
+    mib = 0.  # maximum in between
+
+    ratio = reachability_plot[:-1] / reachability_plot[1:]
+    steep_upward = ratio <= xi_complement
+    steep_downward = ratio >= 1 / xi_complement
+    downward = ratio > 1
+    upward = ratio < 1
+
+    # the following loop is is almost exactly as Figure 19 of the paper.
+    # it jumps over the areas which are not either steep down or up areas
+    for steep_index in iter(np.flatnonzero(steep_upward | steep_downward)):
+        # just continue if steep_index has been a part of a discovered xward
+        # area.
+        if steep_index < index:
+            continue
+
+        mib = max(mib, np.max(reachability_plot[index:steep_index + 1]))
+        index = steep_index
+
+        # check if a steep downward area starts
+        if steep_downward[index]:
+            sdas = _update_filter_sdas(sdas, mib, xi_complement)
+            D_start = index
+            index, end = _extend_downward_region(steep_downward, upward,
+                                                 D_start, min_samples)
+            D = _Area(start=D_start, end=end,
+                      maximum=reachability_plot[D_start], mib=0.)
+            sdas.append(D)
+            mib = reachability_plot[end + 1]
+
+        elif steep_upward[index]:
+            sdas = _update_filter_sdas(sdas, mib, xi_complement)
+            U_start = index
+            index, end = _extend_upward_region(steep_upward, downward, U_start,
+                                               min_samples)
+            U = _Area(start=U_start, end=end, maximum=reachability_plot[end],
+                      mib=-1)
+            mib = reachability_plot[end + 1]
+
+            U_clusters = []
+            for D in sdas:
+                c_start = D.start
+                c_end = min(U.end, n_samples - 1)
+
+                # line (**), sc2*
+                if reachability_plot[c_end + 1] * xi_complement < D.mib:
+                    continue
+
+                # Definition 11: criterion 4
+                if D.maximum * xi_complement >= reachability_plot[c_end + 1]:
+                    # Find the first index from the left side which is almost
+                    # at the same level as the end of the detected cluster.
+                    while (reachability_plot[c_start + 1] >
+                           reachability_plot[c_end + 1]
+                           and c_start < c_end):
+                        c_start += 1
+                elif reachability_plot[c_end] * xi_complement >= D.maximum:
path: sklearn/cluster/optics_.py
text: @user1:
``reachability_plot[c_end]`` -> ``reachability_plot[c_end + 1]``?

diff_hunk: @@ -527,3 +568,382 @@ def cluster_optics_dbscan(reachability, core_distances, ordering, eps=0.5):
     labels[ordering] = np.cumsum(far_reach[ordering] & near_core[ordering]) - 1
     labels[far_reach & ~near_core] = -1
     return labels
+
+
+def cluster_optics_xi(reachability, predecessor, ordering, min_samples,
+                      min_cluster_size=None, xi=0.05,
+                      predecessor_correction=False):
+    """Automatically extract clusters according to the Xi-steep method.
+
+    Parameters
+    ----------
+    reachability : array, shape (n_samples,)
+        Reachability distances calculated by OPTICS (`reachability_`)
+
+    predecessor : array, shape (n_samples,)
+        Predecessors calculated by OPTICS.
+
+    ordering : array, shape (n_samples,)
+        OPTICS ordered point indices (`ordering_`)
+
+    min_samples : integer
+       The same as the min_samples given to OPTICS. Up and down steep
+       regions can't have more then `min_samples` consecutive non-steep
+       points.
+
+    min_cluster_size : int > 1 or float between 0 and 1 (default=None)
+        Minimum number of samples in an OPTICS cluster, expressed as an
+        absolute number or a fraction of the number of samples (rounded to be
+        at least 2). If ``None``, the value of ``min_samples`` is used instead.
+
+    xi : float, between 0 and 1, optional (default=0.05)
+        Determines the minimum steepness on the reachability plot that
+        constitutes a cluster boundary. For example, an upwards point in the
+        reachability plot is defined by the ratio from one point to its
+        successor being at most 1-xi.
+
+    predecessor_correction : bool, optional (default=True)
+        Correct clusters based on the calculated predecessors.
+
+    Returns
+    -------
+    labels : array, shape (n_samples)
+        The labels assigned to samples. Points which are not included
+        in any cluster are labeled as -1.
+
+    clusters : array, shape (n_clusters, 2)
+        The list of clusters in the form of ``[start, end]`` in each row, with
+        all indices inclusive. The clusters are ordered according to ``(end,
+        -start)`` (ascending) so that larger clusters encompassing smaller
+        clusters come after such nested smaller clusters. Since ``labels`` does
+        not reflect the hierarchy, usually ``len(clusters) >
+        np.unique(labels)``.
+    """
+    if min_cluster_size is None:
+        min_cluster_size = min_samples
+
+    clusters = _xi_cluster(reachability[ordering], predecessor[ordering], xi,
+                           min_samples, min_cluster_size,
+                           predecessor_correction)
+    labels = _extract_xi_labels(ordering, clusters)
+    return labels, clusters
+
+
+class _Area:
+    """An (upward or downward) area.
+
+    Attributes
+    ----------
+    start : integer
+        The start of the region.
+
+    end : integer
+        The end of the region.
+
+    maximum : float
+        The maximum reachability in this region, which is the
+        start of the region for a downward area and the end of
+        the region for an upward area.
+
+    mib : float
+        Maximum in between value, i.e. the maximum value between
+        the end of a steep down area and the current index. It is
+        irrelevant for steep up areas.
+    """
+    def __init__(self, start, end, maximum, mib):
+        self.start = start
+        self.end = end
+        # maximum is the beginning of a downward area
+        # and the end of an upward area
+        self.maximum = maximum
path: sklearn/cluster/optics_.py
text: @user1:
``self.maximum``  seems redundant? (not mentioned in the paper, we can simply rely on ``self.start``?)

diff_hunk: @@ -527,3 +568,382 @@ def cluster_optics_dbscan(reachability, core_distances, ordering, eps=0.5):
     labels[ordering] = np.cumsum(far_reach[ordering] & near_core[ordering]) - 1
     labels[far_reach & ~near_core] = -1
     return labels
+
+
+def cluster_optics_xi(reachability, predecessor, ordering, min_samples,
+                      min_cluster_size=None, xi=0.05,
+                      predecessor_correction=False):
+    """Automatically extract clusters according to the Xi-steep method.
+
+    Parameters
+    ----------
+    reachability : array, shape (n_samples,)
+        Reachability distances calculated by OPTICS (`reachability_`)
+
+    predecessor : array, shape (n_samples,)
+        Predecessors calculated by OPTICS.
+
+    ordering : array, shape (n_samples,)
+        OPTICS ordered point indices (`ordering_`)
+
+    min_samples : integer
+       The same as the min_samples given to OPTICS. Up and down steep
+       regions can't have more then `min_samples` consecutive non-steep
+       points.
+
+    min_cluster_size : int > 1 or float between 0 and 1 (default=None)
+        Minimum number of samples in an OPTICS cluster, expressed as an
+        absolute number or a fraction of the number of samples (rounded to be
+        at least 2). If ``None``, the value of ``min_samples`` is used instead.
+
+    xi : float, between 0 and 1, optional (default=0.05)
+        Determines the minimum steepness on the reachability plot that
+        constitutes a cluster boundary. For example, an upwards point in the
+        reachability plot is defined by the ratio from one point to its
+        successor being at most 1-xi.
+
+    predecessor_correction : bool, optional (default=True)
+        Correct clusters based on the calculated predecessors.
+
+    Returns
+    -------
+    labels : array, shape (n_samples)
+        The labels assigned to samples. Points which are not included
+        in any cluster are labeled as -1.
+
+    clusters : array, shape (n_clusters, 2)
+        The list of clusters in the form of ``[start, end]`` in each row, with
+        all indices inclusive. The clusters are ordered according to ``(end,
+        -start)`` (ascending) so that larger clusters encompassing smaller
+        clusters come after such nested smaller clusters. Since ``labels`` does
+        not reflect the hierarchy, usually ``len(clusters) >
+        np.unique(labels)``.
+    """
+    if min_cluster_size is None:
+        min_cluster_size = min_samples
+
+    clusters = _xi_cluster(reachability[ordering], predecessor[ordering], xi,
+                           min_samples, min_cluster_size,
+                           predecessor_correction)
+    labels = _extract_xi_labels(ordering, clusters)
+    return labels, clusters
+
+
+class _Area:
+    """An (upward or downward) area.
+
+    Attributes
+    ----------
+    start : integer
+        The start of the region.
+
+    end : integer
+        The end of the region.
+
+    maximum : float
+        The maximum reachability in this region, which is the
+        start of the region for a downward area and the end of
+        the region for an upward area.
+
+    mib : float
+        Maximum in between value, i.e. the maximum value between
+        the end of a steep down area and the current index. It is
+        irrelevant for steep up areas.
+    """
+    def __init__(self, start, end, maximum, mib):
+        self.start = start
+        self.end = end
+        # maximum is the beginning of a downward area
+        # and the end of an upward area
+        self.maximum = maximum
+        # maximum in between (ref Sec. 4.3.2)
+        self.mib = mib
+
+    def __repr__(self):
+        return ("start: %s, end: %s, max: %4g, mib: %4g" %
+                (self.start, self.end, self.maximum, self.mib))
+
+
+def _extend_region(steep_point, xward_point, start, min_samples):
+    """Extend the area until it's maximal.
+
+    It's the same function for both upward and downward reagions, depending on
+    the given input parameters. To extend an upward reagion,
+    ``steep_point=steep_upward`` and ``xward_point=downward`` are expected, and
+    to extend a downward region, ``steep_point=steep_downward`` and
+    ``xward_point=upward``.
+
+    Parameters
+    ----------
+    steep_point : bool array, shape (n_samples)
+        True if the point is steep downward (upward).
+
+    xward_point : bool array, shape (n_samples)
+        True if the point is an upward (respectively downward) point.
+
+    start : integer
+        The start of the xward region.
+
+    min_samples : integer
+       The same as the min_samples given to OPTICS. Up and down steep
+       regions can't have more then ``min_samples`` consecutive non-steep
+       points.
+
+    Returns
+    -------
+    index : integer
+        The current index iterating over all the samples, i.e. where we are up
+        to in our search.
+
+    end : integer
+        The end of the region, which can be behind the index. The region
+        includes the ``end`` index.
+    """
+    n_samples = len(steep_point)
+    non_xward_points = 0
+    index = start
+    end = start
+    # find a maximal area
+    while index < n_samples:
+        if steep_point[index]:
+            non_xward_points = 0
+            end = index
+        elif not xward_point[index]:
+            # it's not a steep point, but still goes up.
+            non_xward_points += 1
+            # region should include no more than min_samples consecutive
+            # non steep xward points.
+            if non_xward_points > min_samples:
+                break
+        else:
+            return min(index, n_samples - 1), end
+
+        index += 1
+    return min(index + 1, n_samples - 1), end
+
+
+def _extend_upward_region(steep_upward_mask, downward_mask, start,
+                          min_samples):
+    return _extend_region(steep_upward_mask, downward_mask,
+                          start, min_samples)
+
+
+def _extend_downward_region(steep_downward_mask, upward_mask, start,
+                            min_samples):
+    return _extend_region(steep_downward_mask, upward_mask,
+                          start, min_samples)
+
+
+def _update_filter_sdas(sdas, mib, xi_complement):
+    """Update steep down areas (SDAs) using the new maximum in between (mib)
+    value, and the given complement of xi, i.e. ``1 - xi``.
+    """
+    if np.isinf(mib):
+        return []
+    res = [sda for sda in sdas if mib <= sda.maximum * xi_complement]
+    for sda in res:
+        sda.mib = max(sda.mib, mib)
+    return res
+
+
+def _correct_predecessor(reachability_plot, predecessor, s, e):
+    """Correct for predecessors.
+
+    Applies Algorithm 2 of [1]_.
+
+    Input parameters are ordered by the computer OPTICS ordering.
+
+    .. [1] Schubert, Erich, Michael Gertz.
+       "Improving the Cluster Structure Extracted from OPTICS Plots." Proc. of
+       the Conference "Lernen, Wissen, Daten, Analysen" (LWDA) (2018): 318-329.
+    """
+
+    while s < e:
+        if reachability_plot[s] > reachability_plot[e]:
+            return s, e
+        for i in range(s, e):
+            if predecessor[e] == i:
+                return s, e
+        e -= 1
+    return None, None
+
+
+def _xi_cluster(reachability_plot, predecessor, xi, min_samples,
+                min_cluster_size, predecessor_correction):
+    """Automatically extract clusters according to the Xi-steep method.
+
+    This is rouphly an implementation of Figure 19 of the OPTICS paper.
+
+    Parameters
+    ----------
+    reachability_plot : array, shape (n_samples)
+        The reachability plot, i.e. reachability ordered according to
+        the calculated ordering, all computed by OPTICS.
+
+    predecessor : array, shape (n_samples)
+        Predecessors ordered according to the calculated ordering.
+
+    xi : float, between 0 and 1
+        Determines the minimum steepness on the reachability plot that
+        constitutes a cluster boundary. For example, an upwards point in the
+        reachability plot is defined by the ratio from one point to its
+        successor being at most 1-xi.
+
+    min_samples : integer
+       The same as the min_samples given to OPTICS. Up and down steep
+       regions can't have more then `min_samples` consecutive non-steep
+       points.
+
+    min_cluster_size : int > 1 or float between 0 and 1
+        Minimum number of samples in an OPTICS cluster, expressed as an
+        absolute number or a fraction of the number of samples (rounded
+        to be at least 2).
+
+    predecessor_correction : bool
+        Correct clusters based on the calculated predecessors.
+
+    Returns
+    -------
+    clusters : array, shape (n_clusters, 2)
+        The list of clusters in the form of [start, end] in each row, with all
+        indices inclusive. The clusters are ordered in a way that larger
+        clusters encompassing smaller clusters come after those smaller
+        clusters.
+    """
+
+    # all indices are inclusive (specially at the end)
+    n_samples = len(reachability_plot)
+    # add an inf to the end of reachability plot
+    # but process the data only until the last actual point
+    # this is fine since the last point is considered upward anyway
+    reachability_plot = np.hstack((reachability_plot, np.inf))
+
+    if min_cluster_size <= 1:
+        min_cluster_size = max(2, min_cluster_size * n_samples)
+
+    xi_complement = 1 - xi
+    sdas = []  # steep down areas, introduced in section 4.3.2 of the paper
+    clusters = []
+    index = 0
+    mib = 0.  # maximum in between
+
+    ratio = reachability_plot[:-1] / reachability_plot[1:]
+    steep_upward = ratio <= xi_complement
+    steep_downward = ratio >= 1 / xi_complement
+    downward = ratio > 1
+    upward = ratio < 1
+
+    # the following loop is is almost exactly as Figure 19 of the paper.
+    # it jumps over the areas which are not either steep down or up areas
+    for steep_index in iter(np.flatnonzero(steep_upward | steep_downward)):
+        # just continue if steep_index has been a part of a discovered xward
+        # area.
+        if steep_index < index:
+            continue
+
+        mib = max(mib, np.max(reachability_plot[index:steep_index + 1]))
+        index = steep_index
+
+        # check if a steep downward area starts
+        if steep_downward[index]:
+            sdas = _update_filter_sdas(sdas, mib, xi_complement)
+            D_start = index
+            index, end = _extend_downward_region(steep_downward, upward,
+                                                 D_start, min_samples)
+            D = _Area(start=D_start, end=end,
+                      maximum=reachability_plot[D_start], mib=0.)
+            sdas.append(D)
+            mib = reachability_plot[end + 1]
+
+        elif steep_upward[index]:
+            sdas = _update_filter_sdas(sdas, mib, xi_complement)
path: sklearn/cluster/optics_.py
text: @user1:
so you don't consider the points in the current steep up region (except for the first point)? (i.e., you only check sc1* for ``sD < x <= sU``)?

@author:
according to Figure 19 of the paper, sdas are updated before changing the `index, and that's why I've implemented like this.