Merge pull request #38 from JuliaEarth/refactor

Erase algorithms

Author: juliohm

src/PointPatterns.jl

@@ -23,11 +23,6 @@ export
   UnionProcess,
   ishomogeneous,
 
-  # algorithms
-  PointPatternAlgo,
-  ConstantIntensity,
-  LewisShedler,
-
   # thinning methods
   ThinningMethod,
   RandomThinning,

src/processes.jl

@@ -17,68 +17,28 @@ Tells whether or not the spatial point process `p` is homogeneous.
 ishomogeneous(p::PointProcess) = false
 
 """
-    rand([rng], p, g, n=1; [algo])
+    rand([rng], p, g, n=1)
 
 Generate `n` realizations of spatial point process `p`
-inside geometry or domain `g`. Optionally specify sampling
-algorithm `algo` and random number generator `rng`.
+inside geometry or domain `g`. Optionally specify the
+random number generator `rng`.
 """
-Base.rand(rng::Random.AbstractRNG, p::PointProcess, g, n::Int; algo=defaultalgo(p, g)) =
-  [randsingle(rng, p, g, algo) for i in 1:n]
+Base.rand(rng::Random.AbstractRNG, p::PointProcess, g, n::Int) = [randsingle(rng, p, g) for _ in 1:n]
 
-Base.rand(rng::Random.AbstractRNG, p::PointProcess, g; algo=defaultalgo(p, g)) =
-  randsingle(rng, p, g, algo)
+Base.rand(rng::Random.AbstractRNG, p::PointProcess, g) = randsingle(rng, p, g)
 
-Base.rand(p::PointProcess, g, n::Int; algo=defaultalgo(p, g)) =
-  rand(Random.GLOBAL_RNG, p, g, n; algo=algo)
+Base.rand(p::PointProcess, g, n::Int) = rand(Random.GLOBAL_RNG, p, g, n)
 
-Base.rand(p::PointProcess, g; algo=defaultalgo(p, g)) = rand(Random.GLOBAL_RNG, p, g; algo=algo)
+Base.rand(p::PointProcess, g) = rand(Random.GLOBAL_RNG, p, g)
 
 """
-    randsingle(rng, p, g, algo)
+    randsingle(rng, p, g)
 
-Generate a single realization of spatial point process
-`p` inside geometry or domain `g` with sampling `algo`.
+Generate a single point pattern of point process `p`
+inside geometry or domain `g`.
 """
 function randsingle end
 
-"""
-    defaultalgo(p, g)
-
-Default sampling algorithm for spatial point process `p`
-on geometry or domain `g`.
-"""
-function defaultalgo end
-
-# -------------------------
-# POINT PATTERN ALGORITHMS
-# -------------------------
-
-"""
-    PointPatternAlgo
-
-A method for sampling point patterns.
-"""
-abstract type PointPatternAlgo end
-
-"""
-    LewisShedler(λmax)
-
-Generate sample using Lewis-Shedler algorithm (1979) with
-maximum real value `λmax` of the intensity function.
-"""
-struct LewisShedler{T<:Real} <: PointPatternAlgo
-  λmax::T
-end
-
-"""
-    ConstantIntensity()
-
-Generate sample assuming the intensity is constant over a `Geometry`
-or piecewise constant over a `Domain`.
-"""
-struct ConstantIntensity <: PointPatternAlgo end
-
 #-----------------
 # IMPLEMENTATIONS
 #-----------------

src/processes/binomial.jl

@@ -13,7 +13,5 @@ end
 
 ishomogeneous(p::BinomialProcess) = true
 
-defaultalgo(::BinomialProcess, ::Any) = ConstantIntensity()
-
-randsingle(rng::Random.AbstractRNG, p::BinomialProcess, g, ::ConstantIntensity) =
+randsingle(rng::Random.AbstractRNG, p::BinomialProcess, g) =
   PointSet(sample(rng, g, HomogeneousSampling(p.n)))

src/processes/cluster.jl

@@ -29,9 +29,7 @@ end
 ClusterProcess(p::PointProcess, o::PointProcess, gfun::Function) =
   ClusterProcess(p, parent -> rand(o, gfun(parent)))
 
-defaultalgo(::ClusterProcess, ::Any) = nothing
-
-function randsingle(rng::Random.AbstractRNG, p::ClusterProcess, g, ::Nothing)
+function randsingle(rng::Random.AbstractRNG, p::ClusterProcess, g)
   # retrieve parameters
   Dim = embeddim(g)
   T = coordtype(g)
