Add rng parameter to appropriate methods (#228)

* add `Random` as dependency

* added `rng` parameter to `fuzzy_cmeans`

* added `rng` parameter to `kmeans`

* switch to StableRNGs in `hclust` tests

* use old kwargs syntax

* make k-means result comparison more robust

* generated new test data for `affprop` test

Author: wildart

Project.toml

@@ -10,6 +10,7 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [compat]
 Distances = "0.8, 0.9, 0.10"

src/Clustering.jl

@@ -7,6 +7,7 @@ module Clustering
     using LinearAlgebra
     using SparseArrays
     using Statistics
+    using Random
 
     import Base: show
     import StatsBase: IntegerVector, RealVector, RealMatrix, counts

src/fuzzycmeans.jl

@@ -19,7 +19,7 @@ The output of [`fuzzy_cmeans`](@ref) function.
 struct FuzzyCMeansResult{T<:AbstractFloat}
     centers::Matrix{T}          # cluster centers (d x C)
     weights::Matrix{Float64}    # assigned weights (n x C)
-    iterations::Int             # number of elasped iterations
+    iterations::Int             # number of elapsed iterations
     converged::Bool             # whether the procedure converged
 end
 
@@ -88,14 +88,15 @@ function fuzzy_cmeans(
     maxiter::Int = _fcmeans_default_maxiter,
     tol::Real = _fcmeans_default_tol,
     dist_metric::Metric = Euclidean(),
-    display::Symbol = _fcmeans_default_display
+    display::Symbol = _fcmeans_default_display,
+    rng::AbstractRNG = Random.GLOBAL_RNG
     ) where T<:Real
 
     nrows, ncols = size(data)
     2 <= C < ncols || throw(ArgumentError("C must have 2 <= C < n=$ncols ($C given)"))
     1 < fuzziness || throw(ArgumentError("fuzziness must be greater than 1 ($fuzziness given)"))
 
-    _fuzzy_cmeans(data, C, fuzziness, maxiter, tol, dist_metric, display_level(display))
+    _fuzzy_cmeans(data, C, fuzziness, maxiter, tol, dist_metric, display_level(display),rng)
 
 end
 
@@ -108,13 +109,14 @@ function _fuzzy_cmeans(
     maxiter::Int,                                   # maximum number of iterations
     tol::Real,                                      # tolerance
     dist_metric::Metric,                            # metric to calculate distance
-    displevel::Int                                  # the level of display
+    displevel::Int,                                 # the level of display
+    rng::AbstractRNG                                # RNG object
     ) where T<:Real
 
     nrows, ncols = size(data)
 
     # Initialize weights randomly
-    weights = rand(Float64, ncols, C)
+    weights = rand(rng, Float64, ncols, C)
     weights ./= sum(weights, dims=2)
 
     centers = zeros(T, (nrows, C))

src/kmeans.jl

@@ -49,7 +49,8 @@ function kmeans!(X::AbstractMatrix{<:Real},                # in: data matrix (d
                  maxiter::Integer=_kmeans_default_maxiter, # in: maximum number of iterations
                  tol::Real=_kmeans_default_tol,            # in: tolerance of change at convergence
                  display::Symbol=_kmeans_default_display,  # in: level of display
-                 distance::SemiMetric=SqEuclidean())       # in: function to compute distances
+                 distance::SemiMetric=SqEuclidean(),       # in: function to compute distances
+                 rng::AbstractRNG=Random.GLOBAL_RNG)       # in: RNG object
     d, n = size(X)
     dc, k = size(centers)
     WC = (weights === nothing) ? Int : eltype(weights)
@@ -68,7 +69,7 @@ function kmeans!(X::AbstractMatrix{<:Real},                # in: data matrix (d
         mean!(centers, X)
       end
       return _kmeans!(X, weights, centers, Int(maxiter), Float64(tol),
-                      display_level(display), distance)
+                      display_level(display), distance, rng)
     end
 end
 
