knn predict unreasonably slow b/c of use of scipy.stats.mode
```python
import numpy as np
from sklearn.datasets import make_blobs
from sklearn.neighbors import KNeighborsClassifier

X, y = make_blobs(centers=2, random_state=4, n_samples=30)
knn = KNeighborsClassifier(algorithm='kd_tree').fit(X, y)

x_min, x_max = X[:, 0].min(), X[:, 0].max()
y_min, y_max = X[:, 1].min(), X[:, 1].max()

xx = np.linspace(x_min, x_max, 1000)
# change 100 to 1000 below and wait a long time                                          
yy = np.linspace(y_min, y_max, 100)                                          

X1, X2 = np.meshgrid(xx, yy)                                                  
X_grid = np.c_[X1.ravel(), X2.ravel()]                                        
decision_values = knn.predict(X_grid)
```

spends all it's time in unique within stats.mode, not within the distance calculation. ``mode`` runs ``unique`` for every row.
I'm pretty sure we can replace the call to mode by some call to making a csr matrix and then argmax.

How much is it worth optimizing this? I feel KNN should be fast in low dimensions and people might actually use this. Having the bottleneck in the wrong place just feels wrong to me ;)
WIP PERF Faster KNeighborsClassifier.predict
Closes https://github.com/scikit-learn/scikit-learn/issues/13783

This makes `KNeighborsClassifier.predict` faster by re-writing `scipy.stats.mode` as an argmax of a sparse array as discussed in the parent issue.

Using the example provided in https://github.com/scikit-learn/scikit-learn/issues/13783#issue-440225690,

**On master**
```
%timeit knn.predict(X_grid)                                                 
2.47 s ± 37.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
```
**With this PR**
```
%timeit knn.predict(X_grid)                                                 
793 ms ± 39.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
```

so in this particular case, the `KNeighborsClassifier.predict` is 3.1x faster.

It works in a straightforward way both on weighted an unweighted data making `sklearn.utils.weighted_mode` no longer necessary. 

The downside is that it makes `RadiusNeighborsClassifier.predict` slower by about 30% on the following example,
<details>

```
In [1]: import numpy as np 
   ...: from sklearn.datasets import make_blobs 
   ...: from sklearn.neighbors import RadiusNeighborsClassifier 
   ...:  
   ...: X, y = make_blobs(centers=2, random_state=4, n_samples=30) 
   ...: knn = RadiusNeighborsClassifier(algorithm='kd_tree', radius=2).fit(X, y) 
   ...:  
   ...: x_min, x_max = X[:, 0].min(), X[:, 0].max() 
   ...: y_min, y_max = X[:, 1].min(), X[:, 1].max() 
   ...:  
   ...: xx = np.linspace(x_min, x_max, 100) 
   ...: # change 100 to 1000 below and wait a long time                             
   ...:               
   ...: yy = np.linspace(y_min, y_max, 100)                                         
   ...:   
   ...:  
   ...: X1, X2 = np.meshgrid(xx, yy)                                                
   ...:    
   ...: X_grid = np.c_[X1.ravel(), X2.ravel()]                                      

In [2]: %timeit knn.predict(X_grid)                                                 
1.27 s ± 9.02 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
```
</details>
The difference is mostly that in `RadiusNeighborsClassifier` we call this function repeatedly on 1D arrays, as opposed to a single 2D array.  In worst case, I could revert back to `scipy.stats.mode` and `sklearn.utils.weighted_mode` for `RadiusNeighborsClassifier` but it's annoying to use two different implementations.

**TODO**
 - [ ] fix the remaining 2 test failures on the comparison between the multi-output and single-output case.
 - [ ] ~investigate RadiusNeighborsClassifier.predict performance regression.~