Spatial Hashing

Project.toml

@@ -22,4 +22,4 @@ LinearAlgebra = "1"
 Polyester = "0.7.5"
 Reexport = "1"
 StaticArrays = "1"
-julia = "1.10"
+julia = "1.10"

src/PointNeighbors.jl

@@ -22,7 +22,7 @@ include("gpu.jl")
 
 export foreach_point_neighbor, foreach_neighbor
 export TrivialNeighborhoodSearch, GridNeighborhoodSearch, PrecomputedNeighborhoodSearch
-export DictionaryCellList, FullGridCellList
+export DictionaryCellList, FullGridCellList, SpatialHashingCellList
 export ParallelUpdate, SemiParallelUpdate, SerialIncrementalUpdate, SerialUpdate,
        ParallelIncrementalUpdate
 export requires_update

src/cell_lists/cell_lists.jl

@@ -6,3 +6,4 @@ abstract type AbstractCellList end
 
 include("dictionary.jl")
 include("full_grid.jl")
+include("spatial_hashing.jl")

src/cell_lists/spatial_hashing.jl

@@ -0,0 +1,108 @@
+"""
+    SpatialHashingCellList{NDIMS}(; list_size)
+
+A basic spatial hashing implementation. Similar to [`DictionaryCellList`](@ref), the domain is discretized into cells, 
+and the particles in each cell are stored in a hash map. The hash is computed using the spatial location of each cell, 
+as described by Ihmsen et al. (2011)(@cite Ihmsen2011). By using a hash map that stores entries only for non-empty cells, 
+the domain is effectively infinite. The size of the hash map is recommended to be approximately twice the number of particles 
+to balance memory consumption against the likelihood of hash collisions.
+
+# Arguments
+- `NDIMS::Int`: Number of spatial dimensions (e.g., `2` or `3`).
+- `list_size::Int`: Size of the hash map (e.g., `2 * n_points`) .
+"""
+
+struct SpatialHashingCellList{NDIMS, CL, CI, CF} <: AbstractCellList
+    points     :: CL
+    coords     :: CI
+    collisions :: CF
+    list_size  :: Int
+end
+
+@inline index_type(::SpatialHashingCellList) = Int32
+
+@inline Base.ndims(::SpatialHashingCellList{NDIMS}) where {NDIMS} = NDIMS
+
+function supported_update_strategies(::SpatialHashingCellList)
+    return (SerialUpdate,)
+end
+
+function SpatialHashingCellList{NDIMS}(list_size) where {NDIMS}
+    points = [Int[] for _ in 1:list_size]
+    collisions = [false for _ in 1:list_size]
+    coords = [ntuple(_ -> typemin(Int), NDIMS) for _ in 1:list_size]
+    return SpatialHashingCellList{NDIMS, typeof(points), typeof(coords),
+                                  typeof(collisions)}(points, coords, collisions, list_size)
+end
+
+function Base.empty!(cell_list::SpatialHashingCellList)
+    (; list_size) = cell_list
+    NDIMS = ndims(cell_list)
+
+    Base.empty!.(cell_list.points)
+    cell_list.coords .= [ntuple(_ -> typemin(Int), NDIMS) for _ in 1:list_size]
+    cell_list.collisions .= false
+    return cell_list
+end
+
+# For each entry in the hash table, store the coordinates of the cell of the first point being inserted at this entry.
+# If a point with a different cell coordinate is being added, we have found a collision.
+function push_cell!(cell_list::SpatialHashingCellList, cell, point)
+    (; points, coords, collisions, list_size) = cell_list
+    NDIMS = ndims(cell_list)
+    hash_key = spatial_hash(cell, list_size)
+    push!(points[hash_key], point)
+
+    cell_coord = coords[hash_key]
+    if cell_coord == ntuple(_ -> typemin(Int), NDIMS)
+        # If this cell is not used yet, set cell coordinates
+        coords[hash_key] = cell
+    elseif cell_coord != cell
+        # If it is already used by a different cell, mark as collision
+        collisions[hash_key] = true
+    end
+end
+
+function deleteat_cell!(cell_list::SpatialHashingCellList, cell, i)
+    deleteat!(cell_list[cell], i)
+end
+
+@inline each_cell_index(cell_list::SpatialHashingCellList) = eachindex(cell_list.points)
+
+function copy_cell_list(cell_list::SpatialHashingCellList, search_radius,
+                        periodic_box)
+    (; list_size) = cell_list
+    NDIMS = ndims(cell_list)
+
+    return SpatialHashingCellList{NDIMS}(list_size)
+end
+
+@inline function Base.getindex(cell_list::SpatialHashingCellList, cell::Tuple)
+    return cell_list.points[spatial_hash(cell, length(cell_list.points))]
+end
+
+@inline function Base.getindex(cell_list::SpatialHashingCellList, i::Integer)
+    return cell_list.points[i]
+end
+
+@inline function is_correct_cell(cell_list::SpatialHashingCellList{<:Any, Nothing},
+                                 coords, cell_index::Array)
+    return coords == cell_index
+end
+
+# Hash functions according to Ihmsen et al. (2001)
+function spatial_hash(cell::NTuple{1, Real}, list_size)
+    return mod(cell[1] * 73856093, list_size) + 1
+end
+
+function spatial_hash(cell::NTuple{2, Real}, list_size)
+    i, j = cell
+
+    return mod(xor(i * 73856093, j * 19349663), list_size) + 1
+end
+
+function spatial_hash(cell::NTuple{3, Real}, list_size)
+    i, j, k = cell
+
+    return mod(xor(i * 73856093, j * 19349663, k * 83492791), list_size) + 1
+end

