I've added some tests, but way many more tests can be written for this class than what we have now.

Also, are you comfortable with the init parameters and their descriptions? I feel like for a good portion of them, the user need to read the code to understand what they do. Would that be a good idea to add an experimental tag to the class to be able to easily change the API in a non-backward-compatible way?

I would like to hope that most or all subsequent changes here are small bug fixes rather than substantial changes to the algorithm. We have tended in the past to not release things as experimental. I know we have done so with a couple of things in pretty new territory in this release, but I'm not enthusiastic about doing it here.

I'll give this another review soon.

There are some edge case issues left. This is the latest I encountered:

>>> from sklearn.cluster import OPTICS
>>> import numpy as np
>>> 
>>> np.random.seed(0)
>>> n_points_per_cluster = 3
>>> 
>>> C1 = [-5, -2] + .8 * np.random.randn(n_points_per_cluster, 2)
>>> C2 = [4, -1] + .1 * np.random.randn(n_points_per_cluster, 2)
>>> X = np.vstack((C1, np.array([[100, 200]]), C2))
>>> #X = np.vstack((C1, ))
>>> 
>>> clust = OPTICS(min_samples=2, max_bound=2).fit(X)
/home/adrin/Projects/github.com/sklearn/.venv/lib/python3.7/site-packages/sklearn/cluster/optics_.py:748: RuntimeWarning: invalid value encountered in double_scalars
  (avg_reach2 / reachability_plot[s]) > maxima_ratio):
>>> clust.labels_
array([0, 0, 0, 1, 1, 1, 1])
>>> 
>>> c = clust
>>> 
>>> print(np.hstack((X[c.ordering_, :], clust.labels_[c.ordering_, np.newaxis])))
[[ -3.58875812  -1.67987423   0.        ]
 [ -3.50595361  -2.7818223    0.        ]
 [ -4.21700961  -0.20728544   0.        ]
 [100.         200.           1.        ]
 [  4.09500884  -1.01513572   1.        ]
 [  3.98967811  -0.95894015   1.        ]
 [  4.01440436  -0.85457265   1.        ]]
>>> 
>>> print(np.hstack((c.reachability_[c.ordering_].reshape(-1, 1),
...                  c.core_distances_[c.ordering_].reshape(-1, 1),
...                  c.labels_[c.ordering_].reshape(-1, 1))))
[[           inf 1.10505481e+00 0.00000000e+00]
 [1.10505481e+00 1.10505481e+00 0.00000000e+00]
 [1.60100521e+00 1.60100521e+00 0.00000000e+00]
 [           inf 2.22611307e+02 1.00000000e+00]
 [           inf 1.19383853e-01 1.00000000e+00]
 [1.19383853e-01 1.07256525e-01 1.00000000e+00]
 [1.07256525e-01 1.07256525e-01 1.00000000e+00]]

I've already fixed the fist line of _extract_optics as:

normalization_factor = np.max(reachability[reachability < np.inf])

But the rest of the code doesn't handle situations where there are more than one inf in the reachability distances.

EDIT: fixed.

I haven't looked into OPTICS at all, but right now there's a test failure on 32 bit on master: https://travis-ci.org/MacPython/scikit-learn-wheels/jobs/404543529#L1208

Any input appreciated.

If I'm not mistaken, @espg has worked on that part.

That test checks the reachability_ values calculated by OPTICS. It's basically a data -> get neighbors -> get core dists -> calculate reachability -> extract clusters pipeline, and this array is the reachability one.

Values are set on this line: https://github.com/scikit-learn/scikit-learn/blob/master/sklearn/cluster/optics_.py#L405

 self.reachability_[unproc] = new_reach

I could spend some time and look into it further if I had access to a 32 bit system, but I don't. Sorry, but that's as much as I know @amueller !

@amueller I don't have access to a 32-bit system either, but I'm curious if the issue is in _optics_inner, which is called by quick_sort:

https://github.com/scikit-learn/scikit-learn/blob/master/sklearn/cluster/optics_.py#L411

That code expects 64-bit floats, mainly because I wasn't sure how to write cython code that handled both data types and didn't want to clutter the main class with conditional checks on the distances. That's why OPTICS requires a float array (if you're using manhattan distances or some other integer based distance, they're converted to floats when the X array is passed in). Numpy defaults to using 64-bit floats as well...but I don't know how numpy behaves on 32-bit machines.

I'd have to look at the output of the array and look at the mismatch to be sure, but my guess is that the numbers in the array maybe aren't different, but are perhaps out of order in one or two places where the reachability distances are "small" -- i.e., you may have numbers that are different in 64-bit and sorted from "b a" to "a b" due to 'a' being marginally smaller then 'b'... but in 32-bit those numbers could be truncated so that 'a' == 'b', and the default behavior in quick_sort is that tied numbers are left in the original order, so "b a" would stay "b a".

If this is in fact what is happening, the data in the code for test could get changed to increase the magnitude of the points reference array. Multiplying by a 100 might be enough to ensure that small distances still maintain magnitude ordering; provided of course that this is the actual cause of the test failure.

The error we see on #11878 isn't from the cython code @espg. I moved that code into python, and it's the same issue. To be clear, almost all values are almost identical to the sixth digit precision, there's just one, which is like this:

0.4293346721668264, 0.427348, diff: 0.0019866721668264087

The difference stays the same when I take the quick_scan out of cython.

P.S. All the tests are on a 32 bit debian 9 VM @amueller :P

P.P.S. still trying to figure out what the issue is.

So there's one point out of 1500 that is off by 0.002, but only in 32 bit mode... that's really odd. Do you know if numpy still represents the numbers as 64-bit when running on a 32-bit machine? I'm wondering if maybe this is an issue in pairwise_distances, although I don't know why there would only be one data point impacted...

@jnothman could you please give me a review on this one?

I'll be working on the 32bit issue, but that can kinda be treated independently, since it's coming from a different part of the code anyway.

After you review, I'll remove/change my notes accordingly. You'll see some questions/comments in this PR in the code, on which I'll appreciate your feed back.

@espg, the issue persists if we multiply all inputs by 100. Then we'll have this as a large change:

73.32993374225046 74.76657291338213 -1.4366391711316737

And the issue seems to be that on a 32 bit system, the points are processed in a different order. And if the data is multiplied by 100, they're still processed in a different order, but differently (compared to when they're not scaled up).

So on this line:

point = self._set_reach_dist(point, X, nbrs)

at some point the algorithm on 32 bit diverges from the path taken on a 64 bit system. My debugging are happening here: https://github.com/adrinjalali/scikit-learn/tree/bug/optics32

So I finally managed to figure out why the OPTICS diverges its path on 32 bit from 64 bit.

The issue is that at some point the reachability distances (new_reach to be specific) start being different (and the difference being of the order of 1e-13). Then, in quick_scan, those values are compared to each other (with == and <) and that 1e-13 makes a difference there.

To test, I added a tolerance to all comparisons in quick_scan as done here, and the final reachability results in 32 bit and 64 bit are now almost identical.

great work @adrinjalali -- thank you for tracking that down!

I'm away for the weekend. Most probably I won't be able to open that PR before Monday :(

@jnothman I kinda need help how to proceed here.

Maybe put separate changes like max_bounds -> max_eps in a separate PR so this remains tractable to review.

+1 for max_bounds -> max_eps. It seems more reasonable. (After the change we use max_eps and eps, R dbscan uses eps and eps_cl)
@adrinjalali Could you please put the rename in a seperate PR? Thanks.

@adrinjalali I've been trying to catch up this thread... did you still have a pending question for me?