Restructure and fix hard_core_spatial

Author: guilgautier

src/PartialRejectionSampling.jl

@@ -31,6 +31,7 @@ include("dominated_cftp.jl")
 include("grid_prs.jl")
 
 # Spatial point processes
+include("poisson.jl")
 include("strauss.jl")
 include("hard_core_spatial.jl")
 

src/dominated_cftp.jl

@@ -3,35 +3,35 @@
 #  - Huber, Perfect Simulation
 #  - [Kendall's notes on perfect simulation](https://warwick.ac.uk/fac/sci/statistics/staff/  academic-research/kendall/personal/ppt/428.pdf)
 #
-# Requirement: the target spatial point process must have the followingmethods
+# Requirement: the target spatial point process must have the following methods
 #  - function papangelou_conditional_intensity end
 #  - function upper_bound_papangelou_conditional_intensity end
 #  - function window end
 
 function generate_sample_dcftp(
         pp::AbstractSpatialPointProcess{T};
-        n₀::Int64=1,
+        n₀::Int=1,
         win::Union{Nothing,AbstractWindow}=nothing,
         rng=-1
 )::Vector{T} where {T}
 
     @assert n₀ >= 1
     rng = getRNG(rng)
 
-    win_ = win === nothing ? window(pp) : win
+    window_ = win === nothing ? window(pp) : win
     β = upper_bound_papangelou_conditional_intensity(pp)
-    birth_rate = β * volume(win_)
+    birth_rate = β * volume(window_)
 
     # Dominating process
-    k = rand(rng, Distributions.Poisson(birth_rate))
-    D = Set{T}(rand(win_; rng=rng) for _ in 1:k)
+    hp = HomogeneousPoissonPointProcess(β, window_)
+    D = Set{T}(eachcol(generate_sample(hp; rng=rng)))
 
     M = Float64[]  # Marking process
     R = T[]        # Recording process
 
     steps = -1:-1:-n₀
     while true
-        backward_update!(D, M, R, steps, birth_rate, win_; rng=rng)
+        backward_update!(D, M, R, steps, birth_rate, window_; rng=rng)
         coupling, L = forward_coupling(D, M, R, pp, β)
         coupling && return collect(L)
         steps = (steps.stop-1):-1:(2*steps.stop)
@@ -44,7 +44,7 @@ function backward_update!(
         R::Vector{T},       # Recording process
         steps::StepRange,   # Number of backward steps
         birth_rate::Real,
-        window::AbstractWindow;
+        win::AbstractWindow;
         rng=-1
 ) where {T}
     rng = getRNG(rng)
