Merge pull request #250 from JuliaStats/ast/fix_kmedoids

kmedoids(): fix duplicate medoids case

Author: alyst

src/kmedoids.jl

@@ -182,59 +182,51 @@ end
 
 
 # update assignments and related quantities
-function _kmed_update_assignments!(dist::AbstractMatrix{T},      # in: (n, n)
+# returns the total cost and the number of assignment changes
+function _kmed_update_assignments!(dist::AbstractMatrix{<:Real}, # in: (n, n)
                                    medoids::AbstractVector{Int}, # in: (k,)
                                    assignments::Vector{Int},     # out: (n,)
                                    groups::Vector{Vector{Int}},  # out: (k,)
-                                   costs::Vector{T},             # out: (n,)
-                                   isinit::Bool) where T                 # in
+                                   costs::AbstractVector{<:Real},# out: (n,)
+                                   initial::Bool)                # in
     n = size(dist, 1)
     k = length(medoids)
-    ch = 0
 
-    if !isinit
-        for i = 1:k
-            empty!(groups[i])
-        end
-    end
+    # reset cluster groups (note: assignments are not touched yet)
+    initial || foreach(empty!, groups)
 
     tcost = 0.0
+    ch = 0
     for j = 1:n
-        p = 1
+        p = 1 # initialize the closest medoid for j
         mv = dist[medoids[1], j]
 
-        for i = 2:k
-            v = dist[medoids[i], j]
-            if v < mv
+        # find the closest medoid for j
+        @inbounds for i = 2:k
+            m = medoids[i]
+            v = dist[m, j]
+            # assign if current medoid is closer or if it is j itself
+            if (v < mv) || (m == j)
+                (v <= mv) || throw(ArgumentError("sample #$j reassigned from medoid[$p]=#$(medoids[p]) (distance=$mv) to medoid[$i]=#$m (distance=$v); check the distance matrix correctness"))
                 p = i
                 mv = v
             end
         end
 
-        if isinit
-            assignments[j] = p
-        else
-            a = assignments[j]
-            if p != a
-                ch += 1
-            end
-            assignments[j] = p
-        end
-
+        ch += !initial && (p != assignments[j])
+        assignments[j] = p
         costs[j] = mv
         tcost += mv
         push!(groups[p], j)
     end
 
-    return (tcost, ch)::Tuple{Float64, Int}
+    return (tcost, ch)
 end
 
 
 # find medoid for a given group
-#
-# TODO: faster way without creating temporary arrays
-function _find_medoid(dist::AbstractMatrix, grp::Vector{Int})
+function _find_medoid(dist::AbstractMatrix, grp::AbstractVector{Int})
     @assert !isempty(grp)
-    p = argmin(sum(dist[grp, grp], dims=2))
-    return grp[p]::Int
+    p = argmin(sum(view(dist, grp, grp), dims=2))
+    return grp[p]
 end

test/kmedoids.jl

@@ -9,11 +9,15 @@ include("test_helpers.jl")
     Random.seed!(34568)
     @test_throws ArgumentError kmedoids(randn(2, 3), 1)
     @test_throws ArgumentError kmedoids(randn(2, 3), 4)
-    dist = inv.(max.(pairwise(Euclidean(), randn(2, 3), dims=2), 0.1))
-    @test kmedoids(dist, 2) isa KmedoidsResult
+    dist = max.(pairwise(Euclidean(), randn(2, 3), dims=2), 0.1)
+    @test @inferred(kmedoids(dist, 2)) isa KmedoidsResult
+    # incorrect distance matrix
+    invdist = inv.(max.(pairwise(Euclidean(), randn(2, 3), dims=2), 0.1))
+    @test_throws ArgumentError kmedoids(invdist, 2)
+
     @test_throws ArgumentError kmedoids(dist, 2, display=:mylog)
     for disp in keys(Clustering.DisplayLevels)
-        @test kmedoids(dist, 2, display=disp) isa KmedoidsResult
+        @test @inferred(kmedoids(dist, 2, display=disp)) isa KmedoidsResult
     end
 end
 
@@ -28,7 +32,7 @@ dist = pairwise(SqEuclidean(), X, dims=2)
 @assert size(dist) == (n, n)
 
 Random.seed!(34568)  # reset seed again to known state
