ClusterProcess (#36)

* start with PoissonClusterProcess

* Update src/processes/clustered.jl

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

* Apply suggestions from code review

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

* Update src/processes/clustered.jl

* Update src/processes/clustered.jl

* Rename to ClusterProcess

* Update src/processes/cluster.jl

* more general cluster process

* sampling algorithm

* fix pointset

* Apply suggestions from code review

* Update src/processes/cluster.jl

* Update src/processes/cluster.jl

* Update src/processes/cluster.jl

* Refactor cluster.jl

---------

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

Author: ErickChacon

src/PointPatterns.jl

@@ -19,6 +19,7 @@ export
   PointProcess,
   BinomialProcess,
   PoissonProcess,
+  ClusterProcess,
   UnionProcess,
   ishomogeneous,
 

src/processes.jl

@@ -85,6 +85,7 @@ struct ConstantIntensity <: PointPatternAlgo end
 
 include("processes/binomial.jl")
 include("processes/poisson.jl")
+include("processes/cluster.jl")
 
 # ----------------
 # UNION PROCESSES

src/processes/cluster.jl

@@ -0,0 +1,53 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+"""
+    ClusterProcess(p, ofun)
+
+A cluster process with parent process `p` and offsprings generated
+with `ofun`. It is a function that takes a parent point and returns
+a point pattern from another point process.
+
+    ClusterProcess(p, o, gfun)
+
+Alternatively, specify the parent process `p`, the offspring process `o`
+and the geometry function `gfun`. It is a function that takes a parent
+point and returns a geometry or domain for the offspring process.
+
+## Examples
+
+```julia
+# TODO
+```
+"""
+struct ClusterProcess{P<:PointProcess,F<:Function} <: PointProcess
+  p::P
+  ofun::F
+end
+
+ClusterProcess(p::PointProcess, o::PointProcess, gfun::Function) =
+  ClusterProcess(p, parent -> rand(o, gfun(parent)))
+
+default_sampling_algorithm(p::ClusterProcess, g) = nothing
+
+function rand_single(rng::Random.AbstractRNG, p::ClusterProcess, g, ::Nothing)
+  # retrieve parameters
+  Dim = embeddim(g)
+  T = coordtype(g)
+
+  # generate parents
+  parents = rand_single(rng, p.p, g)
+
+  # generate offsprings
+  offsprings = p.ofun.(parents)
+
+  # combine offsprings into single set
+  points = mapreduce(vcat, offsprings) do pset
+    isnothing(pset) ? Point{Dim,T}[] : collect(view(pset, g))
+  end
+
+  # return point pattern
+  PointSet(points)
+end
+

src/processes/poisson.jl

@@ -65,8 +65,7 @@ function rand_single(rng::Random.AbstractRNG, p::PoissonProcess{<:Function}, g,
   pp = rand_single(rng, PoissonProcess(algo.λmax), g, ConstantIntensity())
 
   # thin point pattern
-  points = thin(pp, RandomThinning(x -> p.λ(x) / algo.λmax))
-  PointSet(collect(points))
+  isnothing(pp) ? nothing : PointSet(collect(thin(pp, RandomThinning(x -> p.λ(x) / algo.λmax))))
 end
 
 rand_single(rng::Random.AbstractRNG, p::PoissonProcess{<:Function}, d::Domain, algo::ConstantIntensity) =