@@ -50,4 +48,3 @@ function randsingle(rng::Random.AbstractRNG, p::ClusterProcess, g, ::Nothing)
   # return point pattern
   PointSet(points)
 end
-

src/processes/poisson.jl

@@ -21,21 +21,11 @@ ishomogeneous(p::PoissonProcess{<:Function}) = false
 
 ishomogeneous(p::PoissonProcess{<:AbstractVector}) = false
 
-defaultalgo(::PoissonProcess, ::Any) = ConstantIntensity()
-
-defaultalgo(p::PoissonProcess{<:Function}, g) = LewisShedler(maxintensity(p, g))
-
-function maxintensity(p::PoissonProcess{<:Function}, g)
-  points = sample(g, HomogeneousSampling(10000))
-  λmin, λmax = extrema(p.λ, points)
-  λmax + 0.05 * (λmax - λmin)
-end
-
 #------------------
 # HOMOGENEOUS CASE
 #------------------
 
-function randsingle(rng::Random.AbstractRNG, p::PoissonProcess{<:Real}, g, ::ConstantIntensity)
+function randsingle(rng::Random.AbstractRNG, p::PoissonProcess{<:Real}, g)
   # simulate number of points
   λ = p.λ
   V = measure(g)
@@ -49,23 +39,31 @@ end
 # INHOMOGENEOUS CASE
 #--------------------
 
-function randsingle(rng::Random.AbstractRNG, p::PoissonProcess{<:Function}, g, algo::LewisShedler)
+function randsingle(rng::Random.AbstractRNG, p::PoissonProcess{<:Function}, g)
+  # upper bound for intensity
+  λmax = maxintensity(p, g)
+
   # simulate a homogeneous process
-  pp = randsingle(rng, PoissonProcess(algo.λmax), g, ConstantIntensity())
+  pset = randsingle(rng, PoissonProcess(λmax), g)
 
   # thin point pattern
-  isnothing(pp) ? nothing : PointSet(collect(thin(pp, RandomThinning(x -> p.λ(x) / algo.λmax))))
+  isnothing(pset) ? nothing : thin(pset, RandomThinning(x -> p.λ(x) / λmax))
 end
 
-randsingle(rng::Random.AbstractRNG, p::PoissonProcess{<:Function}, d::Domain, algo::ConstantIntensity) =
-  randsingle(rng, PoissonProcess(p.λ.(centroid.(d))), d, algo)
+randsingle(rng::Random.AbstractRNG, p::PoissonProcess{<:Function}, d::Domain) =
+  randsingle(rng, PoissonProcess(p.λ.(centroid.(d))), d)
 
-function randsingle(rng::Random.AbstractRNG, p::PoissonProcess{<:AbstractVector}, d::Domain, algo::ConstantIntensity)
+function randsingle(rng::Random.AbstractRNG, p::PoissonProcess{<:AbstractVector}, d::Domain)
   # simulate number of points
-  λ = p.λ
-  V = measure.(d)
-  n = rand(rng, Poisson(sum(λ .* V)))
+  λ = p.λ .* measure(d)
+  n = rand(rng, Poisson(sum(λ)))
 
-  # simulate n points
-  iszero(n) ? nothing : PointSet(sample(rng, d, HomogeneousSampling(n, λ .* V)))
+  # simulate point pattern
+  iszero(n) ? nothing : PointSet(sample(rng, d, HomogeneousSampling(n, λ)))
+end
+
+function maxintensity(p::PoissonProcess{<:Function}, g)
+  points = sample(g, HomogeneousSampling(10000))
+  λmin, λmax = extrema(p.λ, points)
+  λmax + 0.05 * (λmax - λmin)
 end

src/processes/union.jl

@@ -14,9 +14,7 @@ end
 
 ishomogeneous(p::UnionProcess) = ishomogeneous(p.p₁) && ishomogeneous(p.p₂)
 
-defaultalgo(::UnionProcess, ::Any) = nothing
-
-function randsingle(rng::Random.AbstractRNG, p::UnionProcess, g, ::Nothing)
+function randsingle(rng::Random.AbstractRNG, p::UnionProcess, g)
   pp₁ = rand(rng, p.p₁, g)
   pp₂ = rand(rng, p.p₂, g)
   PointSet([collect(pp₁); collect(pp₂)])

test/processes.jl

