Set warnings when using prs with non provable guaranty for hard core graph and spatial

Author: guilgautier

src/graph/hard_core.jl

@@ -62,7 +62,7 @@ function HardCoreGraph(
     β::Real
 ) where {T<:Integer}
     @assert β >= 0
-    return HardCoreGraph{T}(graph, float(β))
+    return HardCoreGraph{T}(graph, β)
 end
 
 # Sampling
@@ -114,6 +114,12 @@ function generate_sample_prs(
     pp::HardCoreGraph{T}
 )::Vector{T} where {T}
 
+    β_max = 1 / (2 * sqrt(exp(1)) * LG.Δ(pp.graph) - 1)
+    if pp.β > β_max
+        @warn "The arguments do not satisfy Theorem 35 of Guo, Jerrum and Liu (2019): partial rejection sampling may not be efficient.
+        Given `graph`, consider choosing β ≤ $(β_max)"
+    end
+
     proba = pp.β / (one(pp.β) + pp.β)
 
     adj = LG.adjacency_matrix(pp.graph)

src/grid_prs.jl

@@ -240,10 +240,10 @@ function weighted_interaction_graph(
 )::SWG.SimpleWeightedGraph
     window_ = window(spp)
     if allequal(window_.w)
-        nb_cells_x = ceil(Int, window_.w[1] / spp.r)
+        nb_cells_x = cld(window_.w[1], spp.r)
         return uniform_weighted_graph(rng, king_graph(nb_cells_x))
     else
-        nb_cells_x, nb_cells_y = ceil.(Int, window_.w ./ spp.r)
+        nb_cells_x, nb_cells_y = cld.(Int, window_.w, spp.r)
         return uniform_weighted_graph(rng, king_graph(nb_cells_x, nb_cells_y))
     end
 end
@@ -264,9 +264,9 @@ function initialize_cells(
     window_ = window(spp)
     nb_cells_x = nb_cells_y = 1
     if allequal(window_.w)  # SquareWindow
-        nb_cells_x = nb_cells_y = ceil(Int, window_.w[1] / spp.r)
+        nb_cells_x = nb_cells_y = cld(window_.w[1], spp.r)
     else
-        nb_cells_x, nb_cells_y = ceil.(Int, window_.w ./ spp.r)
+        nb_cells_x, nb_cells_y = cld.(Int, window_.w, spp.r)
     end
 
     cells = Vector{SpatialCellGridPRS{T}}(undef, nb_cells_x * nb_cells_y)

src/spatial/hard_core.jl

@@ -197,9 +197,16 @@ function generate_sample_prs(
     pp::HardCorePointProcess{T},
     win::Union{Nothing,AbstractWindow}=nothing
 )::Vector{T} where {T}
-    win_ = win === nothing ? window(pp) : win
 
     β, r = intensity(pp), interaction_range(pp)
+    win_ = win === nothing ? window(pp) : win
+
+    β_max = 0.21027 / volume(BallWindow(zeros(dimension(win_)), r/2))
+    if β > β_max
+        @warn "The arguments do not satisfy Lemma 6 of Guo and Jerrum (2018): partial rejection sampling may not be efficient.
+        Given the interaction range `pp.r`, consider choosing β ≤ $(β_max)"
+    end
+
     points = generate_sample(rng, HomogeneousPoissonPointProcess(β, win_))
     while true
         bad = vec(any(pairwise_distances(points) .<= r, dims=2))