diff --git a/Project.toml b/Project.toml index a90aea82..ca802606 100644 --- a/Project.toml +++ b/Project.toml @@ -22,4 +22,4 @@ LinearAlgebra = "1" Polyester = "0.7.5" Reexport = "1" StaticArrays = "1" -julia = "1.10" +julia = "1.10" \ No newline at end of file diff --git a/src/PointNeighbors.jl b/src/PointNeighbors.jl index f6458409..6e7ad93e 100644 --- a/src/PointNeighbors.jl +++ b/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 diff --git a/src/cell_lists/cell_lists.jl b/src/cell_lists/cell_lists.jl index fef8c299..cde404ed 100644 --- a/src/cell_lists/cell_lists.jl +++ b/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") diff --git a/src/cell_lists/spatial_hashing.jl b/src/cell_lists/spatial_hashing.jl new file mode 100644 index 00000000..4822cb2c --- /dev/null +++ b/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 diff --git a/src/nhs_grid.jl b/src/nhs_grid.jl index 08139ba9..ba613ecd 100644 --- a/src/nhs_grid.jl +++ b/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) diff --git a/test/cell_lists/spatial_hashing.jl b/test/cell_lists/spatial_hashing.jl new file mode 100644 index 00000000..2c83f611 --- /dev/null +++ b/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 diff --git a/test/neighborhood_search.jl b/test/neighborhood_search.jl index ac1da6a6..a1ef9b3b 100644 --- a/test/neighborhood_search.jl +++ b/test/neighborhood_search.jl @@ -58,7 +58,11 @@ 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 = [ @@ -66,7 +70,8 @@ "`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)