@@ -1,77 +1,57 @@
 @testset "Processes" begin
+  # geometries and domains
   seg = Segment((0.0, 0.0), (11.3, 11.3))
   tri = Triangle((0.0, 0.0), (5.65, 0.0), (5.65, 5.65))
   quad = Quadrangle((0.0, 0.0), (0.0, 4.0), (4.0, 4.0), (4.0, 0.0))
   box = Box((0.0, 0.0), (4.0, 4.0))
   ball = Ball((1.0, 1.0), 2.25)
-  outer = Point2[(0, -4), (4, -1), (4, 1.5), (0, 3)]
-  hole1 = Point2[(0.2, -0.2), (1.4, -0.2), (1.4, 0.6), (0.2, 0.6)]
-  hole2 = Point2[(2, -0.2), (3, -0.2), (3, 0.4), (2, 0.4)]
+  outer = [(0, -4), (4, -1), (4, 1.5), (0, 3)]
+  hole1 = [(0.2, -0.2), (1.4, -0.2), (1.4, 0.6), (0.2, 0.6)]
+  hole2 = [(2, -0.2), (3, -0.2), (3, 0.4), (2, 0.4)]
   poly = PolyArea(outer, [hole1, hole2])
   grid = CartesianGrid((0, 0), (4, 4), dims=(10, 10))
   points = Point2[(0, 0), (4.5, 0), (0, 4.2), (4, 4.3), (1.5, 1.5)]
   connec = connect.([(1, 2, 5), (2, 4, 5), (4, 3, 5), (3, 1, 5)], Triangle)
   mesh = SimpleMesh(points, connec)
+  geoms = [seg, tri, quad, box, ball, poly, grid, mesh]
+
+  # point processes
+  λ(s) = sum(coordinates(s) .^ 2)
+  binom = BinomialProcess(100)
+  poisson1 = PoissonProcess(100.0)
+  poisson2 = PoissonProcess(λ)
+  procs = [binom, poisson1, poisson2]
 
   @testset "Basic" begin
-    for p in [BinomialProcess(100), PoissonProcess(100.0)]
-      b = Box((0.0, 1.0), (1.0, 2.0))
-      pp = rand(p, b)
-      xs = coordinates.(pp)
-      @test all(0 .≤ first.(xs) .≤ 1)
-      @test all(1 .≤ last.(xs) .≤ 2)
+    for p in procs, g in geoms
+      pp = rand(p, g)
+      @test all(∈(g), pp)
     end
   end
 
   @testset "Binomial" begin
     p = BinomialProcess(10)
-    for g in [seg, tri, quad, box, ball, poly, grid, mesh]
+    for g in geoms
       pp = rand(p, g)
       @test nelements(pp) == 10
-      @test all(∈(g), pp)
     end
   end
 
   @testset "Poisson" begin
-    # helper function
-    λ(s::Point2) = sum(coordinates(s) .^ 2)
-
-    # homogeneous
-    p = PoissonProcess(10.0)
-    for g in [seg, tri, quad, box, ball, poly, grid, mesh]
-      pp = rand(p, g)
-      @test all(∈(g), pp)
-    end
-
-    # inhomogeneous with intensity function
-    p = PoissonProcess(λ)
-    for g in [seg, tri, quad, box, ball, poly, grid, mesh]
-      pp = rand(p, g)
-      @test all(∈(g), pp)
-
-      pp = rand(p, g, algo=LewisShedler(λ(Point2(12.0, 12.0))))
-      @test all(∈(g), pp)
-
-      g = discretize(g)
-      pp = rand(p, g, algo=ConstantIntensity())
-      @test all(∈(g), g)
-    end
-
     # inhomogeneous with piecewise constant intensity
     for g in [grid, mesh]
       λvec = λ.(centroid.(g))
       p = PoissonProcess(λvec)
-      # discretizedsampling by default
       pp = rand(p, g)
       @test all(∈(g), pp)
     end
 
     # empty pointsets
-    for g in [seg, tri, quad, box, ball, poly, grid, mesh]
+    for g in geoms
       @test isnothing(rand(PoissonProcess(0.0), seg))
     end
-    ps = PointSet(rand(Point2, 10))
-    @test isnothing(rand(PoissonProcess(100.0), ps))
+    pp = PointSet(rand(Point2, 10))
+    @test isnothing(rand(PoissonProcess(100.0), pp))
   end
 
   @testset "Union" begin