Refactor rand to include multiple algorithms

Author: juliohm

src/processes.jl

@@ -17,23 +17,35 @@ Tells whether or not the spatial point process `p` is homogeneous.
 ishomogeneous(p::PointProcess) = false
 
 """
-    rand(p, r, n=1)
+    rand(p, r, n=1; [algo])
 
 Generate `n` realizations of spatial point process `p`
-inside spatial region `r`.
+inside spatial region `r`. Optionally specify sampling
+algorithm `algo`.
 """
-Base.rand(p::PointProcess, r::AbstractRegion, n::Int) =
-  [rand_single(p, r) for i in 1:n]
+Base.rand(p::PointProcess, r::AbstractRegion, n::Int;
+          algo=default_sampling_algorithm(p)) =
+  [rand_single(p, r, algo) for i in 1:n]
 
-Base.rand(p::PointProcess, r::AbstractRegion) = rand_single(p, r)
+Base.rand(p::PointProcess, r::AbstractRegion;
+          algo=default_sampling_algorithm(p)) =
+  rand_single(p, r, algo)
 
 """
-    rand_single(p, r)
+    rand_single(p, r, algo)
 
 Generate a single realization of spatial point process
-`p` inside spatial region `r`.
+`p` inside spatial region `r` with sampling `algo`.
 """
-rand_single(p::PointProcess, r::AbstractRegion) = @error "not implemented"
+rand_single(p::PointProcess, r::AbstractRegion, algo) =
+  @error "not implemented"
+
+"""
+    default_sampling_algorithm(p)
+
+Default sampling algorithm for spatial point process `p`.
+"""
+default_sampling_algorithm(p::PointProcess) = @error "not implemented"
 
 """
     p₁ ∪ p₂

src/processes/binomial.jl

@@ -13,7 +13,12 @@ end
 
 ishomogeneous(p::BinomialProcess) = true
 
-function rand_single(p::BinomialProcess, r::RectangleRegion{T,N}) where {N,T}
+struct BinomialSampling end
+
+default_sampling_algorithm(::BinomialProcess) = BinomialSampling()
+
+function rand_single(p::BinomialProcess, r::RectangleRegion{T,N},
+                     algo::BinomialSampling) where {N,T}
   # region configuration
   lo = lowerleft(r)
   up = upperright(r)

src/processes/poisson.jl

@@ -12,9 +12,20 @@ struct PoissonProcess{L<:Union{Real,Function}} <: PointProcess
 end
 
 ishomogeneous(p::PoissonProcess{<:Real}) = true
+ishomogeneous(p::PoissonProcess{<:Function}) = false
 
+struct ProductSampling end
+struct ThinnedSampling end
+struct DiscretizedSampling end
+
+default_sampling_algorithm(::PoissonProcess{<:Real}) = ProductSampling()
+default_sampling_algorithm(::PoissonProcess{<:Function}) = DiscretizedSampling()
+
+# -----------------
 # homogeneous case
-function rand_single(p::PoissonProcess{<:Real}, r::RectangleRegion{T,N}) where {N,T}
+# -----------------
+function rand_single(p::PoissonProcess{<:Real}, r::RectangleRegion{T,N},
+                     algo::ProductSampling) where {N,T}
   # region configuration
   lo = lowerleft(r)
   up = upperright(r)
@@ -30,9 +41,35 @@ function rand_single(p::PoissonProcess{<:Real}, r::RectangleRegion{T,N}) where {
   PointPattern(rand(U, n))
 end
 
+# -------------------
 # inhomogeneous case
-function rand_single(p::PoissonProcess{<:Function}, r::RectangleRegion{T,N}) where {N,T}
+# -------------------
+function rand_single(p::PoissonProcess{<:Function}, r::RectangleRegion{T,N},
+                     algo::DiscretizedSampling) where {N,T}
+  # discretize region
+  # TODO
+
+  # discretize retention
+  # TODO
+
+  # sample each element
   # TODO
 end
 
-Base.union(p₁::PoissonProcess, p₂::PoissonProcess) = PoissonProcess(p₁.λ + p₂.λ)
+function rand_single(p::PoissonProcess{<:Function}, r::RectangleRegion{T,N},
+                     algo::ThinnedSampling) where {N,T}
+  # Lewis-Shedler algorithm
+  # TODO
+end
+
+Base.union(p₁::PoissonProcess{<:Real}, p₂::PoissonProcess{<:Real}) =
+  PoissonProcess(p₁.λ + p₂.λ)
+
+Base.union(p₁::PoissonProcess{<:Function}, p₂::PoissonProcess{<:Function}) =
+  PoissonProcess(x -> p₁.λ(x) + p₂.λ(x))
+
+Base.union(p₁::PoissonProcess{<:Real}, p₂::PoissonProcess{<:Function}) =
+  PoissonProcess(x -> p₁.λ + p₂.λ(x))
+
+Base.union(p₁::PoissonProcess{<:Function}, p₂::PoissonProcess{<:Real}) =
+  PoissonProcess(x -> p₁.λ(x) + p₂.λ)

src/processes/union.jl

@@ -14,7 +14,11 @@ end
 
 ishomogeneous(p::UnionProcess) = ishomogeneous(p.p₁) && ishomogeneous(p.p₂)
 
-function rand_single(p::UnionProcess, r::AbstractRegion)
+struct UnionSampling end
+
+default_sampling_algorithm(::UnionProcess) = UnionSampling()
+
+function rand_single(p::UnionProcess, r::AbstractRegion, algo::UnionSampling)
   pp₁ = rand(p.p₁, r)
   pp₂ = rand(p.p₂, r)
 

src/thinning/random.jl

@@ -18,13 +18,13 @@ thin(p::PoissonProcess{<:Real}, t::RandomThinning{<:Real}) =
   PoissonProcess(t.p * p.λ)
 
 thin(p::PoissonProcess{<:Function}, t::RandomThinning{<:Function}) =
-  PoissonProcess(u -> t.p(u) * p.λ(u))
+  PoissonProcess(x -> t.p(x) * p.λ(x))
 
 thin(p::PoissonProcess{<:Real}, t::RandomThinning{<:Function}) =
-  PoissonProcess(u -> t.p(u) * p.λ)
+  PoissonProcess(x -> t.p(x) * p.λ)
 
 thin(p::PoissonProcess{<:Function}, t::RandomThinning{<:Real}) =
-  PoissonProcess(u -> t.p * p.λ(u))
+  PoissonProcess(x -> t.p * p.λ(x))
 
 # -----------------------
 # thinning point pattern