Fix poisson in union of disk for hard core spatial

Author: guilgautier

src/display.jl

@@ -19,7 +19,7 @@ function plot(
         Plots.plot!(x[1] .+ circ_x,
                     x[2] .+ circ_y,
                     color="black",
-                    linewidth=0.2)
+                    linewidth=0.5)
     end
 
     win = pp.window
@@ -33,7 +33,6 @@ function plot(
         ising::Ising,
         state
 )
-
     pos = collect(Iterators.product(1:ising.dims[1], 1:ising.dims[2]))[:]
     locs_x, locs_y = map(x->x[1], pos), map(x->x[2], pos)
 

src/hard_core_spatial.jl

@@ -67,6 +67,9 @@ end
 
 """
 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};
@@ -115,32 +118,33 @@ function generate_sample_poisson_union_disks(
     for (i, c) in enumerate(eachcol(centers))
         n = rand(rng, 𝒫)
         n == 0 && continue
-        proposed_points = generate_sample_uniform_in_disk(n, r, c; rng=rng)
+        proposed_points = generate_sample_uniform_in_disk(n, c, r; rng=rng)
         if i == 1
             points = hcat(points, proposed_points)
             continue
         end
-        centers_1i = @view centers[:, 1:i]
-        accept = vec(all(pairwise_distances(centers_1i, proposed_points) .> r, dims=1))
+        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; rng=-1)
+    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,
-        radius::Real,
-        center::AbstractVector;
+        center::AbstractVector,
+        radius::Real;
         rng=-1
 )::Matrix{Float64}
-    @assert n > 0
-    @assert length(center) == 2
     @assert radius > 0
+    @assert length(center) == 2
+
+    n == 0 && return Matrix{Float64}(undef, 2, 0)
     rng = getRNG(rng)
 
     sample = repeat(center, 1, n)
@@ -152,22 +156,6 @@ function generate_sample_uniform_in_disk(
     return sample
 end
 
-"""
-    pairwise_distances(X, Y)
-
-Pairwise distance matrix between columns of `X` and `Y`.
-Equivalent to `[norm(x - y) for x in X, y in Y]`.
-"""
-function pairwise_distances(X, Y)
-    return Distances.pairwise(Distances.Euclidean(1e-8), X, Y; dims=2)
-end
-
-function pairwise_distances(X)
-    dist = Distances.pairwise(Distances.Euclidean(1e-8), X; dims=2)
-    dist[LA.diagind(dist)] .= Inf
-    return dist
-end
-
 # Default sample generate_sample
 
 """

src/utils.jl

@@ -55,3 +55,20 @@ function random_edge_orientation(
     rng = getRNG(rng)
     return LG.SimpleDiGraphFromIterator(rand(rng) < p ? e : reverse(e) for e in LG.edges(g))
 end
+
+
+"""
+    pairwise_distances(X, Y)
+
+Pairwise distance matrix between columns of `X` and `Y`.
+Equivalent to `[norm(x - y) for x in X, y in Y]`.
+"""
+function pairwise_distances(X, Y)
+    return Distances.pairwise(Distances.Euclidean(1e-8), X, Y; dims=2)
+end
+
+function pairwise_distances(X)
+    dist = Distances.pairwise(Distances.Euclidean(1e-8), X; dims=2)
+    dist[LA.diagind(dist)] .= Inf
+    return dist
+end

src/window.jl

@@ -1,44 +1,54 @@
 abstract type AbstractWindow end
 
+# Discrete
+
 abstract type AbstractDiscreteWindow{T} <: AbstractWindow end
 
 struct GraphNode{T<:Int64} <: AbstractDiscreteWindow{T}
     # Corner
     idx::T
 end
+function GraphNode(idx::Integer)
+    return GraphNode{typeof(idx)}(idx)
+end
+
+# Spatial
 
 abstract type AbstractSpatialWindow{T<:Float64} <: AbstractWindow end
+
+abstract type AbstractRectangleWindow{T} <: AbstractSpatialWindow{T} end
 # Πᵢ [cᵢ, cᵢ + wᵢ]
-struct RectangleWindow{T} <: AbstractSpatialWindow{T}
+struct RectangleWindow{T} <: AbstractRectangleWindow{T}
     # Corner
     c::Vector{T}
     # Width
     w::Vector{T}
 end
 
+function RectangleWindow(c::Vector, w::Vector)
+    @assert length(c) == length(w)
+    return RectangleWindow{Float64}(c, w)
+end
+
 # Πᵢ [cᵢ, cᵢ + w]
-struct SquareWindow{T} <: AbstractSpatialWindow{T}
+struct SquareWindow{T} <: AbstractRectangleWindow{T}
     # Corner
     c::Vector{T}
     # Width
     w::T
 end
 
-function GraphNode(idx::Integer)
-    return GraphNode{typeof(idx)}(idx)
-end
-
-function RectangleWindow(c::Vector, w::Vector)
-    @assert length(c) == length(w)
-    return RectangleWindow{Float64}(float.(c), float.(w))
+function SquareWindow(c::Vector, w::Real)
+    return SquareWindow{Float64}(c, w)
 end
 
-function SquareWindow(c::Vector, w::Real)
-    return SquareWindow{Float64}(float.(c), float(w))
+struct BallWindow{T} <: AbstractSpatialWindow{T}
+    center::Vector{T}
+    radius::T
 end
 
-function spatial_window(c::Vector, w::Union{Real,Vector})
-    return allequal(w) ? SquareWindow(c, w[1]) : RectangleWindow(c, w)
+function BallWindow(center::Vector, radius::Real)
+    return BallWindow{Float64}(center, radius)
 end
 
 dimension(win::GraphNode) = 0
@@ -47,11 +57,47 @@ dimension(win::AbstractSpatialWindow) = length(win.c)
 volume(win::GraphNode) = 0
 volume(win::RectangleWindow) = prod(win.w)
 volume(win::SquareWindow) = win.w^dimension(win)
+function volume(win::BallWindow)
+    d = dimension(win)
+    return π^(d/2) * win.radius^d / gamma(d/2 + 1)
+end
+
+function Base.in(
+    x::AbstractVector,
+    win::AbstractRectangleWindow
+)
+    return all(win.c .<= x) && all(x .<= win.c .+ win.w)
+end
+
+function Base.in(
+    x::AbstractVector,
+    win::BallWindow
+)
+    return Distances.euclidean(x, win.center) <= win.radius
+end
 
+"""
+Sample uniformly in `RectangleWindow`
+"""
 function Base.rand(
-    win::AbstractSpatialWindow{T};
+    win::AbstractRectangleWindow{T};
     rng=-1
 )::Vector{T} where {T}
     rng = getRNG(rng)
     return win.c .+ win.w .* rand(rng, T, dimension(win))
 end
+
+"""
+Sample uniformly in `BallWindow`
+"""
+function Base.rand(
+    win::BallWindow{T};
+    rng=-1
+)::Vector{T} where {T}
+    rng = getRNG(rng)
+    d = dimension(win)
+    x = zeros(T, d+2)
+    randn!(rng, x)
+    normalize!(x)
+    return win.radius .* x[1:d] .+ win.center
+end

test/hard_core_spatial.jl

@@ -2,7 +2,7 @@ using PartialRejectionSampling
 const PRS = PartialRejectionSampling
 
 β₀ = 0.1
-r = 0.01  # interaction range = 2*radius
+r = 0.05  # interaction range = 2*radius
 b = β₀ / (π * (r/2)^2)
 
 c, w = [0.0, 0.0], 1.0