@@ -99,19 +100,20 @@ function kmeans(X::AbstractMatrix{<:Real},                # in: data matrix (d x
                 maxiter::Integer=_kmeans_default_maxiter, # in: maximum number of iterations
                 tol::Real=_kmeans_default_tol,            # in: tolerance  of change at convergence
                 display::Symbol=_kmeans_default_display,  # in: level of display
-                distance::SemiMetric=SqEuclidean())       # in: function to calculate distance with
+                distance::SemiMetric=SqEuclidean(),       # in: function to calculate distance with
+                rng::AbstractRNG=Random.GLOBAL_RNG)       # in: RNG object
     d, n = size(X)
     (1 <= k <= n) || throw(ArgumentError("k must be from 1:n (n=$n), k=$k given."))
 
     # initialize the centers using a type wide enough so that the updates
     # centers[i, cj] += X[i, j] * wj will occur without loss of precision through rounding
     T = float(weights === nothing ? eltype(X) : promote_type(eltype(X), eltype(weights)))
-    iseeds = initseeds(init, X, k)
+    iseeds = initseeds(init, X, k, rng=rng)
     centers = copyseeds!(Matrix{T}(undef, d, k), X, iseeds)
 
     kmeans!(X, centers;
             weights=weights, maxiter=Int(maxiter), tol=Float64(tol),
-            display=display, distance=distance)
+            display=display, distance=distance, rng=rng)
 end
 
 #### Core implementation
@@ -123,7 +125,8 @@ function _kmeans!(X::AbstractMatrix{<:Real},                # in: data matrix (d
                   maxiter::Int,                             # in: maximum number of iterations
                   tol::Float64,                             # in: tolerance of change at convergence
                   displevel::Int,                           # in: the level of display
-                  distance::SemiMetric)                     # in: function to calculate distance
+                  distance::SemiMetric,                     # in: function to calculate distance
+                  rng::AbstractRNG)                         # in: RNG object
     d, n = size(X)
     k = size(centers, 2)
     to_update = Vector{Bool}(undef, k) # whether a center needs to be updated
@@ -161,7 +164,7 @@ function _kmeans!(X::AbstractMatrix{<:Real},                # in: data matrix (d
         update_centers!(X, weights, assignments, to_update, centers, wcounts)
 
         if !isempty(unused)
-            repick_unused_centers(X, costs, centers, unused, distance)
+            repick_unused_centers(X, costs, centers, unused, distance, rng)
             to_update[unused] .= true
         end
 
@@ -372,14 +375,15 @@ function repick_unused_centers(X::AbstractMatrix{<:Real}, # in: the data matrix
                                costs::Vector{<:Real},     # in: the current assignment costs (n)
                                centers::AbstractMatrix{<:AbstractFloat}, # out: the centers (d x k)
                                unused::Vector{Int},       # in: indices of centers to be updated
-                               distance::SemiMetric)      # in: function to calculate the distance with
+                               distance::SemiMetric,      # in: function to calculate the distance with
+                               rng::AbstractRNG)          # in: RNG object
     # pick new centers using a scheme like kmeans++
     ds = similar(costs)
     tcosts = copy(costs)
     n = size(X, 2)
 
     for i in unused
-        j = wsample(1:n, tcosts)
+        j = wsample(rng, 1:n, tcosts)
         tcosts[j] = 0
         v = view(X, :, j)
         centers[:, i] = v

src/seeding.jl

@@ -5,14 +5,15 @@
 #   Let alg be an instance of such an algorithm, then it should
 #   support the following usage:
 #
-#       initseeds!(iseeds, alg, X)
-#       initseeds_by_costs!(iseeds, alg, costs)
+#       initseeds!(iseeds, alg, X; kwargs...)
+#       initseeds_by_costs!(iseeds, alg, costs; kwargs...)
 #
 #   Here:
 #       - iseeds:   a vector of resultant indexes of the chosen seeds
 #       - alg:      the seeding algorithm instance
 #       - X:        the data matrix, each column being a data point
 #       - costs:    pre-computed pairwise cost matrix.
+#       - kwargs:   additional kw-arguments, i.e. `rng`
 #
 #   This function returns iseeds
 #
@@ -39,8 +40,8 @@ name of the algorithm.
 
 Returns the vector of `k` seed indices.
 """
-initseeds(alg::SeedingAlgorithm, X::AbstractMatrix{<:Real}, k::Integer) =
-    initseeds!(Vector{Int}(undef, k), alg, X)
+initseeds(alg::SeedingAlgorithm, X::AbstractMatrix{<:Real}, k::Integer; kwargs...) =
+    initseeds!(Vector{Int}(undef, k), alg, X; kwargs...)
 
 """
     initseeds_by_costs(alg::Union{SeedingAlgorithm, Symbol},
@@ -54,8 +55,8 @@ between the points as the cost.
 
 Returns the vector of `k` seed indices.
 """
-initseeds_by_costs(alg::SeedingAlgorithm, costs::AbstractMatrix{<:Real}, k::Integer) =
-    initseeds_by_costs!(Vector{Int}(undef, k), alg, costs)
+initseeds_by_costs(alg::SeedingAlgorithm, costs::AbstractMatrix{<:Real}, k::Integer; kwargs...) =
+    initseeds_by_costs!(Vector{Int}(undef, k), alg, costs; kwargs...)
 
 seeding_algorithm(s::Symbol) =
     s == :rand ? RandSeedAlg() :
@@ -71,14 +72,14 @@ end
 check_seeding_args(X::AbstractMatrix, iseeds::AbstractVector) =
     check_seeding_args(size(X, 2), length(iseeds))
 
-initseeds(algname::Symbol, X::AbstractMatrix{<:Real}, k::Integer) =
-    initseeds(seeding_algorithm(algname), X, k)::Vector{Int}
+initseeds(algname::Symbol, X::AbstractMatrix{<:Real}, k::Integer; kwargs...) =
+    initseeds(seeding_algorithm(algname), X, k; kwargs...)::Vector{Int}
 
-initseeds_by_costs(algname::Symbol, costs::AbstractMatrix{<:Real}, k::Integer) =
-    initseeds_by_costs(seeding_algorithm(algname), costs, k)
+initseeds_by_costs(algname::Symbol, costs::AbstractMatrix{<:Real}, k::Integer; kwargs...) =
+    initseeds_by_costs(seeding_algorithm(algname), costs, k; kwargs...)
 
 # use specified vector of seeds
-function initseeds(iseeds::AbstractVector{<:Integer}, X::AbstractMatrix{<:Real}, k::Integer)
+function initseeds(iseeds::AbstractVector{<:Integer}, X::AbstractMatrix{<:Real}, k::Integer; kwargs...)
     length(iseeds) == k ||
         throw(ArgumentError("The length of seeds vector ($(length(iseeds))) differs from the number of seeds requested ($k)"))
     check_seeding_args(X, iseeds)
@@ -90,8 +91,8 @@ function initseeds(iseeds::AbstractVector{<:Integer}, X::AbstractMatrix{<:Real},
     # NOTE no duplicate checks are done, should we?
     convert(Vector{Int}, iseeds)
 end
-initseeds_by_costs(iseeds::AbstractVector{<:Integer}, costs::AbstractMatrix{<:Real}, k::Integer) =
-    initseeds(iseeds, costs, k) # NOTE: passing costs as X, but should be fine since only size(X, 2) is used
+initseeds_by_costs(iseeds::AbstractVector{<:Integer}, costs::AbstractMatrix{<:Real}, k::Integer; kwargs...) =
+    initseeds(iseeds, costs, k; kwargs...) # NOTE: passing costs as X, but should be fine since only size(X, 2) is used
 
 function copyseeds!(S::AbstractMatrix{<:AbstractFloat},
                     X::AbstractMatrix{<:Real},
@@ -119,9 +120,10 @@ struct RandSeedAlg <: SeedingAlgorithm end
 Initialize `iseeds` with the indices of cluster seeds for the `X` data matrix
 using the `alg` seeding algorithm.
 """
-function initseeds!(iseeds::IntegerVector, alg::RandSeedAlg, X::AbstractMatrix{<:Real})
+function initseeds!(iseeds::IntegerVector, alg::RandSeedAlg, X::AbstractMatrix{<:Real};
+                    rng::AbstractRNG=Random.GLOBAL_RNG)
     check_seeding_args(X, iseeds)
-    sample!(1:size(X, 2), iseeds; replace=false)
+    sample!(rng, 1:size(X, 2), iseeds; replace=false)
 end
 
 """
@@ -135,9 +137,9 @@ Here, `costs[i, j]` is the cost of assigning points ``i`` and ``j``
 to the same cluster. One may, for example, use the squared Euclidean distance
 between the points as the cost.
 """
-function initseeds_by_costs!(iseeds::IntegerVector, alg::RandSeedAlg, X::AbstractMatrix{<:Real})
+function initseeds_by_costs!(iseeds::IntegerVector, alg::RandSeedAlg, X::AbstractMatrix{<:Real}; rng::AbstractRNG=Random.GLOBAL_RNG)
     check_seeding_args(X, iseeds)
-    sample!(1:size(X,2), iseeds; replace=false)
+    sample!(rng, 1:size(X,2), iseeds; replace=false)
 end
 
 """
@@ -157,13 +159,14 @@ struct KmppAlg <: SeedingAlgorithm end
 
 function initseeds!(iseeds::IntegerVector, alg::KmppAlg,
                     X::AbstractMatrix{<:Real},
-                    metric::PreMetric = SqEuclidean())
+                    metric::PreMetric = SqEuclidean();
+                    rng::AbstractRNG=Random.GLOBAL_RNG)
     n = size(X, 2)
     k = length(iseeds)
     check_seeding_args(n, k)
 
     # randomly pick the first center
-    p = rand(1:n)
+    p = rand(rng, 1:n)
     iseeds[1] = p
 
     if k > 1
@@ -173,7 +176,7 @@ function initseeds!(iseeds::IntegerVector, alg::KmppAlg,
         # pick remaining (with a chance proportional to mincosts)
         tmpcosts = zeros(n)
         for j = 2:k
-            p = wsample(1:n, mincosts)
+            p = wsample(rng, 1:n, mincosts)
             iseeds[j] = p
 
             # update mincosts
@@ -188,13 +191,14 @@ function initseeds!(iseeds::IntegerVector, alg::KmppAlg,
 end
 
 function initseeds_by_costs!(iseeds::IntegerVector, alg::KmppAlg,
-                             costs::AbstractMatrix{<:Real})
+                             costs::AbstractMatrix{<:Real};
+                             rng::AbstractRNG=Random.GLOBAL_RNG)
     n = size(costs, 1)
     k = length(iseeds)
     check_seeding_args(n, k)
 
     # randomly pick the first center
-    p = rand(1:n)
+    p = rand(rng, 1:n)
     iseeds[1] = p
 
     if k > 1
@@ -203,7 +207,7 @@ function initseeds_by_costs!(iseeds::IntegerVector, alg::KmppAlg,
 
         # pick remaining (with a chance proportional to mincosts)
         for j = 2:k
-            p = wsample(1:n, mincosts)
+            p = wsample(rng, 1:n, mincosts)
             iseeds[j] = p
 
             # update mincosts
@@ -230,7 +234,7 @@ Choose the ``k`` points with the highest *centrality* as seeds.
 struct KmCentralityAlg <: SeedingAlgorithm end
 
 function initseeds_by_costs!(iseeds::IntegerVector, alg::KmCentralityAlg,
-                             costs::AbstractMatrix{<:Real})
+                             costs::AbstractMatrix{<:Real}; kwargs...)
 
     n = size(costs, 1)
     k = length(iseeds)
@@ -255,5 +259,5 @@ function initseeds_by_costs!(iseeds::IntegerVector, alg::KmCentralityAlg,
 end
 
 initseeds!(iseeds::IntegerVector, alg::KmCentralityAlg, X::AbstractMatrix{<:Real},
-           metric::PreMetric = SqEuclidean()) =
-    initseeds_by_costs!(iseeds, alg, pairwise(metric, X, dims=2))
+           metric::PreMetric = SqEuclidean(); kwargs...) =
+    initseeds_by_costs!(iseeds, alg, pairwise(metric, X, dims=2); kwargs...)