Implement rand()

Project.toml

@@ -6,17 +6,20 @@ version = "0.5.12"
 CompositeTypes = "b152e2b5-7a66-4b01-a709-34e65c35f657"
 IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
 CompositeTypes = "0.1.2"
 IntervalSets = "0.5, 0.6, 0.7"
 StaticArrays = "0.12.2, 1"
+StableRNGs = "1"
 julia = "1.6"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 
 [targets]
-test = ["Test"]
+test = ["Test", "StableRNGs"]

src/DomainSets.jl

@@ -3,6 +3,8 @@ module DomainSets
 using StaticArrays
 using LinearAlgebra, Statistics
 import LinearAlgebra: cross, ×, pinv
+import Random
+using Random: AbstractRNG
 
 using IntervalSets
 using CompositeTypes, CompositeTypes.Display, CompositeTypes.Indexing
@@ -216,6 +218,7 @@ include("domains/indicator.jl")
 include("domains/boundingbox.jl")
 
 include("applications/coordinates.jl")
+include("applications/random.jl")
 include("applications/rotation.jl")
 
 end # module

src/applications/random.jl

@@ -0,0 +1,38 @@
+Random.gentype(::Type{<:Domain{T}}) where T = T
+
+Base.rand(rng::AbstractRNG, s::Random.SamplerTrivial{<:ProductDomain}) = toexternalpoint(s[], map(i->(rand(rng, i)), factors(s[])))
+
+Base.rand(rng::AbstractRNG, s::Random.SamplerTrivial{<:SimpleLazyDomain}) = toexternalpoint(s[], rand(rng, superdomain(s[])))
+
+function Base.rand(rng::AbstractRNG, s::Random.SamplerTrivial{<:Ball})
+    # Technical details: http://extremelearning.com.au/how-to-generate-uniformly-random-points-on-n-spheres-and-n-balls/
+
+    b = s[]
+
+    # for low-dimensional balls, use rejection sampling - acceptance rate is at least 52%
+    if dimension(b) <= 3
+        while true
+            r = rand(rng, boundingbox(b))
+            if r in b
+                return r
+            end
+        end
+
+    # for higher dimensional balls, use the "Mueller" method
+    else
+        u = randn_dimension(rng, eltype(b), dimension(b))
+        r = radius(b)*rand(rng)^(1/dimension(b))
+        return (r/norm(u))*u + center(b)
+    end
+end
+
+randn_dimension(rng::AbstractRNG, t::Type{<:StaticVector}, d) = randn(rng, t)
+randn_dimension(rng::AbstractRNG, t::Type{<:Vector}, d) = randn(rng, eltype(t), d)
+
+# Implementation notes
+# ====================
+#
+# The methods implemented above are the easy ones
+# Unions and intersections could be implemented with rejection sampling, but it might be inefficient
+# Sphere will require some decisions because `rand(sphere) in sphere` will usually only be approximately satisfied
+# Maps may be difficult because the map could distort the distribution so that it is not uniform.

test/runtests.jl

@@ -1,4 +1,4 @@
-using Test, LinearAlgebra, StaticArrays
+using Test, LinearAlgebra, StaticArrays, Random, StableRNGs
 
 using DomainSets
 using CompositeTypes.Indexing

test/test_applications.jl

@@ -48,11 +48,84 @@ function test_cart_polar_map(T)
     @test inverse(m2) == m1
 end
 
+function test_rand(T)
+    r = Rectangle(T[-1, 2], T[3, 4])
+    @test @inferred(Random.gentype(r)) == Vector{T}
+    @test typeof(rand(r)) == Random.gentype(r)
+    @test @inferred(rand(r)) in r
 
+    r = Rectangle(SA[T(-1), T(2)], SA[T(3), T(4)])
+    @test @inferred(Random.gentype(r)) == SVector{2, T}
+    @test typeof(rand(r)) == Random.gentype(r)
+    @test @inferred(rand(r)) in r
+
+    hybrid_product = ProductDomain(["a", "b"], T(1.0)..T(2.0))
+    @test @inferred(Random.gentype(hybrid_product)) == Tuple{String, T}
+    @test typeof(rand(hybrid_product)) == Random.gentype(hybrid_product)
+    @test @inferred(rand(hybrid_product)) in hybrid_product
+
+    b = Ball(2.0, SA[T(1.0), T(2.0)])
+    @test @inferred(Random.gentype(b)) == SVector{2, T}
+    if T == BigFloat
+        @test_broken typeof(rand(b)) == Random.gentype(b) # all rand tests for BigFloat StaticVector Balls are broken due to issue #114
+    else
+        @test typeof(rand(b)) == Random.gentype(b)
+        @test @inferred(rand(b)) in b
+        @test all(p in b for p in rand(b, 100))
+        test_rng_consistency(b)
+    end
+
+    b = Ball(2.0, [T(1.0), T(2.0)])
+    @test @inferred(Random.gentype(b)) == Vector{T}
+    @test typeof(rand(b)) == Random.gentype(b)
+    @test @inferred(rand(b)) in b
+    @test all(p in b for p in rand(b, 100))
+    test_rng_consistency(b)
+
+    b = Ball(2.0, T(1.0))
+    @test @inferred(Random.gentype(b)) == T
+    @test_broken typeof(rand(b)) == Random.gentype(b) # broken due to issue #113
+    @test_broken @inferred(rand(b)) in b # broken due to issue #113
+    # test_rng_consistency(b) # uncomment when issue #113 is fixed
+
+    # Higher dimension
+    b = Ball(2.0, SA[T(1.0), T(1.0), T(1.0), T(1.0)])
+    @test @inferred(Random.gentype(b)) == SVector{4, T}
+    if T == BigFloat
+        @test_broken typeof(rand(b)) == Random.gentype(b) # all rand tests for BigFloat StaticVector Balls are broken due to issue #114
+    else
+        @test typeof(rand(b)) == Random.gentype(b)
+        @test @inferred(rand(b)) in b
+        @test all(p in b for p in rand(b, 100))
+        test_rng_consistency(b)
+    end
+
+    if T == Float64
+        # Test numerical accuracy - two rectangles of the same size should have the same number of points
+        # Only works for Float64 since StableRNG doesn't generate BigFloats
+        b = Ball(1.5, SA[1.0, -1.0, 2.0, -3.0])
+        rng = StableRNG(1)
+        n = 1_000_000
+        region_1 = Rectangle([0.0, -0.3, 0.0, -0.3], [0.3, 0.0, 0.3, 0.0]) .+ center(b)
+        lower_corner = radius(b)*fill(-0.5, 4)
+        region_2 = Rectangle(lower_corner, lower_corner.+0.3) .+ center(b)
+        rs = rand(rng, b, n)
+        n_1 = sum(r in region_1 for r in rs)
+        n_2 = sum(r in region_2 for r in rs)
+        @test isapprox(n_1, n_2, rtol=0.1) 
+    end
+end
+
+function test_rng_consistency(set)
+    rng1 = Random.MersenneTwister(1)
+    rng2 = Random.MersenneTwister(1)
+    @test rand(rng1, set) == rand(rng2, set)
+end
 
 function test_applications(T)
     test_rotation_map(T)
     test_cart_polar_map(T)
+    test_rand(T)
 end
 
 @testset "applications" begin