-R = kmedoids(dist, k)
+R = @inferred(kmedoids(dist, k))
 @test isa(R, KmedoidsResult)
 @test nclusters(R) == k
 @test length(R.medoids) == length(unique(R.medoids))
@@ -48,36 +52,54 @@ R = kmedoids(dist, k)
     end
 end
 
-# k=1 and k=n cases
-x = pairwise(SqEuclidean(), [1 2 3; .1 .2 .3; 4 5.6 7], dims=2)
-kmed1 = kmedoids(x, 1)
-@test nclusters(kmed1) == 1
-@test assignments(kmed1) == [1, 1, 1]
-@test kmed1.medoids == [2]
-kmed3 = kmedoids(x, 3)
-@test nclusters(kmed3) == 3
-@test sort(assignments(kmed3)) == [1, 2, 3]
-@test sort(kmed3.medoids) == [1, 2, 3]
-
-
-# this data set has three obvious groups:
-# group 1: [1, 3, 4], values: [1, 2, 3]
-# group 2: [2, 5, 7], values: [6, 7, 8]
-# group 3: [6, 8, 9], values: [21, 20, 22]
-#
-
-X = reshape(map(Float64, [1, 6, 2, 3, 7, 21, 8, 20, 22]), 1, 9)
-dist = pairwise(SqEuclidean(), X, dims=2)
+@testset "Duplicated points (#231)" begin
+    pts = [0.0 0.0]
+    dists = pairwise(SqEuclidean(), pts, dims=2)
+    dupmed = kmedoids(dists, 2)
+    @test nclusters(dupmed) == 2
+    @test sort(dupmed.medoids) == [1, 2]
+    @test sort(dupmed.assignments) == [1, 2]
+end
 
-R = kmedoids!(dist, [1, 2, 6])
-@test isa(R, KmedoidsResult)
-@test nclusters(R) == 3
-@test R.medoids == [3, 5, 6]
-@test R.assignments == [1, 2, 1, 1, 2, 3, 2, 3, 3]
-@test counts(R) == [3, 3, 3]
-@test wcounts(R) == counts(R)
-@test R.costs ≈ [1, 1, 0, 1, 0, 0, 1, 1, 1]
-@test R.totalcost ≈ 6.0
-@test R.converged
+@testset "Toy example #1" begin
+    pts = [1 2 3; .1 .2 .3; 4 5.6 7]
+    # k=1 and k=n cases
+    dists = pairwise(SqEuclidean(), pts, dims=2)
+
+    @testset "k=1" begin
+        kmed1 = @inferred(kmedoids(dists, 1))
+        @test nclusters(kmed1) == 1
+        @test assignments(kmed1) == [1, 1, 1]
+        @test kmed1.medoids == [2]
+    end
+
+    @testset "k=3" begin
+        kmed3 = @inferred(kmedoids(dists, 3))
+        @test nclusters(kmed3) == 3
+        @test sort(assignments(kmed3)) == [1, 2, 3]
+        @test sort(kmed3.medoids) == [1, 2, 3]
+    end
+end
+
+@testset "Toy example #2" begin
+    pts = reshape(map(Float64, [1, 6, 2, 3, 7, 21, 8, 20, 22]), 1, 9)
+    # this data set has three obvious groups:
+    # group 1: [1, 3, 4], values: [1, 2, 3]
+    # group 2: [2, 5, 7], values: [6, 7, 8]
+    # group 3: [6, 8, 9], values: [21, 20, 22]
+
+    dists = pairwise(SqEuclidean(), pts, dims=2)
+
+    R = @inferred(kmedoids!(dists, [1, 2, 6]))
+    @test isa(R, KmedoidsResult)
+    @test nclusters(R) == 3
+    @test R.medoids == [3, 5, 6]
+    @test R.assignments == [1, 2, 1, 1, 2, 3, 2, 3, 3]
+    @test counts(R) == [3, 3, 3]
+    @test wcounts(R) == counts(R)
+    @test R.costs ≈ [1, 1, 0, 1, 0, 0, 1, 1, 1]
+    @test R.totalcost ≈ 6.0
+    @test R.converged
+end
 
 end