Restructure backend grid prs to work on square rectangle

Author: guilgautier

src/display.jl

@@ -69,7 +69,8 @@ function plot(
         locs_x, locs_y = map(x->x[1], pos), map(x->x[2], pos)
 
         p = GraphPlot.gplot(hcg.g,
-                locs_x, locs_y,
+                locs_x,
+                reverse(locs_y),
                 nodelabel=LG.vertices(hcg.g),
                 nodefillc=col_nodes
                 )
@@ -99,7 +100,8 @@ function plot(
         locs_x, locs_y = map(x->x[1], pos), map(x->x[2], pos)
 
         p = GraphPlot.gplot(sample,
-            locs_x, locs_y,
+            locs_x,
+            reverse(locs_y),
             nodelabel=LG.vertices(sample),
             nodefillc=col_nodes,
             edgestrokec=col_edges,
@@ -129,7 +131,8 @@ function plot(
         locs_x, locs_y = map(x->x[1], pos), map(x->x[2], pos)
 
         p = GraphPlot.gplot(sample,
-            locs_x, locs_y,
+            locs_x,
+            reverse(locs_y),
             nodelabel=LG.vertices(sample),
             nodefillc=col_nodes,
             arrowlengthfrac=0.05,

src/grid_prs.jl

@@ -25,9 +25,9 @@ function generate_sample_grid_prs(
         rng=-1
 )::Vector{T} where {T}
     rng = getRNG(rng)
-    g = weighted_interaction_graph(pp; rng=rng)
 
-    cells = initialize_cells(pp, LG.nv(g))
+    g = weighted_interaction_graph(pp; rng=rng)
+    cells = initialize_cells(pp)
     resample_indices = Set(1:length(cells))
 
     while !isempty(resample_indices)
@@ -47,7 +47,17 @@ function weighted_interaction_graph(
         rng=-1
 )::SWG.SimpleWeightedGraph
     rng = getRNG(rng)
-    g = SWG.SimpleWeightedGraph(king_graph(ceil(Int, inv(pp.r))))
+
+    window_ = window(pp)
+    g = SWG.SimpleWeightedGraph()
+    nb_cells_x = nb_cells_y = 1
+    if allequal(window_.w)
+        nb_cells_x = ceil(Int, window_.w[1] / pp.r)
+        g = SWG.SimpleWeightedGraph(king_graph(nb_cells_x))
+    else
+        nb_cells_x, nb_cells_y = ceil.(Int, window_.w ./ pp.r)
+        g = SWG.SimpleWeightedGraph(king_graph(nb_cells_x, nb_cells_y))
+    end
     for e in LG.edges(g)
         i, j = Tuple(e)
         @inbounds g.weights[i, j] = g.weights[j, i] = rand(rng, weighttype(g))
@@ -144,14 +154,21 @@ end
 
 function initialize_cells(
         spp::AbstractSpatialPointProcess{T},
-        size::Integer
 )::Vector{SpatialCellGridPRS{T}} where {T}
-    cells = Vector{SpatialCellGridPRS{T}}(undef, size)
-    k = ceil(Int, inv(spp.r))
+    dimension(spp) != 2 && throw(DomainError(spp, "must be 2D point process for now"))
+    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)
+    else
+        nb_cells_x, nb_cells_y = ceil.(Int, window_.w ./ spp.r)
+    end
+
+    cells = Vector{SpatialCellGridPRS{T}}(undef, nb_cells_x * nb_cells_y)
     for i in eachindex(cells)
-        c_y, c_x = divrem(i-1, k)
-        c = spp.window.c + spp.r .* [c_x, c_y]
-        win_i = rectangle_square_window(c, min.(spp.r, spp.window.w .- c))
+        c_ij = divrem(i-1, nb_cells_x)  # coordinate (i, j) of the cell in the grid
+        c_xy = @. window_.c + spp.r * c_ij
+        win_i = rectangle_square_window(c_xy, @. min(spp.r, window_.w - c_xy))
         cells[i] = SpatialCellGridPRS(win_i, T[])
     end
     return cells

src/ising/ising_sampling_grid_prs.jl

@@ -14,10 +14,9 @@ function weighted_interaction_graph(
 end
 
 function initialize_cells(
-        ising::Ising{T},
-        size::Integer
+        ising::Ising{T}
 )::Vector{GraphCellGridPRS{T}} where {T}
-    return [GraphCellGridPRS(GraphNode(i), zero(T)) for i in 1:size]
+    return [GraphCellGridPRS(GraphNode(i), zero(T)) for i in 1:LG.nv(ising.g)]
 end
 
 function generate_sample!(

src/poisson.jl

@@ -19,9 +19,9 @@ function generate_sample(
         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)
+    window_ = win === nothing ? window(hp) : win
+    n = rand(rng, Distributions.Poisson(hp.β * volume(window_)))
+    return rand(window_, n; rng=rng)
 end
 
 """

src/utils.jl

@@ -25,13 +25,35 @@ end
 edgemap(g::LG.AbstractGraph) = Dict(LG.edges(g) .=> 1:LG.ne(g))
 sink_nodes(g::LG.SimpleDiGraph) = [v for v in LG.vertices(g) if LG.outdegree(g, v) == 0]
 
-
 """
-Return the [strong product](https://en.wikipedia.org/wiki/Strong_product_of_graphs) of a graph with itself
+Return the [strong product](https://en.wikipedia.org/wiki/Strong_product_of_graphs) of two graphs
 
 See also:
 Section 5 of [alternative text](https://drops.dagstuhl.de/opus/volltexte/2019/11538/pdf/LIPIcs-ISAAC-2019-42.pdf)
 """
+function strong_product(g::G, h::G)::G where {G<:LG.AbstractGraph}
+    sp_gh = G(LG.nv(g) * LG.nv(h))
+    id(i, j) = (i - 1) * LG.nv(h) + j
+    for e1 in LG.edges(g)
+        i1, i2 = Tuple(e1)
+        for j = 1:LG.nv(h)
+            LG.add_edge!(sp_gh, id(i1, j), id(i2, j))
+        end
+        for e2 in LG.edges(h)
+            j1, j2 = Tuple(e2)
+            LG.add_edge!(sp_gh, id(i1, j1), id(i2, j2))
+            LG.add_edge!(sp_gh, id(i1, j2), id(i2, j1))
+        end
+    end
+    for e in LG.edges(h)
+        j1, j2 = Tuple(e)
+        for i in LG.vertices(g)
+            LG.add_edge!(sp_gh, id(i, j1), id(i, j2))
+        end
+    end
+    return sp_gh
+end
+
 function strong_product(g::G)::G where {G<:LG.AbstractGraph}
     n = LG.nv(g)
     sp_gg = G(n * n)  # nv = n^2, ne = 2n(2n+1)
@@ -51,6 +73,7 @@ function strong_product(g::G)::G where {G<:LG.AbstractGraph}
     return sp_gg
 end
 
+king_graph(k::Integer, l::Integer) = strong_product(LG.path_graph(l), LG.path_graph(k))
 king_graph(k::Integer) = strong_product(LG.path_graph(k))
 
 function random_edge_orientation(
@@ -63,7 +86,6 @@ function random_edge_orientation(
     return LG.SimpleDiGraphFromIterator(rand(rng) < p ? e : reverse(e) for e in LG.edges(g))
 end
 
-
 """
     pairwise_distances(X, Y)
 

src/window.jl

@@ -44,8 +44,8 @@ function SquareWindow(c::AbstractVector, w::Real=1.0)
     return SquareWindow{Float64}(c, w)
 end
 
-function rectangle_square_window(c::AbstractVector, w::Union{Real, AbstractVector})
-    return w isa Real ? SquareWindow(c, w) : RectangleWindow(c, w)
+function rectangle_square_window(c, w)
+    return allequal(w) ? SquareWindow(c, w[1]) : RectangleWindow(c, w)
 end
 
 struct BallWindow{T<:Float64} <: AbstractSpatialWindow{T}

test/hard_core_spatial.jl

@@ -12,6 +12,6 @@ hc = PRS.HardCorePointProcess(b, r, win)
 
 seed = -1
 @time sample = PRS.generate_sample_prs(hc; rng=seed)
-# @time sample = PRS.generate_sample_grid_prs(hc; rng=-1)
+# @time sample = PRS.generate_sample_grid_prs(hc; rng=seed)
 
 p = PRS.plot(hc, sample; title="")

test/strauss_grid_prs.jl

@@ -5,7 +5,7 @@ const PRS = PartialRejectionSampling
 r = 0.05  # interaction range = 2*radius
 β, γ = β₀ / (π * (r/2)^2), 0.1  # γ = 0 ≡ Hard core
 
-c, w = [0.0, 0.0], 1.0
+c, w = [0.0, 0.0], 2.0
 win = PRS.SquareWindow(c, w)
 
 strauss = PRS.StraussPointProcess(β, γ, r, win)