@@ -58,7 +58,7 @@ function backward_update!(
             pushfirst!(M, rand(rng))
         else
             # forward birth (push) ≡ backward death (pushfirst)
-            x = rand(window; rng)
+            x = rand(win; rng=rng)
             push!(D, x)
             pushfirst!(R, x)
             pushfirst!(M, 0.0)
@@ -69,9 +69,9 @@ end
 function forward_coupling(
         D::Set{T},           # Dominating process
         M::Vector{Float64},  # Marking process
-        R::Vector{T},         # Recording process
+        R::Vector{T},        # Recording process
         pp::AbstractSpatialPointProcess{T},
-        β::Real             # Upper bound on papangelou conditional intensity
+        β::Real              # Upper bound on papangelou conditional intensity
 ) where {T}
     # L ⊆ X ⊆ U ⊆ D, where X is the target process
     L, U = empty(D), copy(D)

src/hard_core_spatial.jl

@@ -30,9 +30,6 @@ function HardCorePointProcess(β, r, c, w)
     return HardCorePointProcess(β, r, spatial_window(c, w))
 end
 
-window(hc::HardCorePointProcess) = hc.window
-dimension(hc::HardCorePointProcess) = dimension(window(hc))
-
 ## Sampling
 
 # Used in dominated CFTP, see dominated_cftp.jl
@@ -62,14 +59,10 @@ function gibbs_interaction(
     return !any(any(Distances.euclidean(x, y) <= hc.r for x in c_i) for y in c_j)
 end
 
-# Partial Rejection Sampling
-# The implementation relies on Distances.pairwise, points are stored as columns of a Matrix.
+# Partial Rejection sampling
 
 """
 Sample from spatial hard core model a.k.a Poisson disk sampling with intensity λ and radius r on [0, 1]^2, using Partial Rejection Sampling [https://arxiv.org/pdf/1801.07342.pdf](H. Guo, M. Jerrum).
-!!! warning "Side effects"
-
-    In the current implementation, the return sample may contain some point outside of the box! (see generate_sample_poisson_union_disks)
 """
 function generate_sample_prs(
         hc::HardCorePointProcess{T};
@@ -80,83 +73,20 @@ function generate_sample_prs(
     rng = getRNG(rng)
     window_ = win === nothing ? window(hc) : win
 
-    n = rand(rng, Distributions.Poisson(hc.β * volume(window_)))
-    points = Matrix{Float64}(undef, dimension(hc), n)
-    for x in eachcol(points)
-        x .= rand(window_; rng=rng)
-    end
+    points = generate_sample(HomogeneousPoissonPointProcess(hc.β, window_); rng=rng)
 
     while true
         bad = vec(any(pairwise_distances(points) .< hc.r, dims=2))
         !any(bad) && break
 
-        resampled = generate_sample_poisson_union_disks(hc.β, hc.r, points[:, bad]; rng=rng)
+        resampled = generate_sample_poisson_union_balls(hc.β, points[:, bad], hc.r;
+                                                        win=window_, rng=rng)
         points = hcat(points[:, .!bad], resampled)
     end
-    return [points[:, i] for i in 1:size(points, 2)]
-end
-
-"""
-Sample from Poisson(λ) on ⋃ᵢ B(cᵢ, r) (union of disks of radius r centered at cᵢ)
-
-Use the independence property of the Poisson point process
-
-  - Sample from Poisson(λ) on B(c₁, r)
-  - Sample from Poisson(λ) on B(c₂, r) ∖ B(c₁, r)
-  - Sample from Poisson(λ) on B(cⱼ, r) ∖ ⋃_i<j B(cᵢ, r)
-  - ...
-"""
-function generate_sample_poisson_union_disks(
-        λ::Real,
-        r::Real,
-        centers::Matrix;
-        rng=-1
-)::Matrix{Float64}
-    rng = getRNG(rng)
-    𝒫 = Distributions.Poisson(λ * π * r^2)
-    points = Matrix{Float64}(undef, 2, 0)
-    for (i, c) in enumerate(eachcol(centers))
-        n = rand(rng, 𝒫)
-        n == 0 && continue
-        proposed_points = generate_sample_uniform_in_disk(n, c, r; rng=rng)
-        if i == 1
-            points = hcat(points, proposed_points)
-            continue
-        end
-        centers_ = @view centers[:, 1:i-1]
-        accept = vec(all(pairwise_distances(centers_, proposed_points) .> r, dims=1))
-        points = hcat(points, proposed_points[:, accept])
-    end
-    return points
-end
-
-"""
-    generate_sample_uniform_in_disk(n, center, radius,
-
-Sample `n` points uniformly in ``D(c, r)`` (disk with `center` and `radius`)
-"""
-function generate_sample_uniform_in_disk(
-        n::Integer,
-        center::AbstractVector,
-        radius::Real;
-        rng=-1
-)::Matrix{Float64}
-    @assert radius > 0
-    @assert length(center) == 2
-
-    n == 0 && return Matrix{Float64}(undef, 2, 0)
-    rng = getRNG(rng)
-
-    sample = repeat(center, 1, n)
-    for x in eachcol(sample)
-        r = radius * sqrt(rand(rng))
-        theta = 2π * rand(rng)
-        x .+= [r * cos(theta), r * sin(theta)]
-    end
-    return sample
+    return collect(T, eachcol(points))
 end
 
-# Default sample generate_sample
+# Default sampler generate_sample
 
 """
 Default sampler for HardCorePointProcess

src/poisson.jl

@@ -0,0 +1,68 @@
+struct HomogeneousPoissonPointProcess{T} <: AbstractSpatialPointProcess{T}
+    β::Float64
+    window::AbstractSpatialWindow{Float64}
+end
+
+function HomogeneousPoissonPointProcess(β::Real, window::AbstractSpatialWindow)
+    @assert β > 0
+    return HomogeneousPoissonPointProcess{Vector{Float64}}(β, window)
+end
+
+## Sampling
+
+"""
+Sample from a homogenous Poisson point process `hp` on window `win` for which the `volume` function is defined.
+"""
+function generate_sample(
+        hp::HomogeneousPoissonPointProcess;
+        win::Union{Nothing,AbstractWindow}=nothing,
+        rng=-1
+)::Matrix{Float64}
+    rng = getRNG(rng)
+    win_ = win === nothing ? window(hp) : win
+    n = rand(rng, Distributions.Poisson(hp.β * volume(win_)))
+    return rand(win_, n; rng=rng)
+end
+
+"""
+Sample from a homogenous Poisson point process Poisson(β) on ⋃ᵢ B(cᵢ, r) (union of balls centered at cᵢ with the same radius r)
+
+Use the independence property of the Poisson point process on disjoint subsets in order to
+
+  - Sample from Poisson(β) on B(c₁, r)
+  - Sample from Poisson(β) on B(c₂, r) ∖ B(c₁, r)
+  - Sample from Poisson(β) on B(cⱼ, r) ∖ ⋃_i<j B(cᵢ, r)
+  - ...
+"""
+function generate_sample_poisson_union_balls(
+        β::Real,
+        centers::Matrix,
+        radius::Real;
+        win::Union{Nothing,AbstractWindow}=nothing,
+        rng=-1
+)::Matrix{Float64}
+    rng = getRNG(rng)
+    d = size(centers, 1)
+    !(win === nothing) && @assert dimension(win) == d
+
+    𝒫 = Distributions.Poisson(β * volume(BallWindow(zeros(d), radius)))
+    points = Matrix{Float64}(undef, d, 0)
+    for (i, c) in enumerate(eachcol(centers))
+        n = rand(rng, 𝒫)
+        n == 0 && continue
+        proposed = rand(BallWindow(c, radius), n; rng=rng)
+        if i == 1
+            points = hcat(points, proposed)
+            continue
+        end
+        centers_ = @view centers[:, 1:i-1]
+        accept = vec(all(pairwise_distances(centers_, proposed) .> radius, dims=1))
+        points = hcat(points, proposed[:, accept])
+    end
+
+    if win === nothing || isempty(points)
+        return points
+    else
+        return points[:, vec(mapslices(x -> x in win, points; dims=1))]
+    end
+end

src/strauss.jl

@@ -23,16 +23,9 @@ function StraussPointProcess(β::Real, γ::Real, r::Real, window::AbstractSpatia
     @assert β > 0
     @assert 0 <= γ <= 1
     @assert r > 0
-    return StraussPointProcess{Vector{Float64}}(float(β), float(γ), float(r), window)
+    return StraussPointProcess{Vector{Float64}}(β, γ, r, window)
 end
 
-function StraussPointProcess(β, γ, r, c, w)
-    return StraussPointProcess(β, γ, r, spatial_window(c, w))
-end
-
-window(strauss::StraussPointProcess) = strauss.window
-dimension(strauss::StraussPointProcess) = dimension(window(strauss))
-
 ## Sampling
 
 # Used in dominated_cftp

src/utils.jl

@@ -1,4 +1,4 @@
-getRNG(seed::Integer = -1) = seed >= 0 ? Random.MersenneTwister(seed) : Random.GLOBAL_RNG
+getRNG(seed::Integer=-1) = seed >= 0 ? Random.MersenneTwister(seed) : Random.GLOBAL_RNG
 getRNG(seed::Union{Random.MersenneTwister,Random._GLOBAL_RNG}) = seed
 
 sigmoid(x) = @. 1 / (1 + exp(-x))
@@ -13,6 +13,13 @@ sigmoid(x) = @. 1 / (1 + exp(-x))
     return true
 end
 
+function normalize_columns!(X::Matrix, p::Real=2)
+    for x in eachcol(X)
+        LA.normalize!(x, p)
+    end
+    return X
+end
+
 ## Graph functions
 
 edgemap(g::LG.AbstractGraph) = Dict(LG.edges(g) .=> 1:LG.ne(g))