Add support for rng

Author: juliohm

Project.toml

@@ -7,6 +7,7 @@ version = "0.3.15"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 GeoStatsBase = "323cb8eb-fbf6-51c0-afd0-f8fba70507b2"
 Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [compat]
 Distributions = "0.21, 0.22, 0.23, 0.24, 0.25"

src/PointPatterns.jl

@@ -9,6 +9,8 @@ using GeoStatsBase
 
 using Distributions
 
+import Random
+
 include("processes.jl")
 include("thinning.jl")
 

src/processes.jl

@@ -17,35 +17,46 @@ Tells whether or not the spatial point process `p` is homogeneous.
 ishomogeneous(p::PointProcess) = false
 
 """
-    rand(p, g, n=1; [algo])
+    rand([rng], p, g, n=1; [algo])
 
 Generate `n` realizations of spatial point process `p`
 inside geometry `g`. Optionally specify sampling
-algorithm `algo`.
+algorithm `algo` and random number generator `rng`.
 """
-Base.rand(p::PointProcess, g::Geometry, n::Int;
+Base.rand(rng::Random.AbstractRNG,
+          p::PointProcess, g::Geometry, n::Int;
           algo=default_sampling_algorithm(p)) =
-  [rand_single(p, g, algo) for i in 1:n]
+  [rand_single(rng, p, g, algo) for i in 1:n]
 
-Base.rand(p::PointProcess, g::Geometry;
+Base.rand(rng::Random.AbstractRNG,
+          p::PointProcess, g::Geometry;
           algo=default_sampling_algorithm(p)) =
-  rand_single(p, g, algo)
+  rand_single(rng, p, g, algo)
+
+Base.rand(p::PointProcess, g::Geometry, n::Int;
+          algo=default_sampling_algorithm(p)) = 
+  rand(Random.GLOBAL_RNG, p, g, n; algo=algo)
+
+Base.rand(p::PointProcess, g::Geometry;
+          algo=default_sampling_algorithm(p)) = 
+  rand(Random.GLOBAL_RNG, p, g; algo=algo)
 
 """
     rand_single(p, g, algo)
 
 Generate a single realization of spatial point process
 `p` inside geometry `g` with sampling `algo`.
 """
-rand_single(::PointProcess, ::Geometry, algo) =
-  @error "not implemented"
+rand_single(::Random.AbstractRNG, ::PointProcess, ::Geometry, algo) =
+  throw(ErrorException("not implemented"))
 
 """
     default_sampling_algorithm(p)
 
 Default sampling algorithm for spatial point process `p`.
 """
-default_sampling_algorithm(::PointProcess) = @error "not implemented"
+default_sampling_algorithm(::PointProcess) =
+  throw(ErrorException("not implemented"))
 
 """
     p₁ ∪ p₂

src/processes/binomial.jl

@@ -17,7 +17,8 @@ struct BinomialSampling end
 
 default_sampling_algorithm(::BinomialProcess) = BinomialSampling()
 
-function rand_single(p::BinomialProcess, b::Box{Dim,T},
+function rand_single(rng::Random.AbstractRNG,
+                     p::BinomialProcess, b::Box{Dim,T},
                      ::BinomialSampling) where {Dim,T}
   # region configuration
   lo, up = coordinates.(extrema(b))
@@ -26,5 +27,5 @@ function rand_single(p::BinomialProcess, b::Box{Dim,T},
   U = product_distribution([Uniform(lo[i], up[i]) for i in 1:Dim])
 
   # return point pattern
-  PointSet(rand(U, p.n))
+  PointSet(rand(rng, U, p.n))
 end

src/processes/poisson.jl

@@ -24,26 +24,30 @@ default_sampling_algorithm(::PoissonProcess{<:Function}) = DiscretizedSampling()
 #------------------
 # HOMOGENEOUS CASE
 #------------------
-function rand_single(p::PoissonProcess{<:Real}, b::Box{Dim,T},
+
+function rand_single(rng::Random.AbstractRNG,
+                     p::PoissonProcess{<:Real}, b::Box{Dim,T},
                      ::ProductSampling) where {Dim,T}
   # region configuration
   lo, up = coordinates.(extrema(b))
 
   # simulate number of points
   λ = p.λ; V = measure(b)
-  n = rand(Poisson(λ*V))
+  n = rand(rng, Poisson(λ*V))
 
   # product of uniform distributions
   U = product_distribution([Uniform(lo[i], up[i]) for i in 1:Dim])
 
   # return point pattern
-  PointSet(rand(U, n))
+  PointSet(rand(rng, U, n))
 end
 
 #--------------------
 # INHOMOGENEOUS CASE
 #--------------------
-function rand_single(p::PoissonProcess{<:Function}, g::Box{Dim,T},
+
+function rand_single(rng::Random.AbstractRNG,
+                     p::PoissonProcess{<:Function}, g::Box{Dim,T},
                      algo::DiscretizedSampling) where {Dim,T}
   # discretize region
   # TODO
@@ -55,7 +59,8 @@ function rand_single(p::PoissonProcess{<:Function}, g::Box{Dim,T},
   # TODO
 end
 
-function rand_single(p::PoissonProcess{<:Function}, b::Box{Dim,T},
+function rand_single(rng::Random.AbstractRNG,
+                     p::PoissonProcess{<:Function}, b::Box{Dim,T},
                      algo::ThinnedSampling) where {Dim,T}
   # Lewis-Shedler algorithm
   # TODO

src/processes/union.jl

@@ -18,12 +18,11 @@ struct UnionSampling end
 
 default_sampling_algorithm(::UnionProcess) = UnionSampling()
 
-function rand_single(p::UnionProcess, g::Geometry, algo::UnionSampling)
-  pp₁ = rand(p.p₁, g)
-  pp₂ = rand(p.p₂, g)
+function rand_single(rng::Random.AbstractRNG,
+                     p::UnionProcess, g::Geometry,
+                     ::UnionSampling)
+  pp₁ = rand(rng, p.p₁, g)
+  pp₂ = rand(rng, p.p₂, g)
 
-  X = [coordinates(pp₁[i]) for i in 1:nelements(pp₁)]
-  Y = [coordinates(pp₂[i]) for i in 1:nelements(pp₂)]
-
-  PointSet([X; Y])
+  PointSet([coordinates.(pp₁); coordinates.(pp₂)])
 end

src/thinning/random.jl

@@ -12,8 +12,9 @@ struct RandomThinning{P<:Union{Real,Function}} <: AbstractThinning
 end
 
 # -----------------------
-# thinning point process
+# THINNING POINT PROCESS
 # -----------------------
+
 thin(p::PoissonProcess{<:Real}, t::RandomThinning{<:Real}) =
   PoissonProcess(t.p * p.λ)
 
@@ -27,8 +28,9 @@ thin(p::PoissonProcess{<:Function}, t::RandomThinning{<:Real}) =
   PoissonProcess(x -> t.p * p.λ(x))
 
 # -----------------------
-# thinning point pattern
+# THINNING POINT PATTERN
 # -----------------------
+
 function thin(pp::PointSet, t::RandomThinning{<:Real})
   draws = rand(Bernoulli(t.p), nelements(pp))
   inds  = findall(isequal(1), draws)

test/processes.jl

@@ -1,9 +1,9 @@
 @testset "Processes" begin
   @testset "Basic" begin
     for p in [BinomialProcess(100), PoissonProcess(100.)]
-      r = Box((0.,1.), (1.,2.))
-      pp = rand(p, r)
-      xs = [coordinates(pp[i]) for i in 1:nelements(pp)]
+      b  = Box((0.,1.), (1.,2.))
+      pp = rand(p, b)
+      xs = coordinates.(pp)
       @test all(0 .≤ first.(xs) .≤ 1)
       @test all(1 .≤ last.(xs) .≤ 2)
     end