Merge pull request #193 from oxinabox/ox/chillaxonthematrix

Make clustering methods accept AbstractMatrix

Author: alyst

src/affprop.jl

@@ -33,7 +33,7 @@ const _afp_default_tol = 1.0e-6
 const _afp_default_display = :none
 
 """
-    affinityprop(S::DenseMatrix; [maxiter=200], [tol=1e-6], [damp=0.5],
+    affinityprop(S::AbstractMatrix; [maxiter=200], [tol=1e-6], [damp=0.5],
                  [display=:none]) -> AffinityPropResult
 
 Perform affinity propagation clustering based on a similarity matrix `S`.
@@ -52,7 +52,7 @@ of the ``i``-th point as an *exemplar*.
 > Brendan J. Frey and Delbert Dueck. *Clustering by Passing Messages
 > Between Data Points.* Science, vol 315, pages 972-976, 2007.
 """
-function affinityprop(S::DenseMatrix{T};
+function affinityprop(S::AbstractMatrix{T};
                       maxiter::Integer=_afp_default_maxiter,
                       tol::Real=_afp_default_tol,
                       damp::Real=_afp_default_damp,
@@ -72,7 +72,7 @@ end
 
 #### Implementation
 
-function _affinityprop(S::DenseMatrix{T},
+function _affinityprop(S::AbstractMatrix{T},
                        maxiter::Int,
                        tol::Real,
                        damp::T,
@@ -134,7 +134,7 @@ end
 
 
 # compute responsibilities
-function _afp_compute_r!(R::Matrix{T}, S::DenseMatrix{T}, A::Matrix{T}) where T
+function _afp_compute_r!(R::Matrix{T}, S::AbstractMatrix{T}, A::Matrix{T}) where T
     n = size(S, 1)
 
     I1 = Vector{Int}(undef, n)  # I1[i] is the column index of the maximum element in (A+S)[i,:]
@@ -245,7 +245,7 @@ function _afp_extract_exemplars(A::Matrix, R::Matrix)
 end
 
 # get assignments
-function _afp_get_assignments(S::DenseMatrix, exemplars::Vector{Int})
+function _afp_get_assignments(S::AbstractMatrix, exemplars::Vector{Int})
     n = size(S, 1)
     k = length(exemplars)
     Se = S[:, exemplars]

src/dbscan.jl

@@ -43,7 +43,7 @@ end
 ## main algorithm
 
 """
-    dbscan(D::DenseMatrix, eps::Real, minpts::Int) -> DbscanResult
+    dbscan(D::AbstractMatrix, eps::Real, minpts::Int) -> DbscanResult
 
 Perform DBSCAN algorithm using the distance matrix `D`.
 
@@ -54,7 +54,7 @@ The following options control which points would be considered
   - `minpts::Int`: the minimum number of neighboring points (including itself)
      to qualify a point as a density point.
 """
-function dbscan(D::DenseMatrix{T}, eps::Real, minpts::Int) where T<:Real
+function dbscan(D::AbstractMatrix{T}, eps::Real, minpts::Int) where T<:Real
     # check arguments
     n = size(D, 1)
     size(D, 2) == n || throw(ArgumentError("D must be a square matrix ($(size(D)) given)."))
@@ -66,7 +66,7 @@ function dbscan(D::DenseMatrix{T}, eps::Real, minpts::Int) where T<:Real
     _dbscan(D, convert(T, eps), minpts, 1:n)
 end
 
-function _dbscan(D::DenseMatrix{T}, eps::T, minpts::Int, visitseq::AbstractVector{Int}) where T<:Real
+function _dbscan(D::AbstractMatrix{T}, eps::T, minpts::Int, visitseq::AbstractVector{Int}) where T<:Real
     n = size(D, 1)
 
     # prepare
@@ -96,7 +96,7 @@ end
 
 ## key steps
 
-function _dbs_region_query(D::DenseMatrix{T}, p::Int, eps::T) where T<:Real
+function _dbs_region_query(D::AbstractMatrix{T}, p::Int, eps::T) where T<:Real
     n = size(D,1)
     nbs = Int[]
     dists = view(D,:,p)
@@ -108,7 +108,7 @@ function _dbs_region_query(D::DenseMatrix{T}, p::Int, eps::T) where T<:Real
     return nbs::Vector{Int}
 end
 
-function _dbs_expand_cluster!(D::DenseMatrix{T},           # distance matrix
+function _dbs_expand_cluster!(D::AbstractMatrix{T},        # distance matrix
                               k::Int,                      # the index of current cluster
                               p::Int,                      # the index of seeding point
                               nbs::Vector{Int},            # eps-neighborhood of p

src/kmedoids.jl

@@ -38,7 +38,7 @@ const _kmed_default_tol = 1.0e-8
 const _kmed_default_display = :none
 
 """
-    kmedoids(dist::DenseMatrix, k::Integer; ...) -> KmedoidsResult
+    kmedoids(dist::AbstractMatrix, k::Integer; ...) -> KmedoidsResult
 
 Perform K-medoids clustering of ``n`` points into `k` clusters,
 given the `dist` matrix (``n×n``, `dist[i, j]` is the distance
@@ -59,7 +59,7 @@ The function implements a *K-means style* algorithm instead of *PAM*
 iterations, but was shown to produce worse (10-20% higher total costs) results
 (see e.g. [Schubert & Rousseeuw (2019)](@ref kmedoid_refs)).
 """
-function kmedoids(dist::DenseMatrix{T}, k::Integer;
+function kmedoids(dist::AbstractMatrix{T}, k::Integer;
                   init=_kmed_default_init,
                   maxiter::Integer=_kmed_default_maxiter,
                   tol::Real=_kmed_default_tol,
@@ -79,7 +79,7 @@ function kmedoids(dist::DenseMatrix{T}, k::Integer;
 end
 
 """
-    kmedoids!(dist::DenseMatrix, medoids::Vector{Int};
+    kmedoids!(dist::AbstractMatrix, medoids::Vector{Int};
               [kwargs...]) -> KmedoidsResult
 
 Update the current cluster `medoids` using the `dist` matrix.
@@ -89,7 +89,7 @@ as `medoids` argument.
 
 See [`kmedoids`](@ref) for the description of optional `kwargs`.
 """
-function kmedoids!(dist::DenseMatrix{T}, medoids::Vector{Int};
+function kmedoids!(dist::AbstractMatrix{T}, medoids::Vector{Int};
                    maxiter::Integer=_kmed_default_maxiter,
                    tol::Real=_kmed_default_tol,
                    display::Symbol=_kmed_default_display) where T<:Real
@@ -110,7 +110,7 @@ end
 #### core algorithm
 
 function _kmedoids!(medoids::Vector{Int},      # initialized medoids
-                    dist::DenseMatrix{T},      # distance matrix
+                    dist::AbstractMatrix{T},   # distance matrix
                     maxiter::Int,              # maximum number of iterations
                     tol::Real,                 # tolerable change of objective
                     displevel::Int) where T<:Real            # level of display
@@ -182,7 +182,7 @@ end
 
 
 # update assignments and related quantities
-function _kmed_update_assignments!(dist::DenseMatrix{T},         # in: (n, n)
+function _kmed_update_assignments!(dist::AbstractMatrix{T},      # in: (n, n)
                                    medoids::AbstractVector{Int}, # in: (k,)
                                    assignments::Vector{Int},     # out: (n,)
                                    groups::Vector{Vector{Int}},  # out: (k,)
@@ -233,7 +233,7 @@ end
 # find medoid for a given group
 #
 # TODO: faster way without creating temporary arrays
-function _find_medoid(dist::DenseMatrix, grp::Vector{Int})
+function _find_medoid(dist::AbstractMatrix, grp::Vector{Int})
     @assert !isempty(grp)
     p = argmin(sum(dist[grp, grp], dims=2))
     return grp[p]::Int

test/affprop.jl

@@ -6,6 +6,7 @@ using Clustering
 using LinearAlgebra
 using Random
 using Statistics
+include("test_helpers.jl")
 
 @testset "affinityprop() (affinity propagation)" begin
     @testset "Argument checks" begin
@@ -46,6 +47,13 @@ using Statistics
         @test R.counts[i] == count(==(i), R.assignments)
     end
 
+    @testset "Support for arrays other than Matrix{T}" begin
+        @testset "$(typeof(M))" for M in equivalent_matrices(S)
+            R2 = affinityprop(M)
+            @test R2.assignments == R.assignments
+        end
+    end
+
     #= compare with python result
     the reference assignments were computed using python sklearn:
     ```julia

test/dbscan.jl

@@ -1,6 +1,7 @@
 using Test
 using Clustering
 using Distances
+include("test_helpers.jl")
 
 @testset "dbscan() (DBSCAN clustering)" begin
 
@@ -36,6 +37,13 @@ for c = 1:k
 end
 @test all(R.counts .>= 180)
 
+@testset "Support for arrays other than Matrix{T}" begin
+    @testset "$(typeof(M))" for M in equivalent_matrices(D)
+        R2 = dbscan(M, 1.0, 10)  # run on complete subarray
+        @test R2.assignments == R.assignments
+    end
+end
+
 @testset "normal points" begin
     Random.seed!(0)
     p0 = randn(3, 1000)
@@ -65,5 +73,4 @@ end
     end
 end
 
-
 end

test/kmedoids.jl

@@ -1,7 +1,7 @@
 using Test
-
 using Distances
 using Clustering
+include("test_helpers.jl")
 
 @testset "kmedoids() (k-medoids)" begin
 
@@ -27,6 +27,7 @@ X = rand(d, n)
 dist = pairwise(SqEuclidean(), X, dims=2)
 @assert size(dist) == (n, n)
 
+Random.seed!(34568)  # reset seed again to known state
 R = kmedoids(dist, k)
 @test isa(R, KmedoidsResult)
 @test nclusters(R) == k
@@ -39,6 +40,14 @@ R = kmedoids(dist, k)
 @test isapprox(sum(R.costs), R.totalcost)
 @test R.converged
 
+@testset "Support for arrays other than Matrix{T}" begin
+    @testset "$(typeof(M))" for M in equivalent_matrices(dist)
+        Random.seed!(34568)  # restore seed as kmedoids is not determantistic
+        R2 = kmedoids(M, k)
+        @test R2.assignments == R.assignments
+    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/test_helpers.jl

@@ -0,0 +1,24 @@
+using LinearAlgebra
+using SparseArrays
+
+"""
+    equivalent_matrices(x::AbstractMatrix)
+
+Returns a collection of matrixes that are equal to the input `x`, but of a different type.
+Useful for testing if things still work on different types of matrix.
+"""
+function equivalent_matrices(x::AbstractMatrix)
+    mats = [
+        Base.PermutedDimsArray(x, (1,2)),  # identity permutation
+        view(x, :, :),
+        view(x, collect.(axes(x))...),  # breaks `strides`
+        sparse(x),
+    ]
+    if issymmetric(x)
+        append!(mats, [
+            Symmetric(x),
+            Transpose(x),
+        ])
+    end
+    return mats
+end