Verified metric support for Hamerly

Author: PyDataBlog

README.md

@@ -66,6 +66,7 @@ ________________________________________________________________________________
 - Support for multi-theading implementation of K-Means clustering algorithm.
 - Kmeans++ initialization for faster and better convergence.
 - Implementation of all the variants of the K-Means algorithm.
+- Supported interface as an [MLJ](https://github.com/alan-turing-institute/MLJ.jl#available-models) model.
 
 _________________________________________________________________________________________________________
 

docs/src/index.md

@@ -51,13 +51,15 @@ pkg> add ParallelKMeans
 The few (and selected) brave ones can simply grab the current experimental features by simply adding the experimental branch to your development environment after invoking the package manager with `]`:
 
 ```julia
-dev git@github.com:PyDataBlog/ParallelKMeans.jl.git
+pkg> add ParallelKMeans#experimental
 ```
 
-Don't forget to checkout the experimental branch and you are good to go with bleeding edge features and breakages!
+You are good to go with bleeding edge features and breakages!
 
-```bash
-git checkout experimental
+To revert to a stable version, you can simply run:
+
+```julia
+pkg> free ParallelKMeans
 ```
 
 ## Features
@@ -78,6 +80,7 @@ git checkout experimental
 - [X] Support of MLJ Random generation hyperparameter.
 - [ ] Support for other distance metrics supported by [Distances.jl](https://github.com/JuliaStats/Distances.jl#supported-distances).
 - [ ] Implementation of [Geometric methods to accelerate k-means algorithm](http://cs.baylor.edu/~hamerly/papers/sdm2016_rysavy_hamerly.pdf).
+- [ ] Support of MLJ Random generation hyperparameter.
 - [ ] Native support for tabular data inputs outside of MLJModels' interface.
 - [ ] Refactoring and finalization of API design.
 - [ ] GPU support.

src/elkan.jl

@@ -266,7 +266,7 @@ function chunk_update_bounds(alg, containers, centroids, r, idx)
     stale = containers.stale
     labels = containers.labels
     T = eltype(centroids)
-
+    # TODO: Add metric support with multiple dispatch
     @inbounds for i in r
         for j in axes(centroids, 2)
             lb[j, i] = lb[j, i] > p[j] ? lb[j, i] + p[j] - T(2)*sqrt(abs(lb[j, i]*p[j])) : zero(T)

src/hamerly.jl

@@ -43,10 +43,10 @@ function kmeans!(alg::Hamerly, containers, X, k, weights;
         collect_containers(alg, containers, n_threads)
 
         J = sum(containers.ub)
-        move_centers(alg, containers, centroids)
+        move_centers(alg, containers, centroids, metric)
 
         r1, r2, pr1, pr2 = double_argmax(p)
-        @parallelize n_threads ncol chunk_update_bounds(alg, containers, r1, r2, pr1, pr2)
+        @parallelize n_threads ncol chunk_update_bounds(alg, containers, r1, r2, pr1, pr2, metric)
 
         if verbose
             # Show progress and terminate if J stops decreasing as specified by the tolerance level.
@@ -241,15 +241,14 @@ end
 Calculates new positions of centers and distance they have moved. Results are stored
 in `centroids` and `p` respectively.
 """
-function move_centers(::Hamerly, containers, centroids)
+function move_centers(::Hamerly, containers, centroids, metric)
     centroids_new = containers.centroids_new[end]
     p = containers.p
     T = eltype(centroids)
 
     @inbounds for i in axes(centroids, 2)
-        d = zero(T)
+        d = distance(metric, centroids, centroids_new, i, i)
         for j in axes(centroids, 1)
-            d += (centroids[j, i] - centroids_new[j, i])^2
             centroids[j, i] = centroids_new[j, i]
         end
         p[i] = d
@@ -258,11 +257,12 @@ end
 
 
 """
-    chunk_update_bounds(alg::Hamerly, containers, r1, r2, pr1, pr2, r, idx)
+    chunk_update_bounds(alg::Hamerly, containers, r1, r2, pr1, pr2, metric::Euclidean, r, idx)
 
-Updates upper and lower bounds of point distance to the centers, with regard to the centers movement.
+Updates upper and lower bounds of point distance to the centers, with regard to the centers movement
+when metric is Euclidean.
 """
-function chunk_update_bounds(alg::Hamerly, containers, r1, r2, pr1, pr2, r, idx)
+function chunk_update_bounds(alg::Hamerly, containers, r1, r2, pr1, pr2, metric::Euclidean, r, idx)
     p = containers.p
     ub = containers.ub
     lb = containers.lb
@@ -297,6 +297,28 @@ function chunk_update_bounds(alg::Hamerly, containers, r1, r2, pr1, pr2, r, idx)
 end
 
 
+"""
+    chunk_update_bounds(alg::Hamerly, containers, r1, r2, pr1, pr2, metric::Metric, r, idx)
+
+Updates upper and lower bounds of point distance to the centers, with regard to the centers movement
+when metric is Euclidean.
+"""
+function chunk_update_bounds(alg::Hamerly, containers, r1, r2, pr1, pr2, metric::Metric, r, idx)
+    p = containers.p
+    ub = containers.ub
+    lb = containers.lb
+    labels = containers.labels
+    T = eltype(containers.ub)
+    # Using notation from original paper, `u` is upper bound and `a` is `labels`, so
+    # `u[i] -> u[i] + p[a[i]]`
+    @inbounds for i in r
+        label = labels[i]
+        ub[i] += p[label]
+        lb[i] -= r1 == label ? pr2 : pr1
+    end
+end
+
+
 """
     double_argmax(p)
 

src/lloyd.jl

@@ -46,7 +46,7 @@ function kmeans!(alg::Lloyd, containers, X, k, weights;
         end
 
         J_previous = J
-        niters += 1
+        niters += 1  # TODO: Investigate the potential bug in number of iterations
     end
     @parallelize n_threads ncol sum_of_squares(containers, X, containers.labels, centroids, weights, metric)
     totalcost = sum(containers.sum_of_squares)

test/test05_hamerly.jl

@@ -118,9 +118,15 @@ end
     Random.seed!(2020)
     X = rand(3, 100)
 
-    res = kmeans(Hamerly(), X, 2, tol = 1e-16, metric=Cityblock())
-    @test res.totalcost ≈ 62.04045252895372
-    @test res.converged
+    baseline = kmeans(Lloyd(), X, 2, tol = 1e-16, metric=Cityblock())
+
+    Random.seed!(2020)
+    X = rand(3, 100)
+
+    res = kmeans(Hamerly(), X, 2; tol = 1e-16, metric=Cityblock())
+    @test res.totalcost ≈ baseline.totalcost
+    @test res.converged == baseline.converged
+    @test res.iterations == baseline.iterations
 end
 
 end # module