Inhomogeneous (#32)

* sampling for inhomogeneous poisson process

* use rand instead of rand_single

* remove all method for box

* export thinnedsampling

* Sampling over a Geometry or Mesh

* remove discretized sampling for inhomogeneous case

* remove completely productsampling

* point process with piecewise constant intensity

* discretizedsampling for inhomogeneous process

* documentation for poisson process

* more methods when lambda is a vector

* fix grammar

* Apply suggestions from code review

Space where to sample without particular type.

Co-authored-by: Júlio Hoffimann <[email protected]>

* fix tests

* return nothing when no points need to be sampled

* discrete sampling with intensity function over mesh

* default sampling depending only on pointprocess

* alternative SimpleThinnedSampling

* fix grammar

* adding tests for sampling PoissonProcess

* remove comment

* lowercase comments

* Update src/processes/binomial.jl

Co-authored-by: Júlio Hoffimann <[email protected]>

* Update src/processes/poisson.jl

Co-authored-by: Júlio Hoffimann <[email protected]>

* Update src/processes/poisson.jl

Co-authored-by: Júlio Hoffimann <[email protected]>

* Update src/processes/poisson.jl

Co-authored-by: Júlio Hoffimann <[email protected]>

* Apply suggestions from code review

Co-authored-by: Júlio Hoffimann <[email protected]>

* remove collect

* remove simplethinnedsampling

* remove export simplethinnedsampling

* empty

* simplify thinnedsampling

* abstractvector

* prepare to not export sampling algorithms

* Do not export sampling methods

Co-authored-by: Júlio Hoffimann <[email protected]>

* default_sampling depends also of geometry

* use rand_single inside rand_single

* documenting DiscretizedSampling

* defaults sampling depends of geometry

* keep same order for sampling methods

---------

Co-authored-by: Júlio Hoffimann <[email protected]>

Author: ErickChacon

src/processes.jl

@@ -23,16 +23,16 @@ Generate `n` realizations of spatial point process `p`
 inside geometry `g`. Optionally specify sampling
 algorithm `algo` and random number generator `rng`.
 """
-Base.rand(rng::Random.AbstractRNG, p::PointProcess, g::Geometry, n::Int; algo=default_sampling_algorithm(p, g)) =
+Base.rand(rng::Random.AbstractRNG, p::PointProcess, g, n::Int; algo=default_sampling_algorithm(p, g)) =
   [rand_single(rng, p, g, algo) for i in 1:n]
 
-Base.rand(rng::Random.AbstractRNG, p::PointProcess, g::Geometry; algo=default_sampling_algorithm(p, g)) =
+Base.rand(rng::Random.AbstractRNG, p::PointProcess, g; algo=default_sampling_algorithm(p, g)) =
   rand_single(rng, p, g, algo)
 
-Base.rand(p::PointProcess, g::Geometry, n::Int; algo=default_sampling_algorithm(p, g)) =
+Base.rand(p::PointProcess, g, n::Int; algo=default_sampling_algorithm(p, g)) =
   rand(Random.GLOBAL_RNG, p, g, n; algo=algo)
 
-Base.rand(p::PointProcess, g::Geometry; algo=default_sampling_algorithm(p, g)) =
+Base.rand(p::PointProcess, g; algo=default_sampling_algorithm(p, g)) =
   rand(Random.GLOBAL_RNG, p, g; algo=algo)
 
 """
@@ -47,17 +47,32 @@ function rand_single end
     default_sampling_algorithm(p, g)
 
 Default sampling algorithm for spatial point process `p`
-on geometry `g`.
+on geometry or domain `g`.
 """
 function default_sampling_algorithm end
 
 # --------------------
 # SAMPLING ALGORITHMS
 # --------------------
 
-struct ProductSampling end
-struct ThinnedSampling end
+"""
+    ThinnedSampling(λmax)
+
+Generate sample using Lewis-Shedler algorithm (1979) with
+maximum real value `λmax` of the intensity function.
+"""
+struct ThinnedSampling{T<:Real}
+  λmax::T
+end
+
+"""
+    DiscretizedSampling()
+
+Generate sample assuming the intensity is constant over a `Geometry`
+or piecewise constant over a `Domain`.
+"""
 struct DiscretizedSampling end
+
 struct UnionSampling end
 
 #-----------------

src/processes/binomial.jl

@@ -13,15 +13,9 @@ end
 
 ishomogeneous(p::BinomialProcess) = true
 
-default_sampling_algorithm(::BinomialProcess, ::Geometry) = ProductSampling()
+default_sampling_algorithm(::BinomialProcess, ::Any) = DiscretizedSampling()
 
-function rand_single(rng::Random.AbstractRNG, p::BinomialProcess, b::Box{Dim,T}, ::ProductSampling) where {Dim,T}
-  # region configuration
-  lo, up = coordinates.(extrema(b))
-
-  # product of uniform distributions
-  U = product_distribution([Uniform(lo[i], up[i]) for i in 1:Dim])
-
-  # return point pattern
-  PointSet(rand(rng, U, p.n))
+function rand_single(rng::Random.AbstractRNG, p::BinomialProcess, g, ::DiscretizedSampling)
+  points = sample(g, HomogeneousSampling(p.n))
+  PointSet(points)
 end

src/processes/poisson.jl

@@ -5,9 +5,14 @@
 """
    PoissonProcess(λ)
 
-A Poisson process with intensity `λ`.
+A Poisson process with intensity `λ`. For a homogeneous process,
+define `λ` as a constant real value, while for an inhomogeneous process,
+define `λ` as a function or vector of values. If `λ` is a vector, it is
+assumed that the process is associated with a `Domain` with the same
+number of elements as `λ`.
 """
-struct PoissonProcess{L<:Union{Real,Function}} <: PointProcess
+
+struct PoissonProcess{L<:Union{Real,Function,AbstractVector}} <: PointProcess
   λ::L
 end
 
@@ -19,61 +24,82 @@ Base.union(p₁::PoissonProcess{<:Real}, p₂::PoissonProcess{<:Function}) = Poi
 
 Base.union(p₁::PoissonProcess{<:Function}, p₂::PoissonProcess{<:Real}) = PoissonProcess(x -> p₁.λ(x) + p₂.λ)
 
+Base.union(p₁::PoissonProcess{<:AbstractVector}, p₂::PoissonProcess{<:AbstractVector}) = PoissonProcess(x -> p₁.λ + p₂.λ)
+
 ishomogeneous(p::PoissonProcess{<:Real}) = true
 ishomogeneous(p::PoissonProcess{<:Function}) = false
+ishomogeneous(p::PoissonProcess{<:AbstractVector}) = false
 
-default_sampling_algorithm(::PoissonProcess, ::Geometry) = DiscretizedSampling()
-default_sampling_algorithm(::PoissonProcess{<:Real}, ::Box) = ProductSampling()
+default_sampling_algorithm(::PoissonProcess, ::Any) = DiscretizedSampling()
+default_sampling_algorithm(p::PoissonProcess{<:Function}, g) = ThinnedSampling(default_lambda_max(p, g))
+
+function default_lambda_max(p::PoissonProcess{<:Function}, g)
+  points = sample(g, HomogeneousSampling(10000))
+  λvec = p.λ.(points)
+  maximum(λvec) + 0.05 * (maximum(λvec) - minimum(λvec))
+end
 
 #------------------
 # HOMOGENEOUS CASE
 #------------------
 
-function rand_single(rng::Random.AbstractRNG, p::PoissonProcess{<:Real}, g::Geometry, ::DiscretizedSampling)
+function rand_single(rng::Random.AbstractRNG, p::PoissonProcess{<:Real}, g, ::DiscretizedSampling)
   # simulate number of points
   λ = p.λ
   V = measure(g)
   n = rand(rng, Poisson(λ * V))
 
-  # simulate homogeneous process
-  pts = sample(g, HomogeneousSampling(n))
-
-  # return point pattern
-  PointSet(collect(pts))
-end
-
-function rand_single(rng::Random.AbstractRNG, p::PoissonProcess{<:Real}, b::Box, ::ProductSampling)
-  # region configuration
-  lo, up = coordinates.(extrema(b))
-
-  # simulate number of points
-  λ = p.λ
-  V = measure(b)
-  n = rand(rng, Poisson(λ * V))
-
-  # product of uniform distributions
-  U = product_distribution([Uniform(lo[i], up[i]) for i in 1:embeddim(b)])
+  if iszero(n)
+    # PointSet{embeddim(g), coordtype(g)}([])
+    nothing
+  else
+    # simulate homogeneous process
+    points = sample(g, HomogeneousSampling(n))
 
-  # return point pattern
-  PointSet(rand(rng, U, n))
+    # return point pattern
+    PointSet(points)
+  end
 end
 
 #--------------------
 # INHOMOGENEOUS CASE
 #--------------------
 
-function rand_single(rng::Random.AbstractRNG, p::PoissonProcess{<:Function}, b::Box, algo::DiscretizedSampling)
-  # discretize region
-  # TODO
+function rand_single(rng::Random.AbstractRNG, p::PoissonProcess{<:AbstractVector}, d::Domain, algo::DiscretizedSampling)
+  # simulate number of points
+  λ = p.λ
+  V = measure.(d)
+  n = rand(rng, Poisson(sum(λ .* V)))
+
+  # simulate n points
+  if iszero(n)
+    nothing
+  else
+    # sample elements with weights proportial to expected number of points
+    w = WeightedSampling(n, λ .* V, replace = true)
+
+    # within each element sample a single point
+    sampler = HomogeneousSampling(1)
+    points = (first(sample(rng, e, sampler)) for e in sample(rng, d, w))
+
+    # return point pattern
+    PointSet(points)
+  end
+end
 
-  # discretize retention
-  # TODO
+function rand_single(rng::Random.AbstractRNG, p::PoissonProcess{<:Function}, g, algo::ThinnedSampling)
+  # simulate a homogeneous process
+  pp = rand_single(rng, PoissonProcess(algo.λmax), g, DiscretizedSampling())
 
-  # sample each element
-  # TODO
+  # thin point pattern
+  thin(pp, RandomThinning(x -> p.λ(x) / algo.λmax))
 end
 
-function rand_single(rng::Random.AbstractRNG, p::PoissonProcess{<:Function}, b::Box, algo::ThinnedSampling)
-  # Lewis-Shedler algorithm
-  # TODO
+function rand_single(rng::Random.AbstractRNG, p::PoissonProcess{<:Function}, d::Domain, algo::DiscretizedSampling)
+  # compute intensity on centroids
+  c = centroid.(d)
+  λvec = p.λ.(c)
+
+  # simulate inhomogeneous process
+  rand_single(rng, PoissonProcess(λvec), d, algo)
 end

src/processes/union.jl

@@ -21,9 +21,9 @@ Base.union(p₁::PointProcess, p₂::PointProcess) = UnionProcess(p₁, p₂)
 
 ishomogeneous(p::UnionProcess) = ishomogeneous(p.p₁) && ishomogeneous(p.p₂)
 
-default_sampling_algorithm(::UnionProcess, ::Geometry) = UnionSampling()
+default_sampling_algorithm(::UnionProcess, ::Any) = UnionSampling()
 
-function rand_single(rng::Random.AbstractRNG, p::UnionProcess, g::Geometry, ::UnionSampling)
+function rand_single(rng::Random.AbstractRNG, p::UnionProcess, g, ::UnionSampling)
   pp₁ = rand(rng, p.p₁, g)
   pp₂ = rand(rng, p.p₂, g)
 

test/processes.jl

@@ -1,4 +1,18 @@
 @testset "Processes" begin
+  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)]
+  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)
+
   @testset "Basic" begin
     for p in [BinomialProcess(100), PoissonProcess(100.0)]
       b = Box((0.0, 1.0), (1.0, 2.0))
@@ -10,22 +24,54 @@
   end
 
   @testset "Binomial" begin
-    # TODO
+    p = BinomialProcess(10)
+    for g in [seg, tri, quad, box, ball, poly, grid, mesh]
+      pp = rand(p, g)
+      @test nelements(pp) == 10
+      @test all(∈(g), pp)
+    end
   end
 
   @testset "Poisson" begin
-    p = PoissonProcess(100.0)
-    t = Triangle((0.0, 0.0), (1.0, 0.0), (1.0, 1.0))
-    pp = rand(p, t)
-    xs = coordinates.(pp)
-    @test all(0 .≤ first.(xs) .≤ 1)
-    @test all(0 .≤ last.(xs))
-
-    q = Quadrangle((0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0))
-    pp = rand(p, q)
-    xs = coordinates.(pp)
-    @test all(0 .≤ first.(xs) .≤ 1)
-    @test all(0 .≤ last.(xs) .≤ 1)
+    # 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
+    λ(s::Point2) = sum(coordinates(s) .^ 2)
+    p = PoissonProcess(λ)
+    for g in [seg, tri, quad, box, ball, poly, grid, mesh]
+      # default thinnedsampling
+      pp = rand(p, g)
+      @test all(∈(g), pp)
+      # custom thinnedsampling using λmax
+      pp = rand(p, g, algo = PointPatterns.ThinnedSampling(λ(Point2(12.0, 12.0))))
+      @test all(∈(g), pp)
+      # discretizedsampling
+      g = discretize(g)
+      pp = rand(p, g, algo = PointPatterns.DiscretizedSampling())
+      @test all(∈(g), g)
+    end
+
+    # inhomogeneous with piecewise constant intensity
+    for g in [grid, mesh]
+      λ(s::Point2) = sum(coordinates(s) .^ 2)
+      λ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]
+      @test isnothing(rand(PoissonProcess(0.0), seg))
+    end
+    ps = PointSet(rand(Point2, 10))
+    @test isnothing(rand(PoissonProcess(100.0), ps))
   end
 
   @testset "Union" begin