src/nhs_grid.jl

@@ -449,19 +449,58 @@ function update_grid!(neighborhood_search::Union{GridNeighborhoodSearch{<:Any,
     initialize_grid!(neighborhood_search, y; parallelization_backend, eachindex_y)
 end
 
+function check_collision(neighbor_cell_, neighbor_coords, cell_list, nhs)
+    # This is only relevant for the `SpatialHashingCellList`
+    return false
+end
+
+# Check if `neighbor_coords` belong to `neighbor_cell`, which might not be the case
+# with the `SpatialHashingCellList` if this cell has a collision.
+function check_collision(neighbor_cell_::CartesianIndex, neighbor_coords,
+                         cell_list::SpatialHashingCellList, nhs)
+    (; list_size, collisions, coords) = cell_list
+    neighbor_cell = periodic_cell_index(Tuple(neighbor_cell_), nhs)
+
+    return neighbor_cell != cell_coords(neighbor_coords, nhs)
+end
+
+function check_cell_collision(neighbor_cell_::CartesianIndex,
+                              cell_list, nhs)
+    # This is only relevant for the `SpatialHashingCellList`
+    return false
+end
+
+# Check if there is a collision in this cell, meaning there is at least one point
+# in this list that doesn't actually belong in this cell.
+function check_cell_collision(neighbor_cell_::CartesianIndex,
+                              cell_list::SpatialHashingCellList, nhs)
+    (; list_size, collisions, coords) = cell_list
+    neighbor_cell = periodic_cell_index(Tuple(neighbor_cell_), nhs)
+    hash = spatial_hash(neighbor_cell, list_size)
+
+    # `collisions[hash] == true` means points from multiple cells are in this list.
+    # `collisions[hash] == false` means points from only one cells are in this list.
+    # We could still have a collision though, if this one cell is not `neighbor_cell`,
+    # which is possible when `neighbor_cell` is empty.
+    return collisions[hash] || coords[hash] != neighbor_cell
+end
+
 # Specialized version of the function in `neighborhood_search.jl`, which is faster
 # than looping over `eachneighbor`.
 @inline function foreach_neighbor(f, neighbor_system_coords,
                                   neighborhood_search::GridNeighborhoodSearch,
                                   point, point_coords, search_radius)
-    (; periodic_box) = neighborhood_search
-
+    (; cell_list, periodic_box) = neighborhood_search
     cell = cell_coords(point_coords, neighborhood_search)
 
     for neighbor_cell_ in neighboring_cells(cell, neighborhood_search)
         neighbor_cell = Tuple(neighbor_cell_)
         neighbors = points_in_cell(neighbor_cell, neighborhood_search)
 
+        # Boolean to indicate if this cell has a collision (only with `SpatialHashingCellList`)
+        cell_collision = check_cell_collision(neighbor_cell_,
+                                              cell_list, neighborhood_search)
+
         for neighbor_ in eachindex(neighbors)
             neighbor = @inbounds neighbors[neighbor_]
 
@@ -480,6 +519,14 @@ end
             if distance2 <= search_radius^2
                 distance = sqrt(distance2)
 
+                # If this cell has a collision, check if this point belongs to this cell
+                # (only with `SpatialHashingCellList`).
+                if cell_collision &&
+                   check_collision(neighbor_cell_, neighbor_coords, cell_list,
+                                   neighborhood_search)
+                    continue
+                end
+
                 # Inline to avoid loss of performance
                 # compared to not using `foreach_point_neighbor`.
                 @inline f(point, neighbor, pos_diff, distance)

test/cell_lists/spatial_hashing.jl

@@ -0,0 +1,69 @@
+@testset verbose=true "SpatialHashingCellList" begin
+    @testset "Collision Handling With Empty Cells" begin
+        # The point is in cell (-1, 0) which has a hash collision with cell (-2, -1)
+        coordinates = [-0.05; 0.05;;]
+        NDIMS = size(coordinates, 1)
+        n_points = size(coordinates, 2)
+        search_radius = 0.1 + 10 * eps()
+        point_index = 1
+
+        nhs = GridNeighborhoodSearch{2}(; search_radius, n_points,
+                                        cell_list = SpatialHashingCellList{NDIMS}(n_points))
+        initialize_grid!(nhs, coordinates)
+
+        @testset "Test For Collision" begin
+            cell1 = (-1, 0)
+            cell2 = (-2, -1)
+            cell1_hash = PointNeighbors.spatial_hash(cell1, n_points)
+            cell2_hash = PointNeighbors.spatial_hash(cell2, n_points)
+            points1 = nhs.cell_list[cell1]
+            points2 = nhs.cell_list[cell2]
+
+            @test points1 == points2 == [1]
+            @test cell1_hash == cell2_hash
+        end
+
+        neighbors = Int[]
+        foreach_neighbor(coordinates, coordinates, nhs,
+                         point_index) do point, neighbor, pos_diff, distance
+            push!(neighbors, neighbor)
+        end
+
+        @test neighbors == [1]
+    end
+
+    @testset "Collision Handling With Non-Empty Cells" begin
+        # Cell (-1, 0) with point 1 has a hash collision with cell (-2, -1) with point 2
+        coordinates = [[-0.05 -0.15]; [0.05 -0.05]]
+        NDIMS = size(coordinates, 1)
+        n_points = size(coordinates, 2)
+        search_radius = 0.1 + 10 * eps()
+        point_index = 1
+
+        nhs = GridNeighborhoodSearch{2}(; search_radius, n_points,
+                                        cell_list = SpatialHashingCellList{NDIMS}(n_points))
+        initialize_grid!(nhs, coordinates)
+
+        @testset "Test For Collision" begin
+            cell1 = (-1, 0)
+            cell2 = (-2, -1)
+            cell1_hash = PointNeighbors.spatial_hash(cell1, n_points)
+            cell2_hash = PointNeighbors.spatial_hash(cell2, n_points)
+            points1 = nhs.cell_list[cell1]
+            points2 = nhs.cell_list[cell2]
+
+            @test points1 == points2 == [1, 2]
+            @test cell1_hash == cell2_hash
+        end
+
+        neighbors = [Int[] for _ in axes(coordinates, 2)]
+        foreach_point_neighbor(coordinates, coordinates, nhs,
+                               points = axes(coordinates, 2)) do point, neighbor, pos_diff,
+                                                                 distance
+            push!(neighbors[point], neighbor)
+        end
+
+        @test neighbors[1] == [1]
+        @test neighbors[2] == [2]
+    end
+end

test/neighborhood_search.jl

@@ -58,15 +58,20 @@
                                                                            search_radius,
                                                                            backend = Vector{Vector{Int32}})),
                 PrecomputedNeighborhoodSearch{NDIMS}(; search_radius, n_points,
-                                                     periodic_box = periodic_boxes[i])
+                                                     periodic_box = periodic_boxes[i]),
+                GridNeighborhoodSearch{NDIMS}(; search_radius, n_points,
+                                              periodic_box = periodic_boxes[i],
+                                              cell_list = SpatialHashingCellList{NDIMS}(2 *
+                                                                                        n_points))
             ]
 
             names = [
                 "`TrivialNeighborhoodSearch`",
                 "`GridNeighborhoodSearch`",
                 "`GridNeighborhoodSearch` with `FullGridCellList` with `DynamicVectorOfVectors`",
                 "`GridNeighborhoodSearch` with `FullGridCellList` with `Vector{Vector}`",
-                "`PrecomputedNeighborhoodSearch`"
+                "`PrecomputedNeighborhoodSearch`",
+                "`GridNeighborhoodSearch` with `SpatialHashingCellList`"
             ]
 
             # Also test copied templates
@@ -80,7 +85,10 @@
                                               cell_list = FullGridCellList(min_corner = periodic_boxes[i].min_corner,
                                                                            max_corner = periodic_boxes[i].max_corner,
                                                                            backend = Vector{Vector{Int32}})),
-                PrecomputedNeighborhoodSearch{NDIMS}(periodic_box = periodic_boxes[i])
+                PrecomputedNeighborhoodSearch{NDIMS}(periodic_box = periodic_boxes[i]),
+                GridNeighborhoodSearch{NDIMS}(periodic_box = periodic_boxes[i],
+                                              cell_list = SpatialHashingCellList{NDIMS}(2 *
+                                                                                        n_points))
             ]
             copied_nhs = copy_neighborhood_search.(template_nhs, search_radius, n_points)
             append!(neighborhood_searches, copied_nhs)
@@ -184,7 +192,10 @@
                                                                            max_corner,
                                                                            search_radius,
                                                                            backend = Vector{Vector{Int}})),
-                PrecomputedNeighborhoodSearch{NDIMS}(; search_radius, n_points)
+                PrecomputedNeighborhoodSearch{NDIMS}(; search_radius, n_points),
+                GridNeighborhoodSearch{NDIMS}(; search_radius, n_points,
+                                              cell_list = SpatialHashingCellList{NDIMS}(2 *
+                                                                                        n_points))
             ]
 
             names = [
@@ -195,7 +206,8 @@
                 "`GridNeighborhoodSearch` with `FullGridCellList` with `DynamicVectorOfVectors` and `ParallelIncrementalUpdate`",
                 "`GridNeighborhoodSearch` with `FullGridCellList` with `DynamicVectorOfVectors` and `SemiParallelUpdate`",
                 "`GridNeighborhoodSearch` with `FullGridCellList` with `Vector{Vector}`",
-                "`PrecomputedNeighborhoodSearch`"
+                "`PrecomputedNeighborhoodSearch`",
+                "`GridNeighborhoodSearch` with `SpatialHashingCellList`"
             ]
 
             # Also test copied templates
@@ -214,7 +226,9 @@
                 GridNeighborhoodSearch{NDIMS}(cell_list = FullGridCellList(; min_corner,
                                                                            max_corner,
                                                                            backend = Vector{Vector{Int32}})),
-                PrecomputedNeighborhoodSearch{NDIMS}()
+                PrecomputedNeighborhoodSearch{NDIMS}(),
+                GridNeighborhoodSearch{NDIMS}(cell_list = SpatialHashingCellList{NDIMS}(2 *
+                                                                                        n_points))
             ]
             copied_nhs = copy_neighborhood_search.(template_nhs, search_radius, n_points)
             append!(neighborhood_searches, copied_nhs)