Merge pull request #54 from kalmarek/mk/poisson_sampler

Add RandomWordSampler for infinite groups

Project.toml

@@ -1,7 +1,7 @@
 name = "GroupsCore"
 uuid = "d5909c97-4eac-4ecc-a3dc-fdd0858a4120"
 authors = ["Marek Kaluba <[email protected]>"]
-version = "0.5.1"
+version = "0.5.2"
 
 [deps]
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

docs/src/groups.md

@@ -98,3 +98,15 @@ information is derivable from the type, there is no need to extend
 Base.isfinite(G::Group)
 istrivial(G::Group)
 ```
+
+## Random elements
+
+We provide two methods for generating random elements of a group or monoid.
+
+```@docs
+GroupsCore.ProductReplacementSampler
+GroupsCore.RandomWordSampler
+```
+
+By default for finite monoids `ProductReplacementSampler` is used and
+`RandomWordSampler` following `Poisson(λ=8)` is employed for inifinite ones.

src/GroupsCore.jl

@@ -13,8 +13,8 @@ export istrivial
 # export one!, inv!, mul!, conj!, commutator!, div_left!, div_right!
 
 include("exceptions.jl")
-include("groups.jl")
-include("group_elements.jl")
+include("monoids_groups.jl")
+include("elements.jl")
 
 include("rand.jl")
 include("extensions.jl")

src/groups.jl → src/monoids_groups.jl

@@ -52,17 +52,17 @@ function Base.iterate(M::Monoid)
     hasgens(M) && throw(
         InterfaceNotImplemented(:Iteration, "Base.iterate(::$(typeof(M)))"),
     )
-    throw(ArgumentError("Group does not have assigned generators."))
+    throw(ArgumentError("Monoid does not have assigned generators."))
 end
 
-function Base.iterate(M::Group, state)
+function Base.iterate(M::Monoid, state)
     hasgens(M) && throw(
         InterfaceNotImplemented(
             :Iteration,
             "Base.iterate(::$(typeof(M)), state)",
         ),
     )
-    throw(ArgumentError("Group does not have assigned generators."))
+    throw(ArgumentError("Monoid does not have assigned generators."))
 end
 
 """

src/rand.jl

@@ -1,41 +1,64 @@
 import Random
 
+function Random.Sampler(
+    RNG::Type{<:Random.AbstractRNG},
+    M::Monoid,
+    repetition::Random.Repetition = Val(Inf),
+)
+    if isfinite(M)
+        return ProductReplacementSampler(RNG(), M)
+    else
+        # for infinite groups/monoids PRA will either return
+        # ridiculously long words or just run out of memory
+        S = gens(M)
+        if M isa Group
+            S = [S; inv.(S)]
+        end
+        return RandomWordSampler(S, Poisson(; λ = 8))
+    end
+end
+
 """
-    PRASampler
-Product Replacement Algorithm sampler for arbitrary group generated by
+    ProductReplacementSampler
+Implements Product Replacement Algorithm for a group or monoid generated by
 an explicit finite set of generators.
+
+Product Replacement Algorithm performs a random walk on the graph of
+`n`-generating tuples of a group (or monoid) with __two accumulators__.
+Each step consists of
+ 1. multiplying the right accumulator by a random generator,
+ 2. replacing one of the generators with the product with the right accumulator,
+ 3. multiplying the left accumulator (on the left) by a random generator.
+
+The left accumulator is returned as a random element.
+
+By default for a group with `k` generators we set `n = 2k + 10` and perform
+`10*n` scrambling steps before returning a random element.
+
+`PRASampler` provides provably uniformly distributed random elements for
+finite groups. For infinite groups [`RandomWordSampler`](@ref) is used.
+
+!!! warning
+    Using `PRASampler` for an infinite group is ill-advised as the exponential
+    growth of words during scrambling will result in excessive memory use and
+    out-of-memory situation.
 """
-mutable struct PRASampler{T} <: Random.Sampler{T}
+mutable struct ProductReplacementSampler{T} <: Random.Sampler{T}
     gentuple::Vector{T}
     right::T
     left::T
 end
 
 # constants taken from GAP
-function PRASampler(
-    M,
-    n::Integer = 2ngens(M) + 10,
-    scramble_time::Integer = 10max(n, 10),
-)
-    return PRASampler(Random.default_rng(), M, n, scramble_time)
-end
-
-function Random.Sampler(
-    RNG::Type{<:Random.AbstractRNG},
-    M::Monoid,
-    repetition::Random.Repetition = Val(Inf),
-)
-    return PRASampler(RNG(), M)
-end
-
-function PRASampler(
+function ProductReplacementSampler(
     rng::Random.AbstractRNG,
     M::Monoid,
     n::Integer = 2ngens(M) + 10,
     scramble_time::Integer = 10max(n, 10),
 )
+    @assert hasgens(M)
     if istrivial(M)
-        return PRASampler(fill(one(M), n), one(M), one(M))
+        return ProductReplacementSampler(fill(one(M), n), one(M), one(M))
     end
     @assert hasgens(M)
     l = max(n, 2ngens(M), 2)
@@ -44,15 +67,15 @@ function PRASampler(
             S = union!(S, inv.(S))
         end
         append!(S, rand(rng, S, l - length(S)))
-        PRASampler(S, one(M), one(M))
+        ProductReplacementSampler(S, one(M), one(M))
     end
     for _ in 1:scramble_time
         _ = rand(rng, sampler)
     end
     return sampler
 end
 
-function Random.rand(rng::Random.AbstractRNG, pra::PRASampler)
+function Random.rand(rng::Random.AbstractRNG, pra::ProductReplacementSampler)
     i = rand(rng, 1:length(pra.gentuple))
 
     pra.right = pra.right * rand(rng, pra.gentuple)
@@ -61,3 +84,42 @@ function Random.rand(rng::Random.AbstractRNG, pra::PRASampler)
 
     return pra.left
 end
+
+struct Poisson
+    λ::Int
+    cdf::Vector{Float64}
+
+    function Poisson(; λ::Integer)
+        cdf = cumsum([λ^k * exp(-λ) / factorial(k) for k in 0:20])
+        return new(λ, cdf)
+    end
+end
+
+cdf(d::Poisson, x) = something(findfirst(≥(x), d.cdf), length(d.cdf) + 1)
+
+"""
+    RandomWordSampler(S, distribution)
+Return elements from a monoid represented by words of length in `S`, obeying `distribution`.
+
+Usually for a monoid (group) `S` is a (symmetric) generating set.
+
+`distribution` object must implement `cdf(distribution, x::Float64)::Integer`.
+
+For finite monoids or groups when uniformity of results is needed
+[`ProductReplacementSampler`](@ref) should be used.
+"""
+struct RandomWordSampler{V,D}
+    gens::V
+    distribution::D
+end
+
+Base.eltype(::Type{<:RandomWordSampler{V}}) where {V} = eltype(V)
+
+function Base.rand(
+    rng::Random.AbstractRNG,
+    rs::Random.SamplerTrivial{<:RandomWordSampler{I}},
+) where {I}
+    rw_sampler = rs[]
+    k = cdf(rw_sampler.distribution, rand(rng))
+    return prod(rand(rng, rw_sampler.gens, k))
+end