From e0c971907c9158df1592ada6aeb2c32ca86f7470 Mon Sep 17 00:00:00 2001 From: RubberLanding Date: Wed, 19 Feb 2025 19:14:33 +0800 Subject: [PATCH 01/27] Added test file for Spatial Hashing. Updated the module file to expose SpatialHashingCellList to the users. --- src/PointNeighbors.jl | 3 +- src/nhs_spatial_hash.jl | 91 ++++++++++++++++++++++++++++++++++++++++ test/nhs_spatial_hash.jl | 69 ++++++++++++++++++++++++++++++ 3 files changed, 162 insertions(+), 1 deletion(-) create mode 100644 src/nhs_spatial_hash.jl create mode 100644 test/nhs_spatial_hash.jl diff --git a/src/PointNeighbors.jl b/src/PointNeighbors.jl index 22ef3893..2154b571 100644 --- a/src/PointNeighbors.jl +++ b/src/PointNeighbors.jl @@ -19,9 +19,10 @@ include("cell_lists/cell_lists.jl") include("nhs_grid.jl") include("nhs_precomputed.jl") include("gpu.jl") +include("nhs_spatial_hash.jl") export foreach_point_neighbor, foreach_neighbor -export TrivialNeighborhoodSearch, GridNeighborhoodSearch, PrecomputedNeighborhoodSearch +export TrivialNeighborhoodSearch, GridNeighborhoodSearch, PrecomputedNeighborhoodSearch, SpatialHashingCellList export DictionaryCellList, FullGridCellList export ParallelUpdate, SemiParallelUpdate, SerialUpdate export initialize!, update!, initialize_grid!, update_grid! diff --git a/src/nhs_spatial_hash.jl b/src/nhs_spatial_hash.jl new file mode 100644 index 00000000..0422d32c --- /dev/null +++ b/src/nhs_spatial_hash.jl @@ -0,0 +1,91 @@ +struct SpatialHashingCellList{CL, CI, CF} + cell_points :: CL + cell_coords :: CI + cell_collision :: CF + size :: Int +end + +function SpatialHashingCellList{NDIMS}(size; handle_collisions = true) where {NDIMS} + cell_points = [Int[] for _ in 1:size] + cell_collision = [false for _ in 1:size] + cell_coords = [Int[] for _ in 1:size] + return SpatialHashingCellList(cell_points, cell_coords, cell_collision, size) +end + +function Base.empty!(cell_list::SpatialHashingCellList) + Base.empty!.(cell_list.cell_points) + Base.empty!.(cell_list.cell_coords) + cell_list.cell_collision .= false + + return cell_list +end + +function push_cell!(cell_list::SpatialHashingCellList, cell, point) + (; cell_points, cell_coords, cell_collision, size) = cell_list + key = spatial_hash(cell, size) + push!(cell_points[key], point) + + if !(cell in cell_coords[key]) + push!(cell_coords[key], cell) + end + + if length(cell_coords[key]) > 1 + cell_collision[key] = true + end +end + +function get_points(nhs, cell, coords_fun) + (; cell_list) = nhs + (; cell_points, cell_coords, cell_collision, size) = cell_list + key = spatial_hash(cell, size) + if cell_collision[key] + points_in_cell = [] + for point in cell_list.cell_points[key] + if cell_coords(coords_fun(point), nhs) == cell + push!(points_in_cell, point) + end + end + return points_in_cell + + else + return cell_list.cell_points[key] + 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.cell_points) + +@inline function Base.getindex(cell_list::SpatialHashingCellList, i::Integer) + return cell_list.cell_points[i] +end + +@inline function is_correct_cell(cell_list::SpatialHashingCellList{<:Any, Nothing}, + cell_coords, cell_index::Array) # What is the correct type for cell_index? It is a list of integers +return cell_coords in cell_list.cell_coords[cell_index] +end + +function spatial_hash(cell::NTuple{1, Real}, size) + return mod(cell[1] * 73856093, size) + 1 +end + +function spatial_hash(cell::NTuple{2, Real}, size::UInt64) + i, j = cell + + return mod(xor(i * 73856093, j * 19349663), size) + 1 +end + +function spatial_hash(cell::NTuple{3, Real}, size) + i, j, k = cell + + return mod(xor(i * 73856093, j * 19349663, k * 83492791), size) + 1 +end + +function spatial_hash(cell::NTuple{3, Real}, index, size) + i, j, k = cell + + return mod(xor(i * 73856093, j * 19349663, k * 83492791, index * 7238423947), size) + 1 +end diff --git a/test/nhs_spatial_hash.jl b/test/nhs_spatial_hash.jl new file mode 100644 index 00000000..78fbcaac --- /dev/null +++ b/test/nhs_spatial_hash.jl @@ -0,0 +1,69 @@ +@testset verbose=true "SpatialHashNeighborhoodSearch" begin + @testset "Rectangular Point Cloud 2D" begin + #### Setup + # Rectangle of equidistantly spaced points + # from (x, y) = (-0.25, -0.25) to (x, y) = (0.35, 0.35). + range = -0.25:0.1:0.35 + coordinates1 = hcat(collect.(Iterators.product(range, range))...) + n_points = size(coordinates1, 2) + + point_position1 = [0.05, 0.05] + search_radius = 0.1 + + @infiltrate + + # Create neighborhood search + nhs1 = GridNeighborhoodSearch{2}(; search_radius, n_points, + cell_list = SpatialHashingCellList{NDIMS}(n_points)) + + initialize_grid!(nhs1, coordinates1) + + # Get each neighbor for `point_position1` + neighbors1 = sort(collect(PointNeighbors.eachneighbor(point_position1, nhs1))) + + # Move points + coordinates2 = coordinates1 .+ [1.4, -3.5] + + # Update neighborhood search + update_grid!(nhs1, coordinates2) + + # Get each neighbor for updated NHS + neighbors2 = sort(collect(PointNeighbors.eachneighbor(point_position1, nhs1))) + + # Change position + point_position2 = point_position1 .+ [1.4, -3.5] + + # Get each neighbor for `point_position2` + neighbors3 = sort(collect(PointNeighbors.eachneighbor(point_position2, nhs1))) + + # Double search radius + nhs2 = GridNeighborhoodSearch{2}(search_radius = 2 * search_radius, + n_points = size(coordinates1, 2)) + initialize!(nhs2, coordinates1, coordinates1) + + # Get each neighbor in double search radius + neighbors4 = sort(collect(PointNeighbors.eachneighbor(point_position1, nhs2))) + + # Move points + coordinates2 = coordinates1 .+ [0.4, -0.4] + + # Update neighborhood search + update!(nhs2, coordinates2, coordinates2) + + # Get each neighbor in double search radius + neighbors5 = sort(collect(PointNeighbors.eachneighbor(point_position1, nhs2))) + + #### Verification against lists of potential neighbors built by hand + @test neighbors1 == [17, 18, 19, 24, 25, 26, 31, 32, 33] + + @test neighbors2 == Int[] + + @test neighbors3 == [17, 18, 19, 24, 25, 26, 31, 32, 33] + + @test neighbors4 == [ + 9, 10, 11, 12, 13, 14, 16, 17, 18, 19, 20, 21, 23, 24, 25, 26, 27, 28, 30, 31, + 32, 33, 34, 35, 37, 38, 39, 40, 41, 42, 44, 45, 46, 47, 48, 49] + + @test neighbors5 == [36, 37, 38, 43, 44, 45] + end +end; From 22586792868c82e032d4d6c5ef498c03ec270d20 Mon Sep 17 00:00:00 2001 From: RubberLanding Date: Thu, 27 Feb 2025 12:54:14 +0800 Subject: [PATCH 02/27] Minor changes. --- src/PointNeighbors.jl | 5 +- src/cell_lists/cell_lists.jl | 1 + src/cell_lists/spatial_hashing.jl | 102 ++++++++++++++++++++++++++++++ test/nhs_spatial_hash.jl | 4 +- 4 files changed, 106 insertions(+), 6 deletions(-) create mode 100644 src/cell_lists/spatial_hashing.jl diff --git a/src/PointNeighbors.jl b/src/PointNeighbors.jl index 2154b571..8596e513 100644 --- a/src/PointNeighbors.jl +++ b/src/PointNeighbors.jl @@ -19,11 +19,10 @@ include("cell_lists/cell_lists.jl") include("nhs_grid.jl") include("nhs_precomputed.jl") include("gpu.jl") -include("nhs_spatial_hash.jl") export foreach_point_neighbor, foreach_neighbor -export TrivialNeighborhoodSearch, GridNeighborhoodSearch, PrecomputedNeighborhoodSearch, SpatialHashingCellList -export DictionaryCellList, FullGridCellList +export TrivialNeighborhoodSearch, GridNeighborhoodSearch, PrecomputedNeighborhoodSearch +export DictionaryCellList, FullGridCellList, SpatialHashingCellList export ParallelUpdate, SemiParallelUpdate, SerialUpdate export initialize!, update!, initialize_grid!, update_grid! export PolyesterBackend, ThreadsDynamicBackend, ThreadsStaticBackend 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..50be7b08 --- /dev/null +++ b/src/cell_lists/spatial_hashing.jl @@ -0,0 +1,102 @@ +struct SpatialHashingCellList{CL, CI, CF} + cell_points :: CL + cell_coords :: CI + cell_collision :: CF + size :: Int +end + +@inline index_type(::SpatialHashingCellList) = Int64 + +# Change the order back to (SemiParallelUpdate, SerialUpdate) when SemiParallelUpdate is implemented +function supported_update_strategies(::SpatialHashingCellList) + return (SerialUpdate, SemiParallelUpdate) +end + +function SpatialHashingCellList{NDIMS}(size; handle_collisions = true) where {NDIMS} + cell_points = [Int[] for _ in 1:size] + cell_collision = [false for _ in 1:size] + cell_coords = [NTuple{NDIMS, Int}[] for _ in 1:size] + return SpatialHashingCellList(cell_points, cell_coords, cell_collision, size) +end + +function Base.empty!(cell_list::SpatialHashingCellList) + Base.empty!.(cell_list.cell_points) + Base.empty!.(cell_list.cell_coords) + cell_list.cell_collision .= false + + return cell_list +end + +function push_cell!(cell_list::SpatialHashingCellList, cell, point) + (; cell_points, cell_coords, cell_collision, size) = cell_list + key = spatial_hash(cell, size) + push!(cell_points[key], point) + + if !(cell in cell_coords[key]) + push!(cell_coords[key], cell) + end + + if length(cell_coords[key]) > 1 + cell_collision[key] = true + end +end + +function get_points(nhs, cell, coords_fun) + (; cell_list) = nhs + (; cell_points, cell_coords, cell_collision, size) = cell_list + key = spatial_hash(cell, size) + if cell_collision[key] + points_in_cell = [] + for point in cell_list.cell_points[key] + if cell_coords(coords_fun(point), nhs) == cell + push!(points_in_cell, point) + end + end + return points_in_cell + + else + return cell_list.cell_points[key] + 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.cell_points) + +@inline function Base.getindex(cell_list::SpatialHashingCellList, i) + if isa(i, Int) + return cell_list.cell_points[i] + elseif isa(i, Tuple) + return cell_list.cell_points[i...] + end +end + +@inline function is_correct_cell(cell_list::SpatialHashingCellList{<:Any, Nothing}, + cell_coords, cell_index::Array) # What is the correct type for cell_index? It is a list of integers +return cell_coords in cell_list.cell_coords[cell_index] +end + +function spatial_hash(cell::NTuple{1, Real}, size) + return mod(cell[1] * 73856093, size) + 1 +end + +function spatial_hash(cell::NTuple{2, Real}, size) + i, j = cell + + return mod(xor(i * 73856093, j * 19349663), size) + 1 +end + +function spatial_hash(cell::NTuple{3, Real}, size) + i, j, k = cell + + return mod(xor(i * 73856093, j * 19349663, k * 83492791), size) + 1 +end + +function spatial_hash(cell::NTuple{3, Real}, index, size) + i, j, k = cell + + return mod(xor(i * 73856093, j * 19349663, k * 83492791, index * 7238423947), size) + 1 +end diff --git a/test/nhs_spatial_hash.jl b/test/nhs_spatial_hash.jl index 78fbcaac..d810f51d 100644 --- a/test/nhs_spatial_hash.jl +++ b/test/nhs_spatial_hash.jl @@ -10,11 +10,9 @@ point_position1 = [0.05, 0.05] search_radius = 0.1 - @infiltrate - # Create neighborhood search nhs1 = GridNeighborhoodSearch{2}(; search_radius, n_points, - cell_list = SpatialHashingCellList{NDIMS}(n_points)) + cell_list = SpatialHashingCellList{2}(n_points)) initialize_grid!(nhs1, coordinates1) From 4b57a71dab9dc6ff8b1b2e4d95469f2be02ffa61 Mon Sep 17 00:00:00 2001 From: RubberLanding Date: Thu, 27 Feb 2025 17:03:38 +0800 Subject: [PATCH 03/27] Spatial Hashing works with initialize!() and empty!(). --- src/cell_lists/spatial_hashing.jl | 66 +++++++++++++++---------------- test/nhs_spatial_hash.jl | 33 +++++++++------- 2 files changed, 51 insertions(+), 48 deletions(-) diff --git a/src/cell_lists/spatial_hashing.jl b/src/cell_lists/spatial_hashing.jl index 50be7b08..bd185592 100644 --- a/src/cell_lists/spatial_hashing.jl +++ b/src/cell_lists/spatial_hashing.jl @@ -1,8 +1,8 @@ struct SpatialHashingCellList{CL, CI, CF} - cell_points :: CL - cell_coords :: CI - cell_collision :: CF - size :: Int + size::Int + cell_points::CL + cell_coords::CI + cell_collision::CF end @inline index_type(::SpatialHashingCellList) = Int64 @@ -23,7 +23,6 @@ function Base.empty!(cell_list::SpatialHashingCellList) Base.empty!.(cell_list.cell_points) Base.empty!.(cell_list.cell_coords) cell_list.cell_collision .= false - return cell_list end @@ -31,52 +30,33 @@ function push_cell!(cell_list::SpatialHashingCellList, cell, point) (; cell_points, cell_coords, cell_collision, size) = cell_list key = spatial_hash(cell, size) push!(cell_points[key], point) - if !(cell in cell_coords[key]) push!(cell_coords[key], cell) end - if length(cell_coords[key]) > 1 cell_collision[key] = true end end -function get_points(nhs, cell, coords_fun) - (; cell_list) = nhs - (; cell_points, cell_coords, cell_collision, size) = cell_list - key = spatial_hash(cell, size) - if cell_collision[key] - points_in_cell = [] - for point in cell_list.cell_points[key] - if cell_coords(coords_fun(point), nhs) == cell - push!(points_in_cell, point) - end - end - return points_in_cell - - else - return cell_list.cell_points[key] - 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.cell_points) -@inline function Base.getindex(cell_list::SpatialHashingCellList, i) - if isa(i, Int) - return cell_list.cell_points[i] - elseif isa(i, Tuple) - return cell_list.cell_points[i...] - end +@inline function Base.getindex(cell_list::SpatialHashingCellList, cell::Tuple) + (; cell_points) = cell_list + + return cell_points[spatial_hash(cell, length(cell_points))] +end + +@inline function Base.getindex(cell_list::SpatialHashingCellList, i::Integer) + return cell_list.cell_points[i] end @inline function is_correct_cell(cell_list::SpatialHashingCellList{<:Any, Nothing}, - cell_coords, cell_index::Array) # What is the correct type for cell_index? It is a list of integers -return cell_coords in cell_list.cell_coords[cell_index] + cell_coords, cell_index::Array) # What is the correct type for cell_index? It is a list of integers + return cell_coords in cell_list.cell_coords[cell_index] end function spatial_hash(cell::NTuple{1, Real}, size) @@ -100,3 +80,21 @@ function spatial_hash(cell::NTuple{3, Real}, index, size) return mod(xor(i * 73856093, j * 19349663, k * 83492791, index * 7238423947), size) + 1 end + +# function get_points(nhs, cell, coords_fun) +# (; cell_list) = nhs +# (; cell_points, cell_coords, cell_collision, size) = cell_list +# key = spatial_hash(cell, size) +# if cell_collision[key] +# points_in_cell = [] +# for point in cell_list.cell_points[key] +# if cell_coords(coords_fun(point), nhs) == cell +# push!(points_in_cell, point) +# end +# end +# return points_in_cell + +# else +# return cell_list.cell_points[key] +# end +# end diff --git a/test/nhs_spatial_hash.jl b/test/nhs_spatial_hash.jl index d810f51d..608976c0 100644 --- a/test/nhs_spatial_hash.jl +++ b/test/nhs_spatial_hash.jl @@ -12,18 +12,17 @@ # Create neighborhood search nhs1 = GridNeighborhoodSearch{2}(; search_radius, n_points, - cell_list = SpatialHashingCellList{2}(n_points)) + cell_list = SpatialHashingCellList{2}(2 * n_points)) initialize_grid!(nhs1, coordinates1) # Get each neighbor for `point_position1` neighbors1 = sort(collect(PointNeighbors.eachneighbor(point_position1, nhs1))) - # Move points + # Move points and update NHS coordinates2 = coordinates1 .+ [1.4, -3.5] - - # Update neighborhood search - update_grid!(nhs1, coordinates2) + empty!(nhs1.cell_list) + initialize_grid!(nhs1, coordinates2) # Get each neighbor for updated NHS neighbors2 = sort(collect(PointNeighbors.eachneighbor(point_position1, nhs1))) @@ -36,7 +35,9 @@ # Double search radius nhs2 = GridNeighborhoodSearch{2}(search_radius = 2 * search_radius, - n_points = size(coordinates1, 2)) + n_points = size(coordinates1, 2), + cell_list = SpatialHashingCellList{2}(2 * n_points)) + initialize!(nhs2, coordinates1, coordinates1) # Get each neighbor in double search radius @@ -44,24 +45,28 @@ # Move points coordinates2 = coordinates1 .+ [0.4, -0.4] - - # Update neighborhood search - update!(nhs2, coordinates2, coordinates2) + empty!(nhs2.cell_list) + initialize_grid!(nhs2, coordinates2) # Get each neighbor in double search radius neighbors5 = sort(collect(PointNeighbors.eachneighbor(point_position1, nhs2))) #### Verification against lists of potential neighbors built by hand - @test neighbors1 == [17, 18, 19, 24, 25, 26, 31, 32, 33] + correct_neighbors1 = [17, 18, 19, 24, 25, 26, 31, 32, 33] + @test all(x -> x in neighbors1, correct_neighbors1) - @test neighbors2 == Int[] + correct_neighbors2 = Int[] + @test all(x -> x in neighbors2, correct_neighbors2) - @test neighbors3 == [17, 18, 19, 24, 25, 26, 31, 32, 33] + correct_neighbors3 = [17, 18, 19, 24, 25, 26, 31, 32, 33] + @test all(x -> x in neighbors3, correct_neighbors3) - @test neighbors4 == [ + correct_neighbors4 = [ 9, 10, 11, 12, 13, 14, 16, 17, 18, 19, 20, 21, 23, 24, 25, 26, 27, 28, 30, 31, 32, 33, 34, 35, 37, 38, 39, 40, 41, 42, 44, 45, 46, 47, 48, 49] + @test all(x -> x in neighbors4, correct_neighbors4) - @test neighbors5 == [36, 37, 38, 43, 44, 45] + correct_neighbors5 = [36, 37, 38, 43, 44, 45] + @test all(x -> x in neighbors5, correct_neighbors5) end end; From ac4d23fb92208268dfcfdbd18a60bb07081613b7 Mon Sep 17 00:00:00 2001 From: RubberLanding Date: Thu, 27 Feb 2025 17:04:18 +0800 Subject: [PATCH 04/27] Minor changes. --- test/nhs_spatial_hash.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/test/nhs_spatial_hash.jl b/test/nhs_spatial_hash.jl index 608976c0..eb948440 100644 --- a/test/nhs_spatial_hash.jl +++ b/test/nhs_spatial_hash.jl @@ -13,7 +13,6 @@ # Create neighborhood search nhs1 = GridNeighborhoodSearch{2}(; search_radius, n_points, cell_list = SpatialHashingCellList{2}(2 * n_points)) - initialize_grid!(nhs1, coordinates1) # Get each neighbor for `point_position1` @@ -37,7 +36,6 @@ nhs2 = GridNeighborhoodSearch{2}(search_radius = 2 * search_radius, n_points = size(coordinates1, 2), cell_list = SpatialHashingCellList{2}(2 * n_points)) - initialize!(nhs2, coordinates1, coordinates1) # Get each neighbor in double search radius From c23840787b701025ecf9ee5d6fc6d457fa219c04 Mon Sep 17 00:00:00 2001 From: RubberLanding Date: Sun, 9 Mar 2025 12:21:26 +0800 Subject: [PATCH 05/27] Working version for collision handling. --- src/PointNeighbors.jl | 2 + src/cell_lists/spatial_hashing.jl | 97 +++++++++++++++++------------- src/nhs_grid.jl | 38 +++++++++++- test/cell_lists/spatial_hashing.jl | 89 +++++++++++++++++++++++++++ 4 files changed, 181 insertions(+), 45 deletions(-) create mode 100644 test/cell_lists/spatial_hashing.jl diff --git a/src/PointNeighbors.jl b/src/PointNeighbors.jl index 8596e513..bd6ba8f8 100644 --- a/src/PointNeighbors.jl +++ b/src/PointNeighbors.jl @@ -11,6 +11,8 @@ using LinearAlgebra: dot using Polyester: Polyester @reexport using StaticArrays: SVector +using Infiltrator: @infiltrate + include("util.jl") include("vector_of_vectors.jl") include("neighborhood_search.jl") diff --git a/src/cell_lists/spatial_hashing.jl b/src/cell_lists/spatial_hashing.jl index bd185592..eeeb1844 100644 --- a/src/cell_lists/spatial_hashing.jl +++ b/src/cell_lists/spatial_hashing.jl @@ -1,7 +1,7 @@ struct SpatialHashingCellList{CL, CI, CF} - size::Int + list_size::Int cell_points::CL - cell_coords::CI + cell_coords::CI # Think about renaming this field, there is a function with the same name cell_collision::CF end @@ -12,32 +12,52 @@ function supported_update_strategies(::SpatialHashingCellList) return (SerialUpdate, SemiParallelUpdate) end -function SpatialHashingCellList{NDIMS}(size; handle_collisions = true) where {NDIMS} - cell_points = [Int[] for _ in 1:size] - cell_collision = [false for _ in 1:size] - cell_coords = [NTuple{NDIMS, Int}[] for _ in 1:size] - return SpatialHashingCellList(cell_points, cell_coords, cell_collision, size) +function SpatialHashingCellList{NDIMS}(list_size) where {NDIMS} + cell_points = [Int[] for _ in 1:list_size] + cell_collision = [false for _ in 1:list_size] + + # cell_choords is used to check if there is a collision and contains + # the coordinates of the first cell added to the hash list + cell_coords = [missing for _ in 1:list_size] + return SpatialHashingCellList(list_size, cell_points, cell_coords, cell_collision) end function Base.empty!(cell_list::SpatialHashingCellList) Base.empty!.(cell_list.cell_points) - Base.empty!.(cell_list.cell_coords) + cell_list.cell_coords .= missing cell_list.cell_collision .= false return cell_list end +""" +cell::NTupl +e{NDIMS, Real}, tuple of cell coordinates +point::Int, index of the point to be added +""" + +@inline function cell_coords(coords, periodic_box::Nothing, cell_list::SpatialHashingCellList, + cell_size) +(; min_corner) = cell_list + +return Tuple(floor_to_int.((coords .- min_corner) ./ cell_size)) .+ 1 +end + + function push_cell!(cell_list::SpatialHashingCellList, cell, point) - (; cell_points, cell_coords, cell_collision, size) = cell_list - key = spatial_hash(cell, size) + (; cell_points, cell_coords, cell_collision, list_size) = cell_list + key = spatial_hash(cell, list_size) push!(cell_points[key], point) - if !(cell in cell_coords[key]) - push!(cell_coords[key], cell) - end - if length(cell_coords[key]) > 1 + + # Check if the a cell has been added at this hash + if ismissing(cell_coords) + cell_coords = cell + # Detect collision + elseif cell != cell_coords cell_collision[key] = true end end +# Implement reset of collision flag, if after the deletion there still is no collision? function deleteat_cell!(cell_list::SpatialHashingCellList, cell, i) deleteat!(cell_list[cell], i) end @@ -54,47 +74,38 @@ end return cell_list.cell_points[i] end +# Naming of cell_coords and cell_index confusing @inline function is_correct_cell(cell_list::SpatialHashingCellList{<:Any, Nothing}, - cell_coords, cell_index::Array) # What is the correct type for cell_index? It is a list of integers - return cell_coords in cell_list.cell_coords[cell_index] + cell_coords, cell_index::Array) + return cell_coords == cell_index end -function spatial_hash(cell::NTuple{1, Real}, size) - return mod(cell[1] * 73856093, size) + 1 +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}, size) +function spatial_hash(cell::NTuple{2, Real}, list_size) i, j = cell - return mod(xor(i * 73856093, j * 19349663), size) + 1 + return mod(xor(i * 73856093, j * 19349663), list_size) + 1 end -function spatial_hash(cell::NTuple{3, Real}, size) +function spatial_hash(cell::NTuple{3, Real}, list_size) i, j, k = cell - return mod(xor(i * 73856093, j * 19349663, k * 83492791), size) + 1 + return mod(xor(i * 73856093, j * 19349663, k * 83492791), list_size) + 1 end -function spatial_hash(cell::NTuple{3, Real}, index, size) +function spatial_hash(cell::NTuple{3, Real}, index, list_size) i, j, k = cell - return mod(xor(i * 73856093, j * 19349663, k * 83492791, index * 7238423947), size) + 1 -end - -# function get_points(nhs, cell, coords_fun) -# (; cell_list) = nhs -# (; cell_points, cell_coords, cell_collision, size) = cell_list -# key = spatial_hash(cell, size) -# if cell_collision[key] -# points_in_cell = [] -# for point in cell_list.cell_points[key] -# if cell_coords(coords_fun(point), nhs) == cell -# push!(points_in_cell, point) -# end -# end -# return points_in_cell - -# else -# return cell_list.cell_points[key] -# end -# end + return mod(xor(i * 73856093, j * 19349663, k * 83492791, index * 7238423947), list_size) + 1 +end + +function spatial_hash(cell::CartesianIndex{2}, list_size) + return spatial_hash(Tuple(cell), list_size) +end + +function spatial_hash(cell::CartesianIndex{3}, list_size) + return spatial_hash(Tuple(cell), list_size) +end diff --git a/src/nhs_grid.jl b/src/nhs_grid.jl index 47011b68..5fc72958 100644 --- a/src/nhs_grid.jl +++ b/src/nhs_grid.jl @@ -358,11 +358,28 @@ end point, point_coords, search_radius) end +function check_collision(neighbor_cell_, neighbor_coords, cell_list, nhs) + return false +end + +function check_collision(neighbor_cell_::CartesianIndex, neighbor_coords, + cell_list::SpatialHashingCellList, nhs) + neighbor_cell_ = Tuple(neighbor_cell_) + (; list_size, cell_collision) = cell_list + + index = spatial_hash(neighbor_cell_, list_size) + collision = cell_collision[index] + if collision + # Check if `neighbor_coords` are in the cell `neighbor_cell_` + cell_coords_ = cell_coords(neighbor_coords, nhs) + return neighbor_cell_ != cell_coords_ + end +end + @inline function __foreach_neighbor(f, system_coords, 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) @@ -377,6 +394,14 @@ end neighbor_coords = extract_svector(neighbor_system_coords, Val(ndims(neighborhood_search)), neighbor) + # Check if `neighbor_coords` are in the cell `neighbor_cell_`. + # For the `SpatialHashingCellList`, this might not be the case + # if we have a collision. + if check_collision(neighbor_cell_, neighbor_coords, cell_list, + neighborhood_search) + continue + end + pos_diff = point_coords - neighbor_coords distance2 = dot(pos_diff, pos_diff) @@ -384,6 +409,9 @@ end search_radius, periodic_box) + is_near = distance2 <= search_radius^2 + print("$neighbor, $distance2, $is_near \n") + if distance2 <= search_radius^2 distance = sqrt(distance2) @@ -436,6 +464,12 @@ end return cell_coords(coords, periodic_box, cell_list, cell_size) end +@inline function cell_coords(coords, neighborhood_search) + (; periodic_box, cell_list, cell_size) = neighborhood_search + + return cell_coords(coords, periodic_box, cell_list, cell_size) +end + @inline function cell_coords(coords, periodic_box::Nothing, cell_list, cell_size) return Tuple(floor_to_int.(coords ./ cell_size)) end diff --git a/test/cell_lists/spatial_hashing.jl b/test/cell_lists/spatial_hashing.jl new file mode 100644 index 00000000..86aef7b2 --- /dev/null +++ b/test/cell_lists/spatial_hashing.jl @@ -0,0 +1,89 @@ +@testset verbose=true "SpatialHashingCellList" begin + @testset "Collisions with foreach_neighbor" begin + #### Setup + # Rectangle of equidistantly spaced points + # from (x, y) = (-0.25, -0.25) to (x, y) = (0.35, 0.35). + range = -0.25:0.1:0.35 + coordinates1 = hcat(collect.(Iterators.product(range, range))...) + n_points = size(coordinates1, 2) + search_radius = 0.1 + + point_position1 = [0.05, 0.05] + point_index1 = findfirst(row -> row == point_position1, eachcol(coordinates1)) + + nhs1 = GridNeighborhoodSearch{2}(; search_radius, n_points, + cell_list = SpatialHashingCellList{2}(2 * n_points)) + initialize_grid!(nhs1, coordinates1) + + # Test if there is a collision at the point we test, else the test does not make much sense + cell = cell_coords(point_position1, nhs1) + index = spatial_hash(cell, nhs.cell_list.list_size) + @test nhs1.cell_list.cell_collision[index] = true + + function test_neighbor_function(point, neighbor, pos_diff, distance) + global found_neighbors1 + push!(found_neighbors1, neighbor) + end + + global found_neighbors1 = Int[] + foreach_neighbor(test_neighbor_function, coordinates1, coordinates1, nhs1, + point_index1) + + correct_neighbors1 = [17, 18, 19, 24, 25, 26, 31, 32, 33] + # @test correct_neighbors1 == found_neighbors1 + + # @infiltrate + + # # Get each neighbor for `point_position1` + # neighbors1 = sort(collect(PointNeighbors.eachneighbor(point_position1, nhs1))) + + # # Move points and update NHS + # coordinates2 = coordinates1 .+ [1.4, -3.5] + # empty!(nhs1.cell_list) + # initialize_grid!(nhs1, coordinates2) + + # # Get each neighbor for updated NHS + # neighbors2 = sort(collect(PointNeighbors.eachneighbor(point_position1, nhs1))) + + # # Change position + # point_position2 = point_position1 .+ [1.4, -3.5] + + # # Get each neighbor for `point_position2` + # neighbors3 = sort(collect(PointNeighbors.eachneighbor(point_position2, nhs1))) + + # # Double search radius + # nhs2 = GridNeighborhoodSearch{2}(search_radius = 2 * search_radius, + # n_points = size(coordinates1, 2), + # cell_list = SpatialHashingCellList{2}(2 * n_points)) + # initialize!(nhs2, coordinates1, coordinates1) + + # # Get each neighbor in double search radius + # neighbors4 = sort(collect(PointNeighbors.eachneighbor(point_position1, nhs2))) + + # # Move points + # coordinates2 = coordinates1 .+ [0.4, -0.4] + # empty!(nhs2.cell_list) + # initialize_grid!(nhs2, coordinates2) + + # # Get each neighbor in double search radius + # neighbors5 = sort(collect(PointNeighbors.eachneighbor(point_position1, nhs2))) + + # #### Verification against lists of potential neighbors built by hand + # correct_neighbors1 = [17, 18, 19, 24, 25, 26, 31, 32, 33] + # @test all(x -> x in neighbors1, correct_neighbors1) + + # correct_neighbors2 = Int[] + # @test all(x -> x in neighbors2, correct_neighbors2) + + # correct_neighbors3 = [17, 18, 19, 24, 25, 26, 31, 32, 33] + # @test all(x -> x in neighbors3, correct_neighbors3) + + # correct_neighbors4 = [ + # 9, 10, 11, 12, 13, 14, 16, 17, 18, 19, 20, 21, 23, 24, 25, 26, 27, 28, 30, 31, + # 32, 33, 34, 35, 37, 38, 39, 40, 41, 42, 44, 45, 46, 47, 48, 49] + # @test all(x -> x in neighbors4, correct_neighbors4) + + # correct_neighbors5 = [36, 37, 38, 43, 44, 45] + # @test all(x -> x in neighbors5, correct_neighbors5) + end +end; From c798f6fc4519840432d13c945241b9fc2a3de400 Mon Sep 17 00:00:00 2001 From: RubberLanding Date: Thu, 13 Mar 2025 17:44:03 +0800 Subject: [PATCH 06/27] Simple test for collisions with spatial hashing and foreach_neighbor(). --- src/cell_lists/spatial_hashing.jl | 10 ++-- src/nhs_grid.jl | 4 +- test/cell_lists/spatial_hashing.jl | 74 ++++-------------------------- 3 files changed, 15 insertions(+), 73 deletions(-) diff --git a/src/cell_lists/spatial_hashing.jl b/src/cell_lists/spatial_hashing.jl index eeeb1844..6db15dba 100644 --- a/src/cell_lists/spatial_hashing.jl +++ b/src/cell_lists/spatial_hashing.jl @@ -35,12 +35,12 @@ e{NDIMS, Real}, tuple of cell coordinates point::Int, index of the point to be added """ -@inline function cell_coords(coords, periodic_box::Nothing, cell_list::SpatialHashingCellList, - cell_size) -(; min_corner) = cell_list +# @inline function cell_coords(coords, periodic_box::Nothing, cell_list::SpatialHashingCellList, +# cell_size) +# (; min_corner) = cell_list -return Tuple(floor_to_int.((coords .- min_corner) ./ cell_size)) .+ 1 -end +# return Tuple(floor_to_int.((coords .- min_corner) ./ cell_size)) .+ 1 +# end function push_cell!(cell_list::SpatialHashingCellList, cell, point) diff --git a/src/nhs_grid.jl b/src/nhs_grid.jl index 5fc72958..e9a49662 100644 --- a/src/nhs_grid.jl +++ b/src/nhs_grid.jl @@ -367,9 +367,7 @@ function check_collision(neighbor_cell_::CartesianIndex, neighbor_coords, neighbor_cell_ = Tuple(neighbor_cell_) (; list_size, cell_collision) = cell_list - index = spatial_hash(neighbor_cell_, list_size) - collision = cell_collision[index] - if collision + if cell_collision[spatial_hash(neighbor_cell_, list_size)] # Check if `neighbor_coords` are in the cell `neighbor_cell_` cell_coords_ = cell_coords(neighbor_coords, nhs) return neighbor_cell_ != cell_coords_ diff --git a/test/cell_lists/spatial_hashing.jl b/test/cell_lists/spatial_hashing.jl index 86aef7b2..80313fd2 100644 --- a/test/cell_lists/spatial_hashing.jl +++ b/test/cell_lists/spatial_hashing.jl @@ -1,24 +1,22 @@ @testset verbose=true "SpatialHashingCellList" begin @testset "Collisions with foreach_neighbor" begin - #### Setup - # Rectangle of equidistantly spaced points - # from (x, y) = (-0.25, -0.25) to (x, y) = (0.35, 0.35). range = -0.25:0.1:0.35 coordinates1 = hcat(collect.(Iterators.product(range, range))...) n_points = size(coordinates1, 2) search_radius = 0.1 - point_position1 = [0.05, 0.05] - point_index1 = findfirst(row -> row == point_position1, eachcol(coordinates1)) + # point_position1 = [0.05, 0.05] + # point_index1 = findfirst(row -> row == point_position1, eachcol(coordinates1)) + point_index1 = 25 nhs1 = GridNeighborhoodSearch{2}(; search_radius, n_points, cell_list = SpatialHashingCellList{2}(2 * n_points)) initialize_grid!(nhs1, coordinates1) - # Test if there is a collision at the point we test, else the test does not make much sense - cell = cell_coords(point_position1, nhs1) - index = spatial_hash(cell, nhs.cell_list.list_size) - @test nhs1.cell_list.cell_collision[index] = true + # TODO: test if there is a collision at the point we test, else the test does not make much sense + # cell = cell_coords(point_position1, nhs1) + # index = spatial_hash(cell, nhs.cell_list.list_size) + # @test nhs1.cell_list.cell_collision[index] = true function test_neighbor_function(point, neighbor, pos_diff, distance) global found_neighbors1 @@ -29,61 +27,7 @@ foreach_neighbor(test_neighbor_function, coordinates1, coordinates1, nhs1, point_index1) - correct_neighbors1 = [17, 18, 19, 24, 25, 26, 31, 32, 33] - # @test correct_neighbors1 == found_neighbors1 - - # @infiltrate - - # # Get each neighbor for `point_position1` - # neighbors1 = sort(collect(PointNeighbors.eachneighbor(point_position1, nhs1))) - - # # Move points and update NHS - # coordinates2 = coordinates1 .+ [1.4, -3.5] - # empty!(nhs1.cell_list) - # initialize_grid!(nhs1, coordinates2) - - # # Get each neighbor for updated NHS - # neighbors2 = sort(collect(PointNeighbors.eachneighbor(point_position1, nhs1))) - - # # Change position - # point_position2 = point_position1 .+ [1.4, -3.5] - - # # Get each neighbor for `point_position2` - # neighbors3 = sort(collect(PointNeighbors.eachneighbor(point_position2, nhs1))) - - # # Double search radius - # nhs2 = GridNeighborhoodSearch{2}(search_radius = 2 * search_radius, - # n_points = size(coordinates1, 2), - # cell_list = SpatialHashingCellList{2}(2 * n_points)) - # initialize!(nhs2, coordinates1, coordinates1) - - # # Get each neighbor in double search radius - # neighbors4 = sort(collect(PointNeighbors.eachneighbor(point_position1, nhs2))) - - # # Move points - # coordinates2 = coordinates1 .+ [0.4, -0.4] - # empty!(nhs2.cell_list) - # initialize_grid!(nhs2, coordinates2) - - # # Get each neighbor in double search radius - # neighbors5 = sort(collect(PointNeighbors.eachneighbor(point_position1, nhs2))) - - # #### Verification against lists of potential neighbors built by hand - # correct_neighbors1 = [17, 18, 19, 24, 25, 26, 31, 32, 33] - # @test all(x -> x in neighbors1, correct_neighbors1) - - # correct_neighbors2 = Int[] - # @test all(x -> x in neighbors2, correct_neighbors2) - - # correct_neighbors3 = [17, 18, 19, 24, 25, 26, 31, 32, 33] - # @test all(x -> x in neighbors3, correct_neighbors3) - - # correct_neighbors4 = [ - # 9, 10, 11, 12, 13, 14, 16, 17, 18, 19, 20, 21, 23, 24, 25, 26, 27, 28, 30, 31, - # 32, 33, 34, 35, 37, 38, 39, 40, 41, 42, 44, 45, 46, 47, 48, 49] - # @test all(x -> x in neighbors4, correct_neighbors4) - - # correct_neighbors5 = [36, 37, 38, 43, 44, 45] - # @test all(x -> x in neighbors5, correct_neighbors5) + correct_neighbors1 = [18, 24, 25, 26, 32] + @test correct_neighbors1 == found_neighbors1 end end; From 614e68e524ef21a785170568b4e347a5ebef20fd Mon Sep 17 00:00:00 2001 From: RubberLanding Date: Thu, 13 Mar 2025 18:25:05 +0800 Subject: [PATCH 07/27] Add update!() for neighborhood_search tests. --- src/cell_lists/spatial_hashing.jl | 21 +++++++++++++++++---- src/nhs_grid.jl | 6 +++--- test/neighborhood_search.jl | 11 +++++++++-- 3 files changed, 29 insertions(+), 9 deletions(-) diff --git a/src/cell_lists/spatial_hashing.jl b/src/cell_lists/spatial_hashing.jl index 6db15dba..ddce591d 100644 --- a/src/cell_lists/spatial_hashing.jl +++ b/src/cell_lists/spatial_hashing.jl @@ -18,7 +18,7 @@ function SpatialHashingCellList{NDIMS}(list_size) where {NDIMS} # cell_choords is used to check if there is a collision and contains # the coordinates of the first cell added to the hash list - cell_coords = [missing for _ in 1:list_size] + cell_coords = [missing for _ in 1:list_size] return SpatialHashingCellList(list_size, cell_points, cell_coords, cell_collision) end @@ -42,7 +42,6 @@ point::Int, index of the point to be added # return Tuple(floor_to_int.((coords .- min_corner) ./ cell_size)) .+ 1 # end - function push_cell!(cell_list::SpatialHashingCellList, cell, point) (; cell_points, cell_coords, cell_collision, list_size) = cell_list key = spatial_hash(cell, list_size) @@ -51,12 +50,25 @@ function push_cell!(cell_list::SpatialHashingCellList, cell, point) # Check if the a cell has been added at this hash if ismissing(cell_coords) cell_coords = cell - # Detect collision + # Detect collision elseif cell != cell_coords cell_collision[key] = true end end +function update!(neighborhood_search::GridNeighborhoodSearch{<:Any, <:Any, <:SpatialHashingCellList}, x::AbstractMatrix, + y::AbstractMatrix; points_moving = (true, true), + parallelization_backend = x) + + # The coordinates of the first set of points are irrelevant for this NHS. + + # Only update when the second set is moving. + + points_moving[2] || return neighborhood_search + + initialize_grid!(neighborhood_search, y) +end + # Implement reset of collision flag, if after the deletion there still is no collision? function deleteat_cell!(cell_list::SpatialHashingCellList, cell, i) deleteat!(cell_list[cell], i) @@ -99,7 +111,8 @@ end function spatial_hash(cell::NTuple{3, Real}, index, list_size) i, j, k = cell - return mod(xor(i * 73856093, j * 19349663, k * 83492791, index * 7238423947), list_size) + 1 + return mod(xor(i * 73856093, j * 19349663, k * 83492791, index * 7238423947), + list_size) + 1 end function spatial_hash(cell::CartesianIndex{2}, list_size) diff --git a/src/nhs_grid.jl b/src/nhs_grid.jl index e9a49662..62f071ce 100644 --- a/src/nhs_grid.jl +++ b/src/nhs_grid.jl @@ -406,9 +406,9 @@ end pos_diff, distance2 = compute_periodic_distance(pos_diff, distance2, search_radius, periodic_box) - - is_near = distance2 <= search_radius^2 - print("$neighbor, $distance2, $is_near \n") + # DEBUG + # is_near = distance2 <= search_radius^2 + # print("$neighbor, $distance2, $is_near \n") if distance2 <= search_radius^2 distance = sqrt(distance2) diff --git a/test/neighborhood_search.jl b/test/neighborhood_search.jl index 673d551e..58c41c91 100644 --- a/test/neighborhood_search.jl +++ b/test/neighborhood_search.jl @@ -155,6 +155,7 @@ min_corner = minimum(coords, dims = 2) .- search_radius max_corner = maximum(coords, dims = 2) .+ search_radius + # Add spatial cell list nhs neighborhood_searches = [ GridNeighborhoodSearch{NDIMS}(; search_radius, n_points, update_strategy = SemiParallelUpdate()), @@ -175,6 +176,9 @@ max_corner, search_radius, backend = Vector{Vector{Int}})), + GridNeighborhoodSearch{NDIMS}(; search_radius, n_points, + cell_list = SpatialHashingCellList{2}(2 * + n_points)), PrecomputedNeighborhoodSearch{NDIMS}(; search_radius, n_points) ] @@ -184,8 +188,8 @@ "`GridNeighborhoodSearch` with `FullGridCellList` with `DynamicVectorOfVectors` and `ParallelUpdate`", "`GridNeighborhoodSearch` with `FullGridCellList` with `DynamicVectorOfVectors` and `SemiParallelUpdate`", "`GridNeighborhoodSearch` with `FullGridCellList` with `Vector{Vector}`", - "`PrecomputedNeighborhoodSearch`" - ] + "`GridNeighborhoodSearch` with `SpatialHashingCellList`", + "`PrecomputedNeighborhoodSearch`"] # Also test copied templates template_nhs = [ @@ -199,6 +203,9 @@ GridNeighborhoodSearch{NDIMS}(cell_list = FullGridCellList(; min_corner, max_corner, backend = Vector{Vector{Int32}})), + GridNeighborhoodSearch{NDIMS}(cell_list = FullGridCellList(; min_corner, + max_corner, + backend = Vector{Vector{Int32}})), PrecomputedNeighborhoodSearch{NDIMS}() ] copied_nhs = copy_neighborhood_search.(template_nhs, search_radius, n_points) From f6b2b12108259119a77776f829722dcbcdd66e99 Mon Sep 17 00:00:00 2001 From: RubberLanding Date: Fri, 21 Mar 2025 16:45:00 +0800 Subject: [PATCH 08/27] Running version for collision handling. Implement changes in spatial_hashing.jl to correct collision handling. Push collision check in __foreach_neighbor() for better performance. Extend test to test on multiple list sizes. Change plot.jl to test SpatialHashingCellList against the other data structures. --- benchmarks/plot.jl | 82 +++++++++++++++++++++++++----- src/PointNeighbors.jl | 2 - src/cell_lists/spatial_hashing.jl | 51 ++++++------------- src/nhs_grid.jl | 44 ++++++++-------- test/cell_lists/spatial_hashing.jl | 47 +++++++++-------- 5 files changed, 134 insertions(+), 92 deletions(-) diff --git a/benchmarks/plot.jl b/benchmarks/plot.jl index 976cf045..a16e4489 100644 --- a/benchmarks/plot.jl +++ b/benchmarks/plot.jl @@ -1,5 +1,6 @@ -using Plots +# using Plots using BenchmarkTools +# using Morton # Generate a rectangular point cloud include("../test/point_cloud.jl") @@ -38,9 +39,22 @@ plot_benchmarks(benchmark_count_neighbors, (10, 10), 3) function plot_benchmarks(benchmark, n_points_per_dimension, iterations; parallel = true, title = "", seed = 1, perturbation_factor_position = 1.0) - neighborhood_searches_names = ["TrivialNeighborhoodSearch";; - "GridNeighborhoodSearch";; - "PrecomputedNeighborhoodSearch"] + neighborhood_searches_names = [#"GridNeighborhoodSearch";; + #"GridNeighborhoodSearch with `PointWithCoordinates`";; + "GNHS with `DictionaryCellList`";; + "GNHS with `FullGridCellList`";; + "GNHS with `SpatialHashingCellList` n_points";; + "GNHS with `SpatialHashingCellList` 2 * n_points";; + "GNHS with `SpatialHashingCellList` 4 * n_points" + # "GNHS with `FullGridCellList` semi-parallel update";; + # "GNHS with `FullGridCellList` parallel incremental update";; + # "GNHS with `FullGridCellList` and `PointWithCoordinates`";;] + # "PrecomputedNeighborhoodSearch";; + # "CellListMap.jl";;] + # "CellListMap.jl with `points_equal_neighbors = true`";;] + # "original";;] + ] + # "z-index";;] # Multiply number of points in each iteration (roughly) by this factor scaling_factor = 4 @@ -54,18 +68,59 @@ function plot_benchmarks(benchmark, n_points_per_dimension, iterations; for iter in 1:iterations coordinates = point_cloud(sizes[iter], seed = seed, perturbation_factor_position = perturbation_factor_position) + # coords_z_index = point_cloud_z_index(sizes[iter], seed = seed, + # perturbation_factor_position = perturbation_factor_position) search_radius = 3.0 NDIMS = size(coordinates, 1) n_particles = size(coordinates, 2) + min_corner = Float32.(minimum(coordinates, dims = 2) .- search_radius) + max_corner = Float32.(maximum(coordinates, dims = 2) .+ search_radius) + neighborhood_searches = [ - TrivialNeighborhoodSearch{NDIMS}(; search_radius, eachpoint = 1:n_particles), - GridNeighborhoodSearch{NDIMS}(; search_radius, n_points = n_particles), - PrecomputedNeighborhoodSearch{NDIMS}(; search_radius, n_points = n_particles) + # TrivialNeighborhoodSearch{NDIMS}(; search_radius, eachpoint = 1:n_particles), + GridNeighborhoodSearch{NDIMS}(; search_radius, n_points = n_particles, + update_strategy = nothing), + GridNeighborhoodSearch{NDIMS}(; + cell_list = FullGridCellList(; search_radius, + min_corner, + max_corner), + search_radius, n_points = n_particles), + # GridNeighborhoodSearch{NDIMS}(; cell_list = DictionaryCellList{PointNeighbors.PointWithCoordinates{NDIMS, Float64}, NDIMS}(), + # search_radius, n_points = n_particles, update_strategy=nothing), + GridNeighborhoodSearch{NDIMS}(; search_radius, n_points = n_particles, + cell_list = SpatialHashingCellList{2}(1 * + n_particles)), + GridNeighborhoodSearch{NDIMS}(; search_radius, n_points = n_particles, + cell_list = SpatialHashingCellList{2}(2 * + n_particles)), + GridNeighborhoodSearch{NDIMS}(; search_radius, n_points = n_particles, + cell_list = SpatialHashingCellList{2}(4 * + n_particles)) + # GridNeighborhoodSearch{NDIMS}(; cell_list = FullGridCellList(; search_radius, + # min_corner, + # max_corner), + # search_radius, n_points = n_particles, + # update_strategy = SemiParallelUpdate()), + # GridNeighborhoodSearch{NDIMS}(; cell_list = FullGridCellList(; search_radius, + # min_corner, + # max_corner), + # search_radius, n_points = n_particles, + # update_strategy = ParallelIncrementalUpdate()), + # GridNeighborhoodSearch{NDIMS}(; cell_list = FullGridCellList(; search_radius, + # min_corner, + # max_corner, backend=PointNeighbors.DynamicVectorOfVectors{PointNeighbors.PointWithCoordinates{3, Float64}}), + # search_radius, n_points = n_particles, update_strategy=nothing), + # PrecomputedNeighborhoodSearch{NDIMS}(; search_radius, n_points = n_particles, update_strategy=nothing), + # CellListMapNeighborhoodSearch(NDIMS; search_radius), + # CellListMapNeighborhoodSearch(NDIMS; search_radius, points_equal_neighbors = true), ] for i in eachindex(neighborhood_searches) + # if i == 2 + # coordinates = coords_z_index + # end neighborhood_search = neighborhood_searches[i] initialize!(neighborhood_search, coordinates, coordinates) @@ -77,10 +132,11 @@ function plot_benchmarks(benchmark, n_points_per_dimension, iterations; end end - plot(n_particles_vec, times, - xaxis = :log, yaxis = :log, - xticks = (n_particles_vec, n_particles_vec), - xlabel = "#particles", ylabel = "Runtime [s]", - legend = :outerright, size = (750, 400), dpi = 200, - label = neighborhood_searches_names, title = title) + return n_particles_vec, times, neighborhood_searches_names + # plot(n_particles_vec, times, + # xaxis = :log, yaxis = :log, + # xticks = (n_particles_vec, n_particles_vec), + # xlabel = "#particles", ylabel = "Runtime [s]", + # legend = :outerright, size = (750, 400), dpi = 200, + # label = neighborhood_searches_names, title = title) end diff --git a/src/PointNeighbors.jl b/src/PointNeighbors.jl index bd6ba8f8..8596e513 100644 --- a/src/PointNeighbors.jl +++ b/src/PointNeighbors.jl @@ -11,8 +11,6 @@ using LinearAlgebra: dot using Polyester: Polyester @reexport using StaticArrays: SVector -using Infiltrator: @infiltrate - include("util.jl") include("vector_of_vectors.jl") include("neighborhood_search.jl") diff --git a/src/cell_lists/spatial_hashing.jl b/src/cell_lists/spatial_hashing.jl index ddce591d..f9d853d3 100644 --- a/src/cell_lists/spatial_hashing.jl +++ b/src/cell_lists/spatial_hashing.jl @@ -1,5 +1,6 @@ struct SpatialHashingCellList{CL, CI, CF} list_size::Int + NDIMS::Int cell_points::CL cell_coords::CI # Think about renaming this field, there is a function with the same name cell_collision::CF @@ -16,60 +17,40 @@ function SpatialHashingCellList{NDIMS}(list_size) where {NDIMS} cell_points = [Int[] for _ in 1:list_size] cell_collision = [false for _ in 1:list_size] - # cell_choords is used to check if there is a collision and contains + # Field cell_choords is used to check if there is a collision and contains # the coordinates of the first cell added to the hash list - cell_coords = [missing for _ in 1:list_size] - return SpatialHashingCellList(list_size, cell_points, cell_coords, cell_collision) + cell_coords = [ntuple(_ -> typemin(Int), NDIMS) for _ in 1:list_size] + return SpatialHashingCellList(list_size, NDIMS, cell_points, cell_coords, + cell_collision) end function Base.empty!(cell_list::SpatialHashingCellList) + (; list_size, NDIMS) = cell_list + Base.empty!.(cell_list.cell_points) - cell_list.cell_coords .= missing + cell_list.cell_coords .= [ntuple(_ -> typemin(Int), NDIMS) for _ in 1:list_size] cell_list.cell_collision .= false return cell_list end -""" -cell::NTupl -e{NDIMS, Real}, tuple of cell coordinates -point::Int, index of the point to be added -""" - -# @inline function cell_coords(coords, periodic_box::Nothing, cell_list::SpatialHashingCellList, -# cell_size) -# (; min_corner) = cell_list - -# return Tuple(floor_to_int.((coords .- min_corner) ./ cell_size)) .+ 1 -# end - function push_cell!(cell_list::SpatialHashingCellList, cell, point) (; cell_points, cell_coords, cell_collision, list_size) = cell_list key = spatial_hash(cell, list_size) - push!(cell_points[key], point) + cell_coord = cell_coords[key] + NDIMS = length(cell) + push!(cell_points[key], point) # Check if the a cell has been added at this hash - if ismissing(cell_coords) - cell_coords = cell - # Detect collision - elseif cell != cell_coords + if cell_coord == ntuple(_ -> typemin(Int), NDIMS) + cell_coords[key] = cell + # Detect collision + elseif cell_coord != cell cell_collision[key] = true end end -function update!(neighborhood_search::GridNeighborhoodSearch{<:Any, <:Any, <:SpatialHashingCellList}, x::AbstractMatrix, - y::AbstractMatrix; points_moving = (true, true), - parallelization_backend = x) - - # The coordinates of the first set of points are irrelevant for this NHS. - - # Only update when the second set is moving. - - points_moving[2] || return neighborhood_search - - initialize_grid!(neighborhood_search, y) -end - # Implement reset of collision flag, if after the deletion there still is no collision? +# Not needed if we do not use update!() but initialize!() function deleteat_cell!(cell_list::SpatialHashingCellList, cell, i) deleteat!(cell_list[cell], i) end diff --git a/src/nhs_grid.jl b/src/nhs_grid.jl index 62f071ce..1a255415 100644 --- a/src/nhs_grid.jl +++ b/src/nhs_grid.jl @@ -202,9 +202,21 @@ end # Update only with neighbor coordinates function update_grid!(neighborhood_search::GridNeighborhoodSearch{NDIMS}, - y::AbstractMatrix; parallelization_backend = y) where {NDIMS} - update_grid!(neighborhood_search, i -> extract_svector(y, Val(NDIMS), i); - parallelization_backend) + y::AbstractMatrix; parallelization_backend = y) where {NDIMS} +update_grid!(neighborhood_search, i -> extract_svector(y, Val(NDIMS), i); +parallelization_backend) +end + +# Needed here for SpatialHashingCellList until IncrementalParallelUpdate is implemented +function update!(neighborhood_search::GridNeighborhoodSearch{<:Any, <:Any, + <:SpatialHashingCellList}, + x::AbstractMatrix, + y::AbstractMatrix; points_moving = (true, true), + parallelization_backend = x) + # The coordinates of the first set of points are irrelevant for this NHS. + # Only update when the second set is moving. + points_moving[2] || return neighborhood_search + initialize_grid!(neighborhood_search, y) end # Serial and semi-parallel update. @@ -372,6 +384,7 @@ function check_collision(neighbor_cell_::CartesianIndex, neighbor_coords, cell_coords_ = cell_coords(neighbor_coords, nhs) return neighbor_cell_ != cell_coords_ end + return false end @inline function __foreach_neighbor(f, system_coords, neighbor_system_coords, @@ -392,27 +405,24 @@ end neighbor_coords = extract_svector(neighbor_system_coords, Val(ndims(neighborhood_search)), neighbor) - # Check if `neighbor_coords` are in the cell `neighbor_cell_`. - # For the `SpatialHashingCellList`, this might not be the case - # if we have a collision. - if check_collision(neighbor_cell_, neighbor_coords, cell_list, - neighborhood_search) - continue - end - pos_diff = point_coords - neighbor_coords distance2 = dot(pos_diff, pos_diff) pos_diff, distance2 = compute_periodic_distance(pos_diff, distance2, search_radius, periodic_box) - # DEBUG - # is_near = distance2 <= search_radius^2 - # print("$neighbor, $distance2, $is_near \n") if distance2 <= search_radius^2 distance = sqrt(distance2) + # Check if `neighbor_coords` are in the cell `neighbor_cell_`. + # For the `SpatialHashingCellList`, this might not be the case + # if we have a collision. + if 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) @@ -462,12 +472,6 @@ end return cell_coords(coords, periodic_box, cell_list, cell_size) end -@inline function cell_coords(coords, neighborhood_search) - (; periodic_box, cell_list, cell_size) = neighborhood_search - - return cell_coords(coords, periodic_box, cell_list, cell_size) -end - @inline function cell_coords(coords, periodic_box::Nothing, cell_list, cell_size) return Tuple(floor_to_int.(coords ./ cell_size)) end diff --git a/test/cell_lists/spatial_hashing.jl b/test/cell_lists/spatial_hashing.jl index 80313fd2..8a0693a2 100644 --- a/test/cell_lists/spatial_hashing.jl +++ b/test/cell_lists/spatial_hashing.jl @@ -1,33 +1,36 @@ @testset verbose=true "SpatialHashingCellList" begin @testset "Collisions with foreach_neighbor" begin - range = -0.25:0.1:0.35 + range = -0.35:0.1:0.25 coordinates1 = hcat(collect.(Iterators.product(range, range))...) n_points = size(coordinates1, 2) - search_radius = 0.1 - - # point_position1 = [0.05, 0.05] - # point_index1 = findfirst(row -> row == point_position1, eachcol(coordinates1)) + search_radius = 0.1 + 10 * eps() point_index1 = 25 + point_position1 = coordinates1[point_index1] - nhs1 = GridNeighborhoodSearch{2}(; search_radius, n_points, - cell_list = SpatialHashingCellList{2}(2 * n_points)) - initialize_grid!(nhs1, coordinates1) - - # TODO: test if there is a collision at the point we test, else the test does not make much sense - # cell = cell_coords(point_position1, nhs1) - # index = spatial_hash(cell, nhs.cell_list.list_size) - # @test nhs1.cell_list.cell_collision[index] = true + @testset verbose=true "List Size $(list_size)" for list_size in [ + 1, + n_points, + 2 * n_points, + 4 * n_points + ] + nhs1 = GridNeighborhoodSearch{2}(; search_radius, n_points, + cell_list = SpatialHashingCellList{2}(list_size)) + initialize_grid!(nhs1, coordinates1) - function test_neighbor_function(point, neighbor, pos_diff, distance) - global found_neighbors1 - push!(found_neighbors1, neighbor) - end + @testset verbose=true "Collision at test point" begin + cell = PointNeighbors.cell_coords(point_position1, nhs1) + index = PointNeighbors.spatial_hash(cell, nhs1.cell_list.list_size) + @test nhs1.cell_list.cell_collision[index] == true + end - global found_neighbors1 = Int[] - foreach_neighbor(test_neighbor_function, coordinates1, coordinates1, nhs1, - point_index1) + found_neighbors1 = Int[] + foreach_neighbor(coordinates1, coordinates1, nhs1, + point_index1) do point, neighbor, pos_diff, distance + push!(found_neighbors1, neighbor) + end - correct_neighbors1 = [18, 24, 25, 26, 32] - @test correct_neighbors1 == found_neighbors1 + correct_neighbors1 = [18, 24, 25, 26, 32] + @test correct_neighbors1 == found_neighbors1 + end end end; From c53ca75ec81d48f2aea95abf2023eef76ab3b55b Mon Sep 17 00:00:00 2001 From: RubberLanding Date: Wed, 9 Apr 2025 11:23:50 +0800 Subject: [PATCH 09/27] Fixed collision handling with empty cells. Clean up the code. --- benchmarks/plot.jl | 62 ++++++---------------- src/cell_lists/spatial_hashing.jl | 85 ++++++++++++++++-------------- src/nhs_grid.jl | 32 +++++++---- test/cell_lists/spatial_hashing.jl | 60 +++++++++++++++++++-- test/neighborhood_search.jl | 15 +++--- 5 files changed, 146 insertions(+), 108 deletions(-) diff --git a/benchmarks/plot.jl b/benchmarks/plot.jl index a16e4489..1a1bd83f 100644 --- a/benchmarks/plot.jl +++ b/benchmarks/plot.jl @@ -1,4 +1,4 @@ -# using Plots +using Plots using BenchmarkTools # using Morton @@ -36,25 +36,16 @@ include("benchmarks/benchmarks.jl") plot_benchmarks(benchmark_count_neighbors, (10, 10), 3) """ + function plot_benchmarks(benchmark, n_points_per_dimension, iterations; parallel = true, title = "", seed = 1, perturbation_factor_position = 1.0) - neighborhood_searches_names = [#"GridNeighborhoodSearch";; - #"GridNeighborhoodSearch with `PointWithCoordinates`";; - "GNHS with `DictionaryCellList`";; + neighborhood_searches_names = ["GNHS with `DictionaryCellList`";; "GNHS with `FullGridCellList`";; "GNHS with `SpatialHashingCellList` n_points";; "GNHS with `SpatialHashingCellList` 2 * n_points";; - "GNHS with `SpatialHashingCellList` 4 * n_points" - # "GNHS with `FullGridCellList` semi-parallel update";; - # "GNHS with `FullGridCellList` parallel incremental update";; - # "GNHS with `FullGridCellList` and `PointWithCoordinates`";;] - # "PrecomputedNeighborhoodSearch";; - # "CellListMap.jl";;] - # "CellListMap.jl with `points_equal_neighbors = true`";;] - # "original";;] + "GNHS with `SpatialHashingCellList` 4 * n_points";; ] - # "z-index";;] # Multiply number of points in each iteration (roughly) by this factor scaling_factor = 4 @@ -68,8 +59,6 @@ function plot_benchmarks(benchmark, n_points_per_dimension, iterations; for iter in 1:iterations coordinates = point_cloud(sizes[iter], seed = seed, perturbation_factor_position = perturbation_factor_position) - # coords_z_index = point_cloud_z_index(sizes[iter], seed = seed, - # perturbation_factor_position = perturbation_factor_position) search_radius = 3.0 NDIMS = size(coordinates, 1) @@ -79,7 +68,6 @@ function plot_benchmarks(benchmark, n_points_per_dimension, iterations; max_corner = Float32.(maximum(coordinates, dims = 2) .+ search_radius) neighborhood_searches = [ - # TrivialNeighborhoodSearch{NDIMS}(; search_radius, eachpoint = 1:n_particles), GridNeighborhoodSearch{NDIMS}(; search_radius, n_points = n_particles, update_strategy = nothing), GridNeighborhoodSearch{NDIMS}(; @@ -87,40 +75,18 @@ function plot_benchmarks(benchmark, n_points_per_dimension, iterations; min_corner, max_corner), search_radius, n_points = n_particles), - # GridNeighborhoodSearch{NDIMS}(; cell_list = DictionaryCellList{PointNeighbors.PointWithCoordinates{NDIMS, Float64}, NDIMS}(), - # search_radius, n_points = n_particles, update_strategy=nothing), GridNeighborhoodSearch{NDIMS}(; search_radius, n_points = n_particles, - cell_list = SpatialHashingCellList{2}(1 * + cell_list = SpatialHashingCellList{NDIMS}(1 * n_particles)), GridNeighborhoodSearch{NDIMS}(; search_radius, n_points = n_particles, - cell_list = SpatialHashingCellList{2}(2 * + cell_list = SpatialHashingCellList{NDIMS}(2 * n_particles)), GridNeighborhoodSearch{NDIMS}(; search_radius, n_points = n_particles, - cell_list = SpatialHashingCellList{2}(4 * + cell_list = SpatialHashingCellList{NDIMS}(4 * n_particles)) - # GridNeighborhoodSearch{NDIMS}(; cell_list = FullGridCellList(; search_radius, - # min_corner, - # max_corner), - # search_radius, n_points = n_particles, - # update_strategy = SemiParallelUpdate()), - # GridNeighborhoodSearch{NDIMS}(; cell_list = FullGridCellList(; search_radius, - # min_corner, - # max_corner), - # search_radius, n_points = n_particles, - # update_strategy = ParallelIncrementalUpdate()), - # GridNeighborhoodSearch{NDIMS}(; cell_list = FullGridCellList(; search_radius, - # min_corner, - # max_corner, backend=PointNeighbors.DynamicVectorOfVectors{PointNeighbors.PointWithCoordinates{3, Float64}}), - # search_radius, n_points = n_particles, update_strategy=nothing), - # PrecomputedNeighborhoodSearch{NDIMS}(; search_radius, n_points = n_particles, update_strategy=nothing), - # CellListMapNeighborhoodSearch(NDIMS; search_radius), - # CellListMapNeighborhoodSearch(NDIMS; search_radius, points_equal_neighbors = true), ] for i in eachindex(neighborhood_searches) - # if i == 2 - # coordinates = coords_z_index - # end neighborhood_search = neighborhood_searches[i] initialize!(neighborhood_search, coordinates, coordinates) @@ -132,11 +98,13 @@ function plot_benchmarks(benchmark, n_points_per_dimension, iterations; end end + p = plot(n_particles_vec, times, + xaxis = :log, yaxis = :log, + xticks = (n_particles_vec, n_particles_vec), + xlabel = "#particles", ylabel = "Runtime [s]", + legend = :outerright, size = (750, 400), dpi = 200, + label = neighborhood_searches_names, title = title) + display(p) + return n_particles_vec, times, neighborhood_searches_names - # plot(n_particles_vec, times, - # xaxis = :log, yaxis = :log, - # xticks = (n_particles_vec, n_particles_vec), - # xlabel = "#particles", ylabel = "Runtime [s]", - # legend = :outerright, size = (750, 400), dpi = 200, - # label = neighborhood_searches_names, title = title) end diff --git a/src/cell_lists/spatial_hashing.jl b/src/cell_lists/spatial_hashing.jl index f9d853d3..11a1dc61 100644 --- a/src/cell_lists/spatial_hashing.jl +++ b/src/cell_lists/spatial_hashing.jl @@ -1,76 +1,88 @@ struct SpatialHashingCellList{CL, CI, CF} list_size::Int NDIMS::Int - cell_points::CL - cell_coords::CI # Think about renaming this field, there is a function with the same name - cell_collision::CF + points::CL + coords::CI + collisions::CF end @inline index_type(::SpatialHashingCellList) = Int64 -# Change the order back to (SemiParallelUpdate, SerialUpdate) when SemiParallelUpdate is implemented -function supported_update_strategies(::SpatialHashingCellList) +function supported_update_strategies(::PointNeighbors.SpatialHashingCellList) return (SerialUpdate, SemiParallelUpdate) end function SpatialHashingCellList{NDIMS}(list_size) where {NDIMS} - cell_points = [Int[] for _ in 1:list_size] - cell_collision = [false for _ in 1:list_size] - - # Field cell_choords is used to check if there is a collision and contains - # the coordinates of the first cell added to the hash list - cell_coords = [ntuple(_ -> typemin(Int), NDIMS) for _ in 1:list_size] - return SpatialHashingCellList(list_size, NDIMS, cell_points, cell_coords, - cell_collision) + 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(list_size, NDIMS, points, coords, + collisions) end function Base.empty!(cell_list::SpatialHashingCellList) (; list_size, NDIMS) = cell_list - Base.empty!.(cell_list.cell_points) - cell_list.cell_coords .= [ntuple(_ -> typemin(Int), NDIMS) for _ in 1:list_size] - cell_list.cell_collision .= false + 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 function push_cell!(cell_list::SpatialHashingCellList, cell, point) - (; cell_points, cell_coords, cell_collision, list_size) = cell_list + (; points, coords, collisions, list_size, NDIMS) = cell_list key = spatial_hash(cell, list_size) - cell_coord = cell_coords[key] - NDIMS = length(cell) - - push!(cell_points[key], point) - # Check if the a cell has been added at this hash + cell_coord = coords[key] + push!(points[key], point) + # Check if a cell has been added at this hash if cell_coord == ntuple(_ -> typemin(Int), NDIMS) - cell_coords[key] = cell + coords[key] = cell # Detect collision elseif cell_coord != cell - cell_collision[key] = true + collisions[key] = true end end -# Implement reset of collision flag, if after the deletion there still is no collision? -# Not needed if we do not use update!() but initialize!() +using Infiltrator + +# function push_cell!(cell_list::SpatialHashingCellList, cell) +# (; coords, collisions, list_size, NDIMS) = cell_list +# key = spatial_hash(cell, list_size) +# cell_coord = coords[key] +# if cell_coord == ntuple(_ -> typemin(Int), NDIMS) +# coords[key] = cell +# elseif cell_coord != cell +# collisions[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.cell_points) +@inline each_cell_index(cell_list::SpatialHashingCellList) = eachindex(cell_list.points) @inline function Base.getindex(cell_list::SpatialHashingCellList, cell::Tuple) - (; cell_points) = cell_list + (; points) = cell_list - return cell_points[spatial_hash(cell, length(cell_points))] + return points[spatial_hash(cell, length(points))] end @inline function Base.getindex(cell_list::SpatialHashingCellList, i::Integer) - return cell_list.cell_points[i] + return cell_list.points[i] end -# Naming of cell_coords and cell_index confusing @inline function is_correct_cell(cell_list::SpatialHashingCellList{<:Any, Nothing}, - cell_coords, cell_index::Array) - return cell_coords == cell_index + coords, cell_index::Array) + return coords == cell_index +end + +function spatial_hash(cell::CartesianIndex{2}, list_size) + return spatial_hash(Tuple(cell), list_size) +end + +function spatial_hash(cell::CartesianIndex{3}, list_size) + return spatial_hash(Tuple(cell), list_size) end function spatial_hash(cell::NTuple{1, Real}, list_size) @@ -96,10 +108,3 @@ function spatial_hash(cell::NTuple{3, Real}, index, list_size) list_size) + 1 end -function spatial_hash(cell::CartesianIndex{2}, list_size) - return spatial_hash(Tuple(cell), list_size) -end - -function spatial_hash(cell::CartesianIndex{3}, list_size) - return spatial_hash(Tuple(cell), list_size) -end diff --git a/src/nhs_grid.jl b/src/nhs_grid.jl index 1a255415..b95c6280 100644 --- a/src/nhs_grid.jl +++ b/src/nhs_grid.jl @@ -202,9 +202,9 @@ end # Update only with neighbor coordinates function update_grid!(neighborhood_search::GridNeighborhoodSearch{NDIMS}, - y::AbstractMatrix; parallelization_backend = y) where {NDIMS} -update_grid!(neighborhood_search, i -> extract_svector(y, Val(NDIMS), i); -parallelization_backend) + y::AbstractMatrix; parallelization_backend = y) where {NDIMS} + update_grid!(neighborhood_search, i -> extract_svector(y, Val(NDIMS), i); + parallelization_backend) end # Needed here for SpatialHashingCellList until IncrementalParallelUpdate is implemented @@ -376,23 +376,35 @@ end function check_collision(neighbor_cell_::CartesianIndex, neighbor_coords, cell_list::SpatialHashingCellList, nhs) - neighbor_cell_ = Tuple(neighbor_cell_) - (; list_size, cell_collision) = cell_list - if cell_collision[spatial_hash(neighbor_cell_, list_size)] - # Check if `neighbor_coords` are in the cell `neighbor_cell_` - cell_coords_ = cell_coords(neighbor_coords, nhs) - return neighbor_cell_ != cell_coords_ + check_collision(Tuple(neighbor_cell_), neighbor_coords, cell_list, nhs) +end + +function check_collision(neighbor_cell_::Tuple, neighbor_coords, + cell_list::SpatialHashingCellList, nhs) + (; list_size, collisions, coords) = cell_list + + hash = spatial_hash(neighbor_cell_, list_size) + if collisions[hash] + # Points from multiple cells are in this list. + # Check if `neighbor_coords` are in the cell `neighbor_cell`. + return neighbor_cell_ != cell_coords(neighbor_coords, nhs) + elseif coords[hash] != neighbor_cell_ + # `cell_collision` is false for this list, meaning only points from one cell are in the list. + # However, the cell of this list is not our `neighbor_cell_`, so `neighbor_coords` + # are not in `neighbor_cell_`. + return true end return false end +using Infiltrator + @inline function __foreach_neighbor(f, system_coords, neighbor_system_coords, neighborhood_search::GridNeighborhoodSearch, point, point_coords, search_radius) (; 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) diff --git a/test/cell_lists/spatial_hashing.jl b/test/cell_lists/spatial_hashing.jl index 8a0693a2..296172ff 100644 --- a/test/cell_lists/spatial_hashing.jl +++ b/test/cell_lists/spatial_hashing.jl @@ -1,4 +1,55 @@ @testset verbose=true "SpatialHashingCellList" begin + @testset verbose=true "Compare Against `TrivialNeighborhoodSearch`" begin + cloud_sizes = [ + (10, 11), + (100, 90), + (3, 3, 3), + (39, 40, 41) + ] + + name(size) = "$(length(size))D with $(prod(size)) Particles" + @testset verbose=true "$(name(cloud_size))" for cloud_size in cloud_sizes + coords = point_cloud(cloud_size, seed = 1) + NDIMS = length(cloud_size) + n_points = size(coords, 2) + search_radius = 2.5 + + # Compute expected neighbor lists by brute-force looping over all points + # as potential neighbors (`TrivialNeighborhoodSearch`). + trivial_nhs = TrivialNeighborhoodSearch{NDIMS}(; search_radius, + eachpoint = axes(coords, 2)) + + neighbors_expected = [Int[] for _ in axes(coords, 2)] + + foreach_point_neighbor(coords, coords, trivial_nhs, + parallel = false) do point, neighbor, + pos_diff, distance + append!(neighbors_expected[point], neighbor) + end + + # Expand the domain by `search_radius`, as we need the neighboring cells of + # the minimum and maximum coordinates as well. + min_corner = minimum(coords, dims = 2) .- search_radius + max_corner = maximum(coords, dims = 2) .+ search_radius + + # Add spatial cell list nhs + nhs = GridNeighborhoodSearch{NDIMS}(; search_radius, n_points, + cell_list = SpatialHashingCellList{NDIMS}(2 * + n_points)) + + initialize!(nhs, coords, coords) + + neighbors = [Int[] for _ in axes(coords, 2)] + + foreach_point_neighbor(coords, coords, nhs, + parallel = false) do point, neighbor, + pos_diff, distance + push!(neighbors[point], neighbor) + end + @test sort.(neighbors) == neighbors_expected + end + end + @testset "Collisions with foreach_neighbor" begin range = -0.35:0.1:0.25 coordinates1 = hcat(collect.(Iterators.product(range, range))...) @@ -8,9 +59,9 @@ point_position1 = coordinates1[point_index1] @testset verbose=true "List Size $(list_size)" for list_size in [ + 2 * n_points, 1, n_points, - 2 * n_points, 4 * n_points ] nhs1 = GridNeighborhoodSearch{2}(; search_radius, n_points, @@ -20,7 +71,7 @@ @testset verbose=true "Collision at test point" begin cell = PointNeighbors.cell_coords(point_position1, nhs1) index = PointNeighbors.spatial_hash(cell, nhs1.cell_list.list_size) - @test nhs1.cell_list.cell_collision[index] == true + @test nhs1.cell_list.collisions[index] == true end found_neighbors1 = Int[] @@ -33,4 +84,7 @@ @test correct_neighbors1 == found_neighbors1 end end -end; + @testset "Collision Handling with empty cells" begin + #TODO: Add explicit test for list behavior with empty cells + end +end diff --git a/test/neighborhood_search.jl b/test/neighborhood_search.jl index 58c41c91..085c3346 100644 --- a/test/neighborhood_search.jl +++ b/test/neighborhood_search.jl @@ -155,7 +155,6 @@ min_corner = minimum(coords, dims = 2) .- search_radius max_corner = maximum(coords, dims = 2) .+ search_radius - # Add spatial cell list nhs neighborhood_searches = [ GridNeighborhoodSearch{NDIMS}(; search_radius, n_points, update_strategy = SemiParallelUpdate()), @@ -176,9 +175,9 @@ max_corner, search_radius, backend = Vector{Vector{Int}})), - GridNeighborhoodSearch{NDIMS}(; search_radius, n_points, - cell_list = SpatialHashingCellList{2}(2 * - n_points)), + # GridNeighborhoodSearch{NDIMS}(; search_radius, n_points, + # cell_list = SpatialHashingCellList{2}(2 * + # n_points)), PrecomputedNeighborhoodSearch{NDIMS}(; search_radius, n_points) ] @@ -188,7 +187,7 @@ "`GridNeighborhoodSearch` with `FullGridCellList` with `DynamicVectorOfVectors` and `ParallelUpdate`", "`GridNeighborhoodSearch` with `FullGridCellList` with `DynamicVectorOfVectors` and `SemiParallelUpdate`", "`GridNeighborhoodSearch` with `FullGridCellList` with `Vector{Vector}`", - "`GridNeighborhoodSearch` with `SpatialHashingCellList`", + # "`GridNeighborhoodSearch` with `SpatialHashingCellList`", "`PrecomputedNeighborhoodSearch`"] # Also test copied templates @@ -203,9 +202,9 @@ GridNeighborhoodSearch{NDIMS}(cell_list = FullGridCellList(; min_corner, max_corner, backend = Vector{Vector{Int32}})), - GridNeighborhoodSearch{NDIMS}(cell_list = FullGridCellList(; min_corner, - max_corner, - backend = Vector{Vector{Int32}})), + # GridNeighborhoodSearch{NDIMS}(cell_list = FullGridCellList(; min_corner, + # max_corner, + # backend = Vector{Vector{Int32}})), PrecomputedNeighborhoodSearch{NDIMS}() ] copied_nhs = copy_neighborhood_search.(template_nhs, search_radius, n_points) From f4fa219754f7d8780ac691486f0a8565f9fdd712 Mon Sep 17 00:00:00 2001 From: RubberLanding Date: Wed, 9 Apr 2025 11:27:24 +0800 Subject: [PATCH 10/27] Merge different tests for spatial hashing together. --- test/cell_lists/spatial_hashing.jl | 91 ++++++++++++++++-------------- 1 file changed, 50 insertions(+), 41 deletions(-) diff --git a/test/cell_lists/spatial_hashing.jl b/test/cell_lists/spatial_hashing.jl index 296172ff..5412e2b6 100644 --- a/test/cell_lists/spatial_hashing.jl +++ b/test/cell_lists/spatial_hashing.jl @@ -1,10 +1,13 @@ + +using Infiltrator @testset verbose=true "SpatialHashingCellList" begin + # General test for 2D and 3D @testset verbose=true "Compare Against `TrivialNeighborhoodSearch`" begin cloud_sizes = [ (10, 11), - (100, 90), - (3, 3, 3), - (39, 40, 41) + # (100, 90), + # (3, 3, 3), + # (39, 40, 41) ] name(size) = "$(length(size))D with $(prod(size)) Particles" @@ -27,18 +30,22 @@ append!(neighbors_expected[point], neighbor) end - # Expand the domain by `search_radius`, as we need the neighboring cells of - # the minimum and maximum coordinates as well. - min_corner = minimum(coords, dims = 2) .- search_radius - max_corner = maximum(coords, dims = 2) .+ search_radius - - # Add spatial cell list nhs nhs = GridNeighborhoodSearch{NDIMS}(; search_radius, n_points, cell_list = SpatialHashingCellList{NDIMS}(2 * n_points)) initialize!(nhs, coords, coords) + # Test if there are any collisions. + @testset verbose=true "Collision Detection" begin + @test any(nhs.cell_list.collisions) + end + + @testset verbose=true "Empty Cell and Collision Detection" begin + # 1. Find empty lists + # 2. Check for collisions with non-empty lists + end + neighbors = [Int[] for _ in axes(coords, 2)] foreach_point_neighbor(coords, coords, nhs, @@ -46,44 +53,46 @@ pos_diff, distance push!(neighbors[point], neighbor) end + @test sort.(neighbors) == neighbors_expected end end - @testset "Collisions with foreach_neighbor" begin - range = -0.35:0.1:0.25 - coordinates1 = hcat(collect.(Iterators.product(range, range))...) - n_points = size(coordinates1, 2) - search_radius = 0.1 + 10 * eps() - point_index1 = 25 - point_position1 = coordinates1[point_index1] - - @testset verbose=true "List Size $(list_size)" for list_size in [ - 2 * n_points, - 1, - n_points, - 4 * n_points - ] - nhs1 = GridNeighborhoodSearch{2}(; search_radius, n_points, - cell_list = SpatialHashingCellList{2}(list_size)) - initialize_grid!(nhs1, coordinates1) - - @testset verbose=true "Collision at test point" begin - cell = PointNeighbors.cell_coords(point_position1, nhs1) - index = PointNeighbors.spatial_hash(cell, nhs1.cell_list.list_size) - @test nhs1.cell_list.collisions[index] == true - end + # # Explicit test, how the data structure handels collisions (very simple) + # @testset "Collisions with foreach_neighbor" begin + # range = -0.35:0.1:0.25 + # coordinates1 = hcat(collect.(Iterators.product(range, range))...) + # n_points = size(coordinates1, 2) + # search_radius = 0.1 + 10 * eps() + # point_index1 = 25 + # point_position1 = coordinates1[point_index1] - found_neighbors1 = Int[] - foreach_neighbor(coordinates1, coordinates1, nhs1, - point_index1) do point, neighbor, pos_diff, distance - push!(found_neighbors1, neighbor) - end + # @testset verbose=true "List Size $(list_size)" for list_size in [ + # 2 * n_points, + # 1, + # n_points, + # 4 * n_points + # ] + # nhs1 = GridNeighborhoodSearch{2}(; search_radius, n_points, + # cell_list = SpatialHashingCellList{2}(list_size)) + # initialize_grid!(nhs1, coordinates1) - correct_neighbors1 = [18, 24, 25, 26, 32] - @test correct_neighbors1 == found_neighbors1 - end - end + # @testset verbose=true "Collision at test point" begin + # cell = PointNeighbors.cell_coords(point_position1, nhs1) + # index = PointNeighbors.spatial_hash(cell, nhs1.cell_list.list_size) + # @test nhs1.cell_list.collisions[index] == true + # end + + # found_neighbors1 = Int[] + # foreach_neighbor(coordinates1, coordinates1, nhs1, + # point_index1) do point, neighbor, pos_diff, distance + # push!(found_neighbors1, neighbor) + # end + + # correct_neighbors1 = [18, 24, 25, 26, 32] + # @test correct_neighbors1 == found_neighbors1 + # end + # end @testset "Collision Handling with empty cells" begin #TODO: Add explicit test for list behavior with empty cells end From ca1054bc711160cec5549683cb28ae0e72f26edd Mon Sep 17 00:00:00 2001 From: RubberLanding Date: Thu, 10 Apr 2025 12:00:58 +0800 Subject: [PATCH 11/27] Clean up. --- src/cell_lists/spatial_hashing.jl | 2 -- src/nhs_grid.jl | 2 -- test/cell_lists/spatial_hashing.jl | 2 -- 3 files changed, 6 deletions(-) diff --git a/src/cell_lists/spatial_hashing.jl b/src/cell_lists/spatial_hashing.jl index 11a1dc61..e74f66f5 100644 --- a/src/cell_lists/spatial_hashing.jl +++ b/src/cell_lists/spatial_hashing.jl @@ -43,8 +43,6 @@ function push_cell!(cell_list::SpatialHashingCellList, cell, point) end end -using Infiltrator - # function push_cell!(cell_list::SpatialHashingCellList, cell) # (; coords, collisions, list_size, NDIMS) = cell_list # key = spatial_hash(cell, list_size) diff --git a/src/nhs_grid.jl b/src/nhs_grid.jl index c4dbac4c..1cb3652d 100644 --- a/src/nhs_grid.jl +++ b/src/nhs_grid.jl @@ -410,8 +410,6 @@ function check_collision(neighbor_cell_::Tuple, neighbor_coords, return false end -using Infiltrator - @inline function __foreach_neighbor(f, system_coords, neighbor_system_coords, neighborhood_search::GridNeighborhoodSearch, point, point_coords, search_radius) diff --git a/test/cell_lists/spatial_hashing.jl b/test/cell_lists/spatial_hashing.jl index 5412e2b6..893ab566 100644 --- a/test/cell_lists/spatial_hashing.jl +++ b/test/cell_lists/spatial_hashing.jl @@ -1,5 +1,3 @@ - -using Infiltrator @testset verbose=true "SpatialHashingCellList" begin # General test for 2D and 3D @testset verbose=true "Compare Against `TrivialNeighborhoodSearch`" begin From 260850cba94d60c840f4e84e250a17b5e9b5a0ee Mon Sep 17 00:00:00 2001 From: RubberLanding Date: Thu, 10 Apr 2025 15:58:09 +0800 Subject: [PATCH 12/27] Add documentation. --- src/cell_lists/spatial_hashing.jl | 22 +++++++++++++++------- test/cell_lists/spatial_hashing.jl | 18 +++++++++--------- 2 files changed, 24 insertions(+), 16 deletions(-) diff --git a/src/cell_lists/spatial_hashing.jl b/src/cell_lists/spatial_hashing.jl index e74f66f5..0286b26f 100644 --- a/src/cell_lists/spatial_hashing.jl +++ b/src/cell_lists/spatial_hashing.jl @@ -1,3 +1,16 @@ +""" + 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. (2001)](@cite Ihmsen2003). By using a hash map, which only stores non-empty cells, +the domain is effectively infinite. The size of the hash map is recommended to be approximately twice the number of particles balance memory consumption +and 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{CL, CI, CF} list_size::Int NDIMS::Int @@ -6,7 +19,7 @@ struct SpatialHashingCellList{CL, CI, CF} collisions::CF end -@inline index_type(::SpatialHashingCellList) = Int64 +@inline index_type(::SpatialHashingCellList) = Int32 function supported_update_strategies(::PointNeighbors.SpatialHashingCellList) return (SerialUpdate, SemiParallelUpdate) @@ -62,7 +75,6 @@ end @inline function Base.getindex(cell_list::SpatialHashingCellList, cell::Tuple) (; points) = cell_list - return points[spatial_hash(cell, length(points))] end @@ -75,11 +87,7 @@ end return coords == cell_index end -function spatial_hash(cell::CartesianIndex{2}, list_size) - return spatial_hash(Tuple(cell), list_size) -end - -function spatial_hash(cell::CartesianIndex{3}, list_size) +function spatial_hash(cell::CartesianIndex, list_size) return spatial_hash(Tuple(cell), list_size) end diff --git a/test/cell_lists/spatial_hashing.jl b/test/cell_lists/spatial_hashing.jl index 893ab566..aca90180 100644 --- a/test/cell_lists/spatial_hashing.jl +++ b/test/cell_lists/spatial_hashing.jl @@ -2,9 +2,9 @@ # General test for 2D and 3D @testset verbose=true "Compare Against `TrivialNeighborhoodSearch`" begin cloud_sizes = [ - (10, 11), + # (10, 11), # (100, 90), - # (3, 3, 3), + (8, 10, 6), # (39, 40, 41) ] @@ -39,10 +39,10 @@ @test any(nhs.cell_list.collisions) end - @testset verbose=true "Empty Cell and Collision Detection" begin - # 1. Find empty lists - # 2. Check for collisions with non-empty lists - end + # @testset verbose=true "Empty Cell and Collision Detection" begin + # # 1. Find empty lists + # # 2. Check for collisions with non-empty lists + # end neighbors = [Int[] for _ in axes(coords, 2)] @@ -91,7 +91,7 @@ # @test correct_neighbors1 == found_neighbors1 # end # end - @testset "Collision Handling with empty cells" begin - #TODO: Add explicit test for list behavior with empty cells - end + # @testset "Collision Handling with empty cells" begin + # #TODO: Add explicit test for list behavior with empty cells + # end end From d79dc256a4a11733916d5bdb2b64f9de62170c34 Mon Sep 17 00:00:00 2001 From: RubberLanding Date: Thu, 10 Apr 2025 16:20:20 +0800 Subject: [PATCH 13/27] Clean up. --- src/nhs_spatial_hash.jl | 91 ----------------------------------------- 1 file changed, 91 deletions(-) delete mode 100644 src/nhs_spatial_hash.jl diff --git a/src/nhs_spatial_hash.jl b/src/nhs_spatial_hash.jl deleted file mode 100644 index 0422d32c..00000000 --- a/src/nhs_spatial_hash.jl +++ /dev/null @@ -1,91 +0,0 @@ -struct SpatialHashingCellList{CL, CI, CF} - cell_points :: CL - cell_coords :: CI - cell_collision :: CF - size :: Int -end - -function SpatialHashingCellList{NDIMS}(size; handle_collisions = true) where {NDIMS} - cell_points = [Int[] for _ in 1:size] - cell_collision = [false for _ in 1:size] - cell_coords = [Int[] for _ in 1:size] - return SpatialHashingCellList(cell_points, cell_coords, cell_collision, size) -end - -function Base.empty!(cell_list::SpatialHashingCellList) - Base.empty!.(cell_list.cell_points) - Base.empty!.(cell_list.cell_coords) - cell_list.cell_collision .= false - - return cell_list -end - -function push_cell!(cell_list::SpatialHashingCellList, cell, point) - (; cell_points, cell_coords, cell_collision, size) = cell_list - key = spatial_hash(cell, size) - push!(cell_points[key], point) - - if !(cell in cell_coords[key]) - push!(cell_coords[key], cell) - end - - if length(cell_coords[key]) > 1 - cell_collision[key] = true - end -end - -function get_points(nhs, cell, coords_fun) - (; cell_list) = nhs - (; cell_points, cell_coords, cell_collision, size) = cell_list - key = spatial_hash(cell, size) - if cell_collision[key] - points_in_cell = [] - for point in cell_list.cell_points[key] - if cell_coords(coords_fun(point), nhs) == cell - push!(points_in_cell, point) - end - end - return points_in_cell - - else - return cell_list.cell_points[key] - 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.cell_points) - -@inline function Base.getindex(cell_list::SpatialHashingCellList, i::Integer) - return cell_list.cell_points[i] -end - -@inline function is_correct_cell(cell_list::SpatialHashingCellList{<:Any, Nothing}, - cell_coords, cell_index::Array) # What is the correct type for cell_index? It is a list of integers -return cell_coords in cell_list.cell_coords[cell_index] -end - -function spatial_hash(cell::NTuple{1, Real}, size) - return mod(cell[1] * 73856093, size) + 1 -end - -function spatial_hash(cell::NTuple{2, Real}, size::UInt64) - i, j = cell - - return mod(xor(i * 73856093, j * 19349663), size) + 1 -end - -function spatial_hash(cell::NTuple{3, Real}, size) - i, j, k = cell - - return mod(xor(i * 73856093, j * 19349663, k * 83492791), size) + 1 -end - -function spatial_hash(cell::NTuple{3, Real}, index, size) - i, j, k = cell - - return mod(xor(i * 73856093, j * 19349663, k * 83492791, index * 7238423947), size) + 1 -end From 4b96172c52438e723749bb1697c99769a2c71021 Mon Sep 17 00:00:00 2001 From: RubberLanding Date: Thu, 17 Apr 2025 15:31:40 +0800 Subject: [PATCH 14/27] Update test for spatial hashing, add test for empty list. --- test/cell_lists/spatial_hashing.jl | 65 ++++++++++++------------------ 1 file changed, 25 insertions(+), 40 deletions(-) diff --git a/test/cell_lists/spatial_hashing.jl b/test/cell_lists/spatial_hashing.jl index aca90180..4127a882 100644 --- a/test/cell_lists/spatial_hashing.jl +++ b/test/cell_lists/spatial_hashing.jl @@ -2,10 +2,10 @@ # General test for 2D and 3D @testset verbose=true "Compare Against `TrivialNeighborhoodSearch`" begin cloud_sizes = [ - # (10, 11), - # (100, 90), + (10, 11), + (100, 90), (8, 10, 6), - # (39, 40, 41) + (39, 40, 41) ] name(size) = "$(length(size))D with $(prod(size)) Particles" @@ -56,42 +56,27 @@ end end - # # Explicit test, how the data structure handels collisions (very simple) - # @testset "Collisions with foreach_neighbor" begin - # range = -0.35:0.1:0.25 - # coordinates1 = hcat(collect.(Iterators.product(range, range))...) - # n_points = size(coordinates1, 2) - # search_radius = 0.1 + 10 * eps() - # point_index1 = 25 - # point_position1 = coordinates1[point_index1] - - # @testset verbose=true "List Size $(list_size)" for list_size in [ - # 2 * n_points, - # 1, - # n_points, - # 4 * n_points - # ] - # nhs1 = GridNeighborhoodSearch{2}(; search_radius, n_points, - # cell_list = SpatialHashingCellList{2}(list_size)) - # initialize_grid!(nhs1, coordinates1) - - # @testset verbose=true "Collision at test point" begin - # cell = PointNeighbors.cell_coords(point_position1, nhs1) - # index = PointNeighbors.spatial_hash(cell, nhs1.cell_list.list_size) - # @test nhs1.cell_list.collisions[index] == true - # end - - # found_neighbors1 = Int[] - # foreach_neighbor(coordinates1, coordinates1, nhs1, - # point_index1) do point, neighbor, pos_diff, distance - # push!(found_neighbors1, neighbor) - # end + # Test list behavior with empty cells + @testset "Collision Handling with empty cells" begin + # The point is in cell (-1, 0) which has a hash collision with cell (-2, -1) + point = [-0.05, 0.05] + coordinates = hcat(point) + 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) + + found_neighbors = Int[] + foreach_neighbor(coordinates, coordinates, nhs, + point_index) do point, neighbor, pos_diff, distance + push!(found_neighbors, neighbor) + end - # correct_neighbors1 = [18, 24, 25, 26, 32] - # @test correct_neighbors1 == found_neighbors1 - # end - # end - # @testset "Collision Handling with empty cells" begin - # #TODO: Add explicit test for list behavior with empty cells - # end + correct_neighbors = [1] + @test correct_neighbors == found_neighbors + end end From ed74323fd113a2f12703b6b8e215e9dd4a0fc7bc Mon Sep 17 00:00:00 2001 From: RubberLanding Date: Thu, 17 Apr 2025 15:33:10 +0800 Subject: [PATCH 15/27] Add copy_cell_list() for spatial hashing. Include SpatialHashingCellList in test/neighborhood_search.jl. --- src/cell_lists/spatial_hashing.jl | 32 ++++++++++++----------------- test/neighborhood_search.jl | 34 +++++++++++++++++++------------ 2 files changed, 34 insertions(+), 32 deletions(-) diff --git a/src/cell_lists/spatial_hashing.jl b/src/cell_lists/spatial_hashing.jl index 0286b26f..f24a8454 100644 --- a/src/cell_lists/spatial_hashing.jl +++ b/src/cell_lists/spatial_hashing.jl @@ -42,37 +42,32 @@ function Base.empty!(cell_list::SpatialHashingCellList) 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, NDIMS) = cell_list - key = spatial_hash(cell, list_size) - cell_coord = coords[key] - push!(points[key], point) - # Check if a cell has been added at this hash + hash_key = spatial_hash(cell, list_size) + cell_coord = coords[hash_key] + push!(points[hash_key], point) if cell_coord == ntuple(_ -> typemin(Int), NDIMS) - coords[key] = cell - # Detect collision + coords[hash_key] = cell elseif cell_coord != cell - collisions[key] = true + collisions[hash_key] = true end end -# function push_cell!(cell_list::SpatialHashingCellList, cell) -# (; coords, collisions, list_size, NDIMS) = cell_list -# key = spatial_hash(cell, list_size) -# cell_coord = coords[key] -# if cell_coord == ntuple(_ -> typemin(Int), NDIMS) -# coords[key] = cell -# elseif cell_coord != cell -# collisions[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) + (; NDIMS, list_size) = cell_list + return SpatialHashingCellList{NDIMS}(list_size) +end + @inline function Base.getindex(cell_list::SpatialHashingCellList, cell::Tuple) (; points) = cell_list return points[spatial_hash(cell, length(points))] @@ -113,4 +108,3 @@ function spatial_hash(cell::NTuple{3, Real}, index, list_size) return mod(xor(i * 73856093, j * 19349663, k * 83492791, index * 7238423947), list_size) + 1 end - diff --git a/test/neighborhood_search.jl b/test/neighborhood_search.jl index 085c3346..30a9d65a 100644 --- a/test/neighborhood_search.jl +++ b/test/neighborhood_search.jl @@ -1,3 +1,5 @@ +using Infiltrator + # This file contains tests for the generic functions in `src/neighborhood_search.jl` and # tests comparing all NHS implementations against the `TrivialNeighborhoodSearch`. @testset verbose=true "All Neighborhood Searches" begin @@ -58,7 +60,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 +72,8 @@ "`GridNeighborhoodSearch`", "`GridNeighborhoodSearch` with `FullGridCellList` with `DynamicVectorOfVectors`", "`GridNeighborhoodSearch` with `FullGridCellList` with `Vector{Vector}`", - "`PrecomputedNeighborhoodSearch`" + "`PrecomputedNeighborhoodSearch`", + "`GridNeighborhoodSearch` with `SpatialHashingCellList`" ] # Also test copied templates @@ -80,7 +87,9 @@ 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) @@ -175,10 +184,10 @@ max_corner, search_radius, backend = Vector{Vector{Int}})), - # GridNeighborhoodSearch{NDIMS}(; search_radius, n_points, - # cell_list = SpatialHashingCellList{2}(2 * - # n_points)), - 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 = [ @@ -187,8 +196,9 @@ "`GridNeighborhoodSearch` with `FullGridCellList` with `DynamicVectorOfVectors` and `ParallelUpdate`", "`GridNeighborhoodSearch` with `FullGridCellList` with `DynamicVectorOfVectors` and `SemiParallelUpdate`", "`GridNeighborhoodSearch` with `FullGridCellList` with `Vector{Vector}`", - # "`GridNeighborhoodSearch` with `SpatialHashingCellList`", - "`PrecomputedNeighborhoodSearch`"] + "`PrecomputedNeighborhoodSearch`", + "`GridNeighborhoodSearch` with `SpatialHashingCellList`", + ] # Also test copied templates template_nhs = [ @@ -202,10 +212,8 @@ GridNeighborhoodSearch{NDIMS}(cell_list = FullGridCellList(; min_corner, max_corner, backend = Vector{Vector{Int32}})), - # 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) From 41bb8710e175de0dbec1900207763faf154676ae Mon Sep 17 00:00:00 2001 From: RubberLanding Date: Thu, 17 Apr 2025 15:34:30 +0800 Subject: [PATCH 16/27] Update check_collision(), the tests are succeeding. Add check_cell_collision() for benchmarking. --- src/nhs_grid.jl | 52 ++++++++++++++++++++++++++++++++++++------------- 1 file changed, 38 insertions(+), 14 deletions(-) diff --git a/src/nhs_grid.jl b/src/nhs_grid.jl index 1cb3652d..c63e0777 100644 --- a/src/nhs_grid.jl +++ b/src/nhs_grid.jl @@ -388,28 +388,50 @@ end function check_collision(neighbor_cell_::CartesianIndex, neighbor_coords, cell_list::SpatialHashingCellList, nhs) - - check_collision(Tuple(neighbor_cell_), neighbor_coords, cell_list, nhs) -end - -function check_collision(neighbor_cell_::Tuple, neighbor_coords, - cell_list::SpatialHashingCellList, nhs) (; list_size, collisions, coords) = cell_list + # Convert CartesianIndex to periodic index + neighbor_cell = periodic_cell_index(Tuple(neighbor_cell_), nhs) + hash_key = spatial_hash(neighbor_cell, list_size) - hash = spatial_hash(neighbor_cell_, list_size) - if collisions[hash] + if collisions[hash_key] # Points from multiple cells are in this list. # Check if `neighbor_coords` are in the cell `neighbor_cell`. - return neighbor_cell_ != cell_coords(neighbor_coords, nhs) - elseif coords[hash] != neighbor_cell_ + return neighbor_cell != cell_coords(neighbor_coords, nhs) + elseif coords[hash_key] != neighbor_cell # `cell_collision` is false for this list, meaning only points from one cell are in the list. - # However, the cell of this list is not our `neighbor_cell_`, so `neighbor_coords` - # are not in `neighbor_cell_`. + # However, the cell of this list is not our `neighbor_cell`, so `neighbor_coords` + # are not in `neighbor_cell`. return true end return false end +function check_collision1(neighbor_cell_, neighbor_coords, cell_list, nhs) + return false +end + +function check_collision1(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) + return false +end + +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) + + return collisions[hash] || coords[hash] != neighbor_cell +end + @inline function __foreach_neighbor(f, system_coords, neighbor_system_coords, neighborhood_search::GridNeighborhoodSearch, point, point_coords, search_radius) @@ -419,6 +441,9 @@ end neighbor_cell = Tuple(neighbor_cell_) neighbors = points_in_cell(neighbor_cell, neighborhood_search) + cell_collision = check_cell_collision(neighbor_cell_, + cell_list, neighborhood_search) + for neighbor_ in eachindex(neighbors) neighbor = @inbounds neighbors[neighbor_] @@ -440,8 +465,7 @@ end # Check if `neighbor_coords` are in the cell `neighbor_cell_`. # For the `SpatialHashingCellList`, this might not be the case # if we have a collision. - if check_collision(neighbor_cell_, neighbor_coords, cell_list, - neighborhood_search) + if check_collision1(neighbor_cell_, neighbor_coords, cell_list, neighborhood_search) && cell_collision continue end From 00822ebbcf3dab108b7bec670dd6b73b7a68f8f2 Mon Sep 17 00:00:00 2001 From: RubberLanding Date: Thu, 17 Apr 2025 19:02:50 +0800 Subject: [PATCH 17/27] Clean up. --- src/cell_lists/spatial_hashing.jl | 7 +++---- src/nhs_grid.jl | 12 ------------ 2 files changed, 3 insertions(+), 16 deletions(-) diff --git a/src/cell_lists/spatial_hashing.jl b/src/cell_lists/spatial_hashing.jl index f24a8454..1683901a 100644 --- a/src/cell_lists/spatial_hashing.jl +++ b/src/cell_lists/spatial_hashing.jl @@ -12,11 +12,11 @@ and the likelihood of hash collisions. - `list_size::Int`: Size of the hash map (e.g., `2 * n_points`) . """ struct SpatialHashingCellList{CL, CI, CF} - list_size::Int - NDIMS::Int points::CL coords::CI collisions::CF + list_size::Int + NDIMS::Int end @inline index_type(::SpatialHashingCellList) = Int32 @@ -29,8 +29,7 @@ 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(list_size, NDIMS, points, coords, - collisions) + return SpatialHashingCellList(points, coords, collisions, list_size, NDIMS) end function Base.empty!(cell_list::SpatialHashingCellList) diff --git a/src/nhs_grid.jl b/src/nhs_grid.jl index c63e0777..d4e890ec 100644 --- a/src/nhs_grid.jl +++ b/src/nhs_grid.jl @@ -210,18 +210,6 @@ function update_grid!(neighborhood_search::GridNeighborhoodSearch{NDIMS}, parallelization_backend) end -# Needed here for SpatialHashingCellList until IncrementalParallelUpdate is implemented -function update!(neighborhood_search::GridNeighborhoodSearch{<:Any, <:Any, - <:SpatialHashingCellList}, - x::AbstractMatrix, - y::AbstractMatrix; points_moving = (true, true), - parallelization_backend = x) - # The coordinates of the first set of points are irrelevant for this NHS. - # Only update when the second set is moving. - points_moving[2] || return neighborhood_search - initialize_grid!(neighborhood_search, y) -end - # Serial and semi-parallel update. # See the warning above. `parallelization_backend = nothing` will use `Polyester.@batch`. function update_grid!(neighborhood_search::Union{GridNeighborhoodSearch{<:Any, From 81c640c811b7efa5394856b3fe61a08ad8fd6840 Mon Sep 17 00:00:00 2001 From: RubberLanding Date: Thu, 17 Apr 2025 19:10:06 +0800 Subject: [PATCH 18/27] Clean up. --- src/nhs_grid.jl | 26 +------------------------- 1 file changed, 1 insertion(+), 25 deletions(-) diff --git a/src/nhs_grid.jl b/src/nhs_grid.jl index d4e890ec..a6233a5b 100644 --- a/src/nhs_grid.jl +++ b/src/nhs_grid.jl @@ -377,30 +377,6 @@ end function check_collision(neighbor_cell_::CartesianIndex, neighbor_coords, cell_list::SpatialHashingCellList, nhs) (; list_size, collisions, coords) = cell_list - # Convert CartesianIndex to periodic index - neighbor_cell = periodic_cell_index(Tuple(neighbor_cell_), nhs) - hash_key = spatial_hash(neighbor_cell, list_size) - - if collisions[hash_key] - # Points from multiple cells are in this list. - # Check if `neighbor_coords` are in the cell `neighbor_cell`. - return neighbor_cell != cell_coords(neighbor_coords, nhs) - elseif coords[hash_key] != neighbor_cell - # `cell_collision` is false for this list, meaning only points from one cell are in the list. - # However, the cell of this list is not our `neighbor_cell`, so `neighbor_coords` - # are not in `neighbor_cell`. - return true - end - return false -end - -function check_collision1(neighbor_cell_, neighbor_coords, cell_list, nhs) - return false -end - -function check_collision1(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) @@ -453,7 +429,7 @@ end # Check if `neighbor_coords` are in the cell `neighbor_cell_`. # For the `SpatialHashingCellList`, this might not be the case # if we have a collision. - if check_collision1(neighbor_cell_, neighbor_coords, cell_list, neighborhood_search) && cell_collision + if check_collision(neighbor_cell_, neighbor_coords, cell_list, neighborhood_search) && cell_collision continue end From 48f9651ec24e4ae4036901b53809116b77b0cc2d Mon Sep 17 00:00:00 2001 From: RubberLanding Date: Thu, 24 Apr 2025 13:02:32 +0800 Subject: [PATCH 19/27] Fix code in foreach_neighbor after merge. --- src/cell_lists/spatial_hashing.jl | 8 +++++--- src/nhs_grid.jl | 4 ++-- test/cell_lists/spatial_hashing.jl | 7 ------- test/neighborhood_search.jl | 2 -- 4 files changed, 7 insertions(+), 14 deletions(-) diff --git a/src/cell_lists/spatial_hashing.jl b/src/cell_lists/spatial_hashing.jl index 1683901a..4d83d46c 100644 --- a/src/cell_lists/spatial_hashing.jl +++ b/src/cell_lists/spatial_hashing.jl @@ -11,7 +11,8 @@ and the likelihood of hash collisions. - `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{CL, CI, CF} + +struct SpatialHashingCellList{CL, CI, CF} <: AbstractCellList points::CL coords::CI collisions::CF @@ -21,8 +22,8 @@ end @inline index_type(::SpatialHashingCellList) = Int32 -function supported_update_strategies(::PointNeighbors.SpatialHashingCellList) - return (SerialUpdate, SemiParallelUpdate) +function supported_update_strategies(::SpatialHashingCellList) + return (SerialUpdate,) end function SpatialHashingCellList{NDIMS}(list_size) where {NDIMS} @@ -45,6 +46,7 @@ end # 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, NDIMS) = cell_list + hash_key = spatial_hash(cell, list_size) cell_coord = coords[hash_key] push!(points[hash_key], point) diff --git a/src/nhs_grid.jl b/src/nhs_grid.jl index 76713a0d..2393b0cc 100644 --- a/src/nhs_grid.jl +++ b/src/nhs_grid.jl @@ -469,15 +469,15 @@ 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, +@inline function foreach_neighbor(f, neighbor_system_coords, neighborhood_search::GridNeighborhoodSearch, point, point_coords, search_radius) (; 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) - cell_collision = check_cell_collision(neighbor_cell_, cell_list, neighborhood_search) diff --git a/test/cell_lists/spatial_hashing.jl b/test/cell_lists/spatial_hashing.jl index 4127a882..90a5c2dd 100644 --- a/test/cell_lists/spatial_hashing.jl +++ b/test/cell_lists/spatial_hashing.jl @@ -7,7 +7,6 @@ (8, 10, 6), (39, 40, 41) ] - name(size) = "$(length(size))D with $(prod(size)) Particles" @testset verbose=true "$(name(cloud_size))" for cloud_size in cloud_sizes coords = point_cloud(cloud_size, seed = 1) @@ -39,13 +38,7 @@ @test any(nhs.cell_list.collisions) end - # @testset verbose=true "Empty Cell and Collision Detection" begin - # # 1. Find empty lists - # # 2. Check for collisions with non-empty lists - # end - neighbors = [Int[] for _ in axes(coords, 2)] - foreach_point_neighbor(coords, coords, nhs, parallel = false) do point, neighbor, pos_diff, distance diff --git a/test/neighborhood_search.jl b/test/neighborhood_search.jl index 2a44f0f8..b2f39b7c 100644 --- a/test/neighborhood_search.jl +++ b/test/neighborhood_search.jl @@ -1,5 +1,3 @@ -using Infiltrator - # This file contains tests for the generic functions in `src/neighborhood_search.jl` and # tests comparing all NHS implementations against the `TrivialNeighborhoodSearch`. @testset verbose=true "All Neighborhood Searches" begin From 0b2d723b6929b440c1a587c0f96c69de9f11da1e Mon Sep 17 00:00:00 2001 From: RubberLanding Date: Tue, 29 Apr 2025 15:28:32 +0800 Subject: [PATCH 20/27] Resolve requested changes: - restore benchmark plotting code - remove NDIMS as a field for SpatialHashingCellList - change evaluation order of cell_collision in foreach_neighbor() to increase efficiency - add small, explicit test for cell collision - clean up --- benchmarks/plot.jl | 38 +++--------- src/cell_lists/spatial_hashing.jl | 28 +++++---- src/nhs_grid.jl | 4 +- test/cell_lists/spatial_hashing.jl | 98 +++++++++++++----------------- test/nhs_spatial_hash.jl | 70 --------------------- 5 files changed, 69 insertions(+), 169 deletions(-) delete mode 100644 test/nhs_spatial_hash.jl diff --git a/benchmarks/plot.jl b/benchmarks/plot.jl index 1a1bd83f..976cf045 100644 --- a/benchmarks/plot.jl +++ b/benchmarks/plot.jl @@ -1,6 +1,5 @@ using Plots using BenchmarkTools -# using Morton # Generate a rectangular point cloud include("../test/point_cloud.jl") @@ -36,16 +35,12 @@ include("benchmarks/benchmarks.jl") plot_benchmarks(benchmark_count_neighbors, (10, 10), 3) """ - function plot_benchmarks(benchmark, n_points_per_dimension, iterations; parallel = true, title = "", seed = 1, perturbation_factor_position = 1.0) - neighborhood_searches_names = ["GNHS with `DictionaryCellList`";; - "GNHS with `FullGridCellList`";; - "GNHS with `SpatialHashingCellList` n_points";; - "GNHS with `SpatialHashingCellList` 2 * n_points";; - "GNHS with `SpatialHashingCellList` 4 * n_points";; - ] + neighborhood_searches_names = ["TrivialNeighborhoodSearch";; + "GridNeighborhoodSearch";; + "PrecomputedNeighborhoodSearch"] # Multiply number of points in each iteration (roughly) by this factor scaling_factor = 4 @@ -64,26 +59,10 @@ function plot_benchmarks(benchmark, n_points_per_dimension, iterations; NDIMS = size(coordinates, 1) n_particles = size(coordinates, 2) - min_corner = Float32.(minimum(coordinates, dims = 2) .- search_radius) - max_corner = Float32.(maximum(coordinates, dims = 2) .+ search_radius) - neighborhood_searches = [ - GridNeighborhoodSearch{NDIMS}(; search_radius, n_points = n_particles, - update_strategy = nothing), - GridNeighborhoodSearch{NDIMS}(; - cell_list = FullGridCellList(; search_radius, - min_corner, - max_corner), - search_radius, n_points = n_particles), - GridNeighborhoodSearch{NDIMS}(; search_radius, n_points = n_particles, - cell_list = SpatialHashingCellList{NDIMS}(1 * - n_particles)), - GridNeighborhoodSearch{NDIMS}(; search_radius, n_points = n_particles, - cell_list = SpatialHashingCellList{NDIMS}(2 * - n_particles)), - GridNeighborhoodSearch{NDIMS}(; search_radius, n_points = n_particles, - cell_list = SpatialHashingCellList{NDIMS}(4 * - n_particles)) + TrivialNeighborhoodSearch{NDIMS}(; search_radius, eachpoint = 1:n_particles), + GridNeighborhoodSearch{NDIMS}(; search_radius, n_points = n_particles), + PrecomputedNeighborhoodSearch{NDIMS}(; search_radius, n_points = n_particles) ] for i in eachindex(neighborhood_searches) @@ -98,13 +77,10 @@ function plot_benchmarks(benchmark, n_points_per_dimension, iterations; end end - p = plot(n_particles_vec, times, + plot(n_particles_vec, times, xaxis = :log, yaxis = :log, xticks = (n_particles_vec, n_particles_vec), xlabel = "#particles", ylabel = "Runtime [s]", legend = :outerright, size = (750, 400), dpi = 200, label = neighborhood_searches_names, title = title) - display(p) - - return n_particles_vec, times, neighborhood_searches_names end diff --git a/src/cell_lists/spatial_hashing.jl b/src/cell_lists/spatial_hashing.jl index 4d83d46c..b9bec546 100644 --- a/src/cell_lists/spatial_hashing.jl +++ b/src/cell_lists/spatial_hashing.jl @@ -12,16 +12,17 @@ and the likelihood of hash collisions. - `list_size::Int`: Size of the hash map (e.g., `2 * n_points`) . """ -struct SpatialHashingCellList{CL, CI, CF} <: AbstractCellList +struct SpatialHashingCellList{NDIMS, CL, CI, CF} <: AbstractCellList points::CL coords::CI collisions::CF list_size::Int - NDIMS::Int end @inline index_type(::SpatialHashingCellList) = Int32 +@inline Base.ndims(::SpatialHashingCellList{NDIMS}) where {NDIMS} = NDIMS + function supported_update_strategies(::SpatialHashingCellList) return (SerialUpdate,) end @@ -30,11 +31,13 @@ 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(points, coords, collisions, list_size, NDIMS) + return SpatialHashingCellList{NDIMS, typeof(points), typeof(coords), + typeof(collisions)}(points, coords, collisions, list_size) end function Base.empty!(cell_list::SpatialHashingCellList) - (; list_size, NDIMS) = cell_list + (; 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] @@ -45,14 +48,17 @@ 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, NDIMS) = cell_list - + (; points, coords, collisions, list_size) = cell_list + NDIMS = ndims(cell_list) hash_key = spatial_hash(cell, list_size) - cell_coord = coords[hash_key] 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 @@ -65,7 +71,9 @@ end function copy_cell_list(cell_list::SpatialHashingCellList, search_radius, periodic_box) - (; NDIMS, list_size) = cell_list + (; list_size) = cell_list + NDIMS = ndims(cell_list) + return SpatialHashingCellList{NDIMS}(list_size) end @@ -83,10 +91,6 @@ end return coords == cell_index end -function spatial_hash(cell::CartesianIndex, list_size) - return spatial_hash(Tuple(cell), list_size) -end - function spatial_hash(cell::NTuple{1, Real}, list_size) return mod(cell[1] * 73856093, list_size) + 1 end diff --git a/src/nhs_grid.jl b/src/nhs_grid.jl index 2393b0cc..3ec16fca 100644 --- a/src/nhs_grid.jl +++ b/src/nhs_grid.jl @@ -458,6 +458,8 @@ function check_cell_collision(neighbor_cell_::CartesianIndex, 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 @@ -501,7 +503,7 @@ end # Check if `neighbor_coords` are in the cell `neighbor_cell_`. # For the `SpatialHashingCellList`, this might not be the case # if we have a collision. - if check_collision(neighbor_cell_, neighbor_coords, cell_list, neighborhood_search) && cell_collision + if cell_collision && check_collision(neighbor_cell_, neighbor_coords, cell_list, neighborhood_search) continue end diff --git a/test/cell_lists/spatial_hashing.jl b/test/cell_lists/spatial_hashing.jl index 90a5c2dd..67a8fc18 100644 --- a/test/cell_lists/spatial_hashing.jl +++ b/test/cell_lists/spatial_hashing.jl @@ -1,75 +1,63 @@ @testset verbose=true "SpatialHashingCellList" begin - # General test for 2D and 3D - @testset verbose=true "Compare Against `TrivialNeighborhoodSearch`" begin - cloud_sizes = [ - (10, 11), - (100, 90), - (8, 10, 6), - (39, 40, 41) - ] - name(size) = "$(length(size))D with $(prod(size)) Particles" - @testset verbose=true "$(name(cloud_size))" for cloud_size in cloud_sizes - coords = point_cloud(cloud_size, seed = 1) - NDIMS = length(cloud_size) - n_points = size(coords, 2) - search_radius = 2.5 - - # Compute expected neighbor lists by brute-force looping over all points - # as potential neighbors (`TrivialNeighborhoodSearch`). - trivial_nhs = TrivialNeighborhoodSearch{NDIMS}(; search_radius, - eachpoint = axes(coords, 2)) - - neighbors_expected = [Int[] for _ in axes(coords, 2)] - - foreach_point_neighbor(coords, coords, trivial_nhs, - parallel = false) do point, neighbor, - pos_diff, distance - append!(neighbors_expected[point], neighbor) - end - - nhs = GridNeighborhoodSearch{NDIMS}(; search_radius, n_points, - cell_list = SpatialHashingCellList{NDIMS}(2 * - n_points)) + @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 - initialize!(nhs, coords, coords) + nhs = GridNeighborhoodSearch{2}(; search_radius, n_points, + cell_list = SpatialHashingCellList{NDIMS}(n_points)) + initialize_grid!(nhs, coordinates) - # Test if there are any collisions. - @testset verbose=true "Collision Detection" begin - @test any(nhs.cell_list.collisions) - end + @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) - neighbors = [Int[] for _ in axes(coords, 2)] - foreach_point_neighbor(coords, coords, nhs, - parallel = false) do point, neighbor, - pos_diff, distance - push!(neighbors[point], neighbor) - end + @test cell1_hash == cell2_hash + end - @test sort.(neighbors) == neighbors_expected + neighbors = Int[] + foreach_neighbor(coordinates, coordinates, nhs, + point_index) do point, neighbor, pos_diff, distance + push!(neighbors, neighbor) end + + @test neighbors == [1] end - # Test list behavior with empty cells - @testset "Collision Handling with empty cells" begin - # The point is in cell (-1, 0) which has a hash collision with cell (-2, -1) - point = [-0.05, 0.05] - coordinates = hcat(point) + @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)) + cell_list = SpatialHashingCellList{NDIMS}(n_points)) initialize_grid!(nhs, coordinates) - found_neighbors = Int[] - foreach_neighbor(coordinates, coordinates, nhs, - point_index) do point, neighbor, pos_diff, distance - push!(found_neighbors, neighbor) + @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) + + @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 - correct_neighbors = [1] - @test correct_neighbors == found_neighbors + @test neighbors[1] == [1] + @test neighbors[2] == [2] end end diff --git a/test/nhs_spatial_hash.jl b/test/nhs_spatial_hash.jl deleted file mode 100644 index eb948440..00000000 --- a/test/nhs_spatial_hash.jl +++ /dev/null @@ -1,70 +0,0 @@ -@testset verbose=true "SpatialHashNeighborhoodSearch" begin - @testset "Rectangular Point Cloud 2D" begin - #### Setup - # Rectangle of equidistantly spaced points - # from (x, y) = (-0.25, -0.25) to (x, y) = (0.35, 0.35). - range = -0.25:0.1:0.35 - coordinates1 = hcat(collect.(Iterators.product(range, range))...) - n_points = size(coordinates1, 2) - - point_position1 = [0.05, 0.05] - search_radius = 0.1 - - # Create neighborhood search - nhs1 = GridNeighborhoodSearch{2}(; search_radius, n_points, - cell_list = SpatialHashingCellList{2}(2 * n_points)) - initialize_grid!(nhs1, coordinates1) - - # Get each neighbor for `point_position1` - neighbors1 = sort(collect(PointNeighbors.eachneighbor(point_position1, nhs1))) - - # Move points and update NHS - coordinates2 = coordinates1 .+ [1.4, -3.5] - empty!(nhs1.cell_list) - initialize_grid!(nhs1, coordinates2) - - # Get each neighbor for updated NHS - neighbors2 = sort(collect(PointNeighbors.eachneighbor(point_position1, nhs1))) - - # Change position - point_position2 = point_position1 .+ [1.4, -3.5] - - # Get each neighbor for `point_position2` - neighbors3 = sort(collect(PointNeighbors.eachneighbor(point_position2, nhs1))) - - # Double search radius - nhs2 = GridNeighborhoodSearch{2}(search_radius = 2 * search_radius, - n_points = size(coordinates1, 2), - cell_list = SpatialHashingCellList{2}(2 * n_points)) - initialize!(nhs2, coordinates1, coordinates1) - - # Get each neighbor in double search radius - neighbors4 = sort(collect(PointNeighbors.eachneighbor(point_position1, nhs2))) - - # Move points - coordinates2 = coordinates1 .+ [0.4, -0.4] - empty!(nhs2.cell_list) - initialize_grid!(nhs2, coordinates2) - - # Get each neighbor in double search radius - neighbors5 = sort(collect(PointNeighbors.eachneighbor(point_position1, nhs2))) - - #### Verification against lists of potential neighbors built by hand - correct_neighbors1 = [17, 18, 19, 24, 25, 26, 31, 32, 33] - @test all(x -> x in neighbors1, correct_neighbors1) - - correct_neighbors2 = Int[] - @test all(x -> x in neighbors2, correct_neighbors2) - - correct_neighbors3 = [17, 18, 19, 24, 25, 26, 31, 32, 33] - @test all(x -> x in neighbors3, correct_neighbors3) - - correct_neighbors4 = [ - 9, 10, 11, 12, 13, 14, 16, 17, 18, 19, 20, 21, 23, 24, 25, 26, 27, 28, 30, 31, - 32, 33, 34, 35, 37, 38, 39, 40, 41, 42, 44, 45, 46, 47, 48, 49] - @test all(x -> x in neighbors4, correct_neighbors4) - - correct_neighbors5 = [36, 37, 38, 43, 44, 45] - @test all(x -> x in neighbors5, correct_neighbors5) - end -end; From 9d15adb6c055f5e7ab48c6408a8acb03402fac25 Mon Sep 17 00:00:00 2001 From: RubberLanding Date: Fri, 2 May 2025 18:37:50 +0800 Subject: [PATCH 21/27] Resolve requested changes: - add comments - adjust version in Project.toml --- src/cell_lists/spatial_hashing.jl | 25 +++++++++---------------- src/nhs_grid.jl | 21 +++++++++++++++------ test/Project.toml | 31 ++++++++++++++++++++++--------- test/neighborhood_search.jl | 12 +++++++----- 4 files changed, 53 insertions(+), 36 deletions(-) diff --git a/src/cell_lists/spatial_hashing.jl b/src/cell_lists/spatial_hashing.jl index b9bec546..60d2c3db 100644 --- a/src/cell_lists/spatial_hashing.jl +++ b/src/cell_lists/spatial_hashing.jl @@ -4,8 +4,8 @@ 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. (2001)](@cite Ihmsen2003). By using a hash map, which only stores non-empty cells, -the domain is effectively infinite. The size of the hash map is recommended to be approximately twice the number of particles balance memory consumption -and the likelihood of hash collisions. +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`). @@ -13,10 +13,10 @@ and the likelihood of hash collisions. """ struct SpatialHashingCellList{NDIMS, CL, CI, CF} <: AbstractCellList - points::CL - coords::CI - collisions::CF - list_size::Int + points :: CL + coords :: CI + collisions :: CF + list_size :: Int end @inline index_type(::SpatialHashingCellList) = Int32 @@ -78,8 +78,7 @@ function copy_cell_list(cell_list::SpatialHashingCellList, search_radius, end @inline function Base.getindex(cell_list::SpatialHashingCellList, cell::Tuple) - (; points) = cell_list - return points[spatial_hash(cell, length(points))] + return cell_list.points[spatial_hash(cell, length(points))] end @inline function Base.getindex(cell_list::SpatialHashingCellList, i::Integer) @@ -91,6 +90,7 @@ end 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 @@ -105,11 +105,4 @@ 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 - -function spatial_hash(cell::NTuple{3, Real}, index, list_size) - i, j, k = cell - - return mod(xor(i * 73856093, j * 19349663, k * 83492791, index * 7238423947), - list_size) + 1 -end +end \ No newline at end of file diff --git a/src/nhs_grid.jl b/src/nhs_grid.jl index 3ec16fca..a9504e42 100644 --- a/src/nhs_grid.jl +++ b/src/nhs_grid.jl @@ -442,9 +442,12 @@ function update_grid!(neighborhood_search::Union{GridNeighborhoodSearch{<:Any, 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 @@ -455,6 +458,7 @@ end function check_cell_collision(neighbor_cell_::CartesianIndex, cell_list, nhs) + # This is only relevant for the `SpatialHashingCellList` return false end @@ -466,11 +470,15 @@ function check_cell_collision(neighbor_cell_::CartesianIndex, 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`. +# than looping over `eachneighbor`. @inline function foreach_neighbor(f, neighbor_system_coords, neighborhood_search::GridNeighborhoodSearch, point, point_coords, search_radius) @@ -480,9 +488,11 @@ end 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_] @@ -500,15 +510,14 @@ end if distance2 <= search_radius^2 distance = sqrt(distance2) - # Check if `neighbor_coords` are in the cell `neighbor_cell_`. - # For the `SpatialHashingCellList`, this might not be the case - # if we have a collision. + # 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`. + # compared to not using `foreach_point_neighbor @inline f(point, neighbor, pos_diff, distance) end end diff --git a/test/Project.toml b/test/Project.toml index c826e5bc..9250ef7d 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,12 +1,25 @@ +name = "PointNeighbors" +uuid = "1c4d5385-0a27-49de-8e2c-43b175c8985c" +authors = ["Erik Faulhaber "] +version = "0.4.12-pre" + [deps] -BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" -Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" -Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -TrixiParticles = "66699cd8-9c01-4e9d-a059-b96c86d16b3a" +Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" +Atomix = "a9b6321e-bd34-4604-b9c9-b65b8de01458" +GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527" +KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588" +Reexport = "189a3867-3050-52da-a836-e630ba90ab69" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [compat] -BenchmarkTools = "1" -Plots = "1" -Test = "1" -TrixiParticles = "0.2" +Adapt = "4" +Atomix = "1" +GPUArraysCore = "0.2" +KernelAbstractions = "0.9" +LinearAlgebra = "1" +Polyester = "0.7.5" +Reexport = "1" +StaticArrays = "1" +julia = "1.10" \ No newline at end of file diff --git a/test/neighborhood_search.jl b/test/neighborhood_search.jl index b2f39b7c..97e36f4e 100644 --- a/test/neighborhood_search.jl +++ b/test/neighborhood_search.jl @@ -87,7 +87,8 @@ backend = Vector{Vector{Int32}})), PrecomputedNeighborhoodSearch{NDIMS}(periodic_box = periodic_boxes[i]), GridNeighborhoodSearch{NDIMS}(periodic_box = periodic_boxes[i], - cell_list = SpatialHashingCellList{NDIMS}(2 * n_points)) + cell_list = SpatialHashingCellList{NDIMS}(2 * + n_points)) ] copied_nhs = copy_neighborhood_search.(template_nhs, search_radius, n_points) append!(neighborhood_searches, copied_nhs) @@ -192,7 +193,7 @@ PrecomputedNeighborhoodSearch{NDIMS}(; search_radius, n_points), GridNeighborhoodSearch{NDIMS}(; search_radius, n_points, cell_list = SpatialHashingCellList{NDIMS}(2 * - n_points)) + n_points)) ] names = [ @@ -204,8 +205,8 @@ "`GridNeighborhoodSearch` with `FullGridCellList` with `DynamicVectorOfVectors` and `SemiParallelUpdate`", "`GridNeighborhoodSearch` with `FullGridCellList` with `Vector{Vector}`", "`PrecomputedNeighborhoodSearch`", - "`GridNeighborhoodSearch` with `SpatialHashingCellList`", - ] + "`GridNeighborhoodSearch` with `SpatialHashingCellList`" + ] # Also test copied templates template_nhs = [ @@ -224,7 +225,8 @@ max_corner, backend = Vector{Vector{Int32}})), PrecomputedNeighborhoodSearch{NDIMS}(), - GridNeighborhoodSearch{NDIMS}(cell_list = SpatialHashingCellList{NDIMS}(2 * n_points)) + 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) From 3a481f5e189026f43e42fcace895e2f9e2b283ec Mon Sep 17 00:00:00 2001 From: RubberLanding Date: Fri, 2 May 2025 19:46:18 +0800 Subject: [PATCH 22/27] Minor changes. --- Project.toml | 4 ++-- src/cell_lists/spatial_hashing.jl | 4 ++-- src/nhs_grid.jl | 12 +++++++----- test/Project.toml | 25 ------------------------- 4 files changed, 11 insertions(+), 34 deletions(-) delete mode 100644 test/Project.toml diff --git a/Project.toml b/Project.toml index 92fb96c6..53010f6c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "PointNeighbors" uuid = "1c4d5385-0a27-49de-8e2c-43b175c8985c" authors = ["Erik Faulhaber "] -version = "0.5" +version = "0.4.9" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" @@ -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/cell_lists/spatial_hashing.jl b/src/cell_lists/spatial_hashing.jl index 60d2c3db..abfe086e 100644 --- a/src/cell_lists/spatial_hashing.jl +++ b/src/cell_lists/spatial_hashing.jl @@ -78,7 +78,7 @@ function copy_cell_list(cell_list::SpatialHashingCellList, search_radius, end @inline function Base.getindex(cell_list::SpatialHashingCellList, cell::Tuple) - return cell_list.points[spatial_hash(cell, length(points))] + return cell_list.points[spatial_hash(cell, length(cell_list.points))] end @inline function Base.getindex(cell_list::SpatialHashingCellList, i::Integer) @@ -105,4 +105,4 @@ 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 \ No newline at end of file +end diff --git a/src/nhs_grid.jl b/src/nhs_grid.jl index a9504e42..c8d9ed1f 100644 --- a/src/nhs_grid.jl +++ b/src/nhs_grid.jl @@ -480,8 +480,8 @@ 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) + neighborhood_search::GridNeighborhoodSearch, + point, point_coords, search_radius) (; cell_list, periodic_box) = neighborhood_search cell = cell_coords(point_coords, neighborhood_search) @@ -491,8 +491,8 @@ end # Boolean to indicate if this cell has a collision (only with `SpatialHashingCellList`) cell_collision = check_cell_collision(neighbor_cell_, - cell_list, neighborhood_search) - + cell_list, neighborhood_search) + for neighbor_ in eachindex(neighbors) neighbor = @inbounds neighbors[neighbor_] @@ -512,7 +512,9 @@ end # 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) + if cell_collision && + check_collision(neighbor_cell_, neighbor_coords, cell_list, + neighborhood_search) continue end diff --git a/test/Project.toml b/test/Project.toml deleted file mode 100644 index 9250ef7d..00000000 --- a/test/Project.toml +++ /dev/null @@ -1,25 +0,0 @@ -name = "PointNeighbors" -uuid = "1c4d5385-0a27-49de-8e2c-43b175c8985c" -authors = ["Erik Faulhaber "] -version = "0.4.12-pre" - -[deps] -Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" -Atomix = "a9b6321e-bd34-4604-b9c9-b65b8de01458" -GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527" -KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c" -LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588" -Reexport = "189a3867-3050-52da-a836-e630ba90ab69" -StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" - -[compat] -Adapt = "4" -Atomix = "1" -GPUArraysCore = "0.2" -KernelAbstractions = "0.9" -LinearAlgebra = "1" -Polyester = "0.7.5" -Reexport = "1" -StaticArrays = "1" -julia = "1.10" \ No newline at end of file From 33f44cf1337f5b96f834b3f2b718e39627b357cb Mon Sep 17 00:00:00 2001 From: RubberLanding Date: Tue, 6 May 2025 10:36:36 +0800 Subject: [PATCH 23/27] Resolve requested changes. --- Project.toml | 2 +- src/nhs_grid.jl | 2 +- test/Project.toml | 12 ++++++++++++ test/cell_lists/spatial_hashing.jl | 8 +++++++- 4 files changed, 21 insertions(+), 3 deletions(-) create mode 100644 test/Project.toml diff --git a/Project.toml b/Project.toml index 53010f6c..3fd7950c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "PointNeighbors" uuid = "1c4d5385-0a27-49de-8e2c-43b175c8985c" authors = ["Erik Faulhaber "] -version = "0.4.9" +version = "0.4.9-dev" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" diff --git a/src/nhs_grid.jl b/src/nhs_grid.jl index c8d9ed1f..55fbfbe3 100644 --- a/src/nhs_grid.jl +++ b/src/nhs_grid.jl @@ -519,7 +519,7 @@ end end # Inline to avoid loss of performance - # compared to not using `foreach_point_neighbor + # compared to not using `foreach_point_neighbor`. @inline f(point, neighbor, pos_diff, distance) end end diff --git a/test/Project.toml b/test/Project.toml new file mode 100644 index 00000000..c826e5bc --- /dev/null +++ b/test/Project.toml @@ -0,0 +1,12 @@ +[deps] +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +TrixiParticles = "66699cd8-9c01-4e9d-a059-b96c86d16b3a" + +[compat] +BenchmarkTools = "1" +Plots = "1" +Test = "1" +TrixiParticles = "0.2" diff --git a/test/cell_lists/spatial_hashing.jl b/test/cell_lists/spatial_hashing.jl index 67a8fc18..16bdb499 100644 --- a/test/cell_lists/spatial_hashing.jl +++ b/test/cell_lists/spatial_hashing.jl @@ -16,7 +16,10 @@ 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 @@ -46,7 +49,10 @@ 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 From 44fab7f080606048f08f85f923ed3b41243bcb4a Mon Sep 17 00:00:00 2001 From: RubberLanding Date: Tue, 13 May 2025 18:55:05 +0800 Subject: [PATCH 24/27] Resolve requested change for docs for SpatialHashingCellList. --- src/cell_lists/spatial_hashing.jl | 6 +++--- src/nhs_grid.jl | 11 ++++------- test/cell_lists/spatial_hashing.jl | 2 +- test/neighborhood_search.jl | 5 +++-- 4 files changed, 11 insertions(+), 13 deletions(-) diff --git a/src/cell_lists/spatial_hashing.jl b/src/cell_lists/spatial_hashing.jl index abfe086e..6fffd28a 100644 --- a/src/cell_lists/spatial_hashing.jl +++ b/src/cell_lists/spatial_hashing.jl @@ -3,9 +3,9 @@ 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. (2001)](@cite Ihmsen2003). By using a hash map, which only stores 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. +[as described by Ihmsen et al. (2001)](@cite Ihmsen2003). 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`). diff --git a/src/nhs_grid.jl b/src/nhs_grid.jl index 2f1c10c3..e480046a 100644 --- a/src/nhs_grid.jl +++ b/src/nhs_grid.jl @@ -312,8 +312,7 @@ end # `each_cell_index(cell_list)` might return a `KeySet`, which has to be `collect`ed # first to be able to be used in a threaded loop. This function takes care of that. - @threaded parallelization_backend for cell_index in - each_cell_index_threadable(cell_list) + @threaded parallelization_backend for cell_index in each_cell_index_threadable(cell_list) mark_changed_cell!(neighborhood_search, cell_index, y) end end @@ -361,13 +360,12 @@ function update_grid!(neighborhood_search::GridNeighborhoodSearch{<:Any, # We can work around this by using the old lengths. # TODO this is hardcoded for the `FullGridCellList` @threaded parallelization_backend for i in eachindex(update_buffer, - cell_list.cells.lengths) + cell_list.cells.lengths) update_buffer[i] = cell_list.cells.lengths[i] end # Add points to new cells - @threaded parallelization_backend for cell_index in - each_cell_index_threadable(cell_list) + @threaded parallelization_backend for cell_index in each_cell_index_threadable(cell_list) for i in 1:update_buffer[cell_index] point = cell_list.cells.backend[i, cell_index] point_coords = extract_svector(y, Val(ndims(neighborhood_search)), point) @@ -381,8 +379,7 @@ function update_grid!(neighborhood_search::GridNeighborhoodSearch{<:Any, end # Remove points from old cells - @threaded parallelization_backend for cell_index in - each_cell_index_threadable(cell_list) + @threaded parallelization_backend for cell_index in each_cell_index_threadable(cell_list) points = cell_list[cell_index] # WARNING!!! diff --git a/test/cell_lists/spatial_hashing.jl b/test/cell_lists/spatial_hashing.jl index 16bdb499..2c83f611 100644 --- a/test/cell_lists/spatial_hashing.jl +++ b/test/cell_lists/spatial_hashing.jl @@ -51,7 +51,7 @@ 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 diff --git a/test/neighborhood_search.jl b/test/neighborhood_search.jl index a1ef9b3b..16d66dcd 100644 --- a/test/neighborhood_search.jl +++ b/test/neighborhood_search.jl @@ -132,10 +132,11 @@ seeds = [1, 2] name(size, - seed) = "$(length(size))D with $(prod(size)) Particles " * - "($(seed == 1 ? "`initialize!`" : "`update!`"))" + seed) = "$(length(size))D with $(prod(size)) Particles " * + "($(seed == 1 ? "`initialize!`" : "`update!`"))" @testset verbose=true "$(name(cloud_size, seed)))" for cloud_size in cloud_sizes, seed in seeds + coords = point_cloud(cloud_size, seed = seed) NDIMS = length(cloud_size) n_points = size(coords, 2) From a3f622db11ad7d32a042802f7e0d102f4799d816 Mon Sep 17 00:00:00 2001 From: RubberLanding Date: Tue, 13 May 2025 19:09:49 +0800 Subject: [PATCH 25/27] Change citation of Ihmsen in doc for SpatialHashingCellList. --- src/cell_lists/spatial_hashing.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/cell_lists/spatial_hashing.jl b/src/cell_lists/spatial_hashing.jl index 6fffd28a..4822cb2c 100644 --- a/src/cell_lists/spatial_hashing.jl +++ b/src/cell_lists/spatial_hashing.jl @@ -2,9 +2,9 @@ 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. (2001)](@cite Ihmsen2003). 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 +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 From 8036f471be6bf95f36e8b9206b99a917d5db1806 Mon Sep 17 00:00:00 2001 From: RubberLanding Date: Thu, 15 May 2025 15:44:33 +0800 Subject: [PATCH 26/27] Fix formatting. --- test/nhs_grid.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/test/nhs_grid.jl b/test/nhs_grid.jl index 59facc83..c1414497 100644 --- a/test/nhs_grid.jl +++ b/test/nhs_grid.jl @@ -190,8 +190,7 @@ 173, 178, 179, 180, 213, 214, 215, 220, 221, 222, 227, 228, 229] update_strategies = (SerialUpdate(), ParallelUpdate()) - @testset verbose=true "eachindex_y $update_strategy" for update_strategy in - update_strategies + @testset verbose=true "eachindex_y $update_strategy" for update_strategy in update_strategies # Test that `eachindex_y` is passed correctly to the neighborhood search. # This requires `SerialUpdate` or `ParallelUpdate`. min_corner = min.(minimum(coordinates1, dims = 2), From 7ce370c20c2d5e5cf665c0205d3f8dff87f30a1a Mon Sep 17 00:00:00 2001 From: RubberLanding Date: Thu, 15 May 2025 16:09:58 +0800 Subject: [PATCH 27/27] Fix formatting. --- src/nhs_grid.jl | 8 +++++--- test/neighborhood_search.jl | 5 ++--- test/nhs_grid.jl | 3 ++- 3 files changed, 9 insertions(+), 7 deletions(-) diff --git a/src/nhs_grid.jl b/src/nhs_grid.jl index 3e38bdd4..ba613ecd 100644 --- a/src/nhs_grid.jl +++ b/src/nhs_grid.jl @@ -397,12 +397,13 @@ function update_grid!(neighborhood_search::GridNeighborhoodSearch{<:Any, # We can work around this by using the old lengths. # TODO this is hardcoded for the `FullGridCellList` @threaded parallelization_backend for i in eachindex(update_buffer, - cell_list.cells.lengths) + cell_list.cells.lengths) update_buffer[i] = cell_list.cells.lengths[i] end # Add points to new cells - @threaded parallelization_backend for cell_index in each_cell_index_threadable(cell_list) + @threaded parallelization_backend for cell_index in + each_cell_index_threadable(cell_list) for i in 1:update_buffer[cell_index] point = cell_list.cells.backend[i, cell_index] point_coords = extract_svector(y, Val(ndims(neighborhood_search)), point) @@ -416,7 +417,8 @@ function update_grid!(neighborhood_search::GridNeighborhoodSearch{<:Any, end # Remove points from old cells - @threaded parallelization_backend for cell_index in each_cell_index_threadable(cell_list) + @threaded parallelization_backend for cell_index in + each_cell_index_threadable(cell_list) points = cell_list[cell_index] # WARNING!!! diff --git a/test/neighborhood_search.jl b/test/neighborhood_search.jl index 16d66dcd..a1ef9b3b 100644 --- a/test/neighborhood_search.jl +++ b/test/neighborhood_search.jl @@ -132,11 +132,10 @@ seeds = [1, 2] name(size, - seed) = "$(length(size))D with $(prod(size)) Particles " * - "($(seed == 1 ? "`initialize!`" : "`update!`"))" + seed) = "$(length(size))D with $(prod(size)) Particles " * + "($(seed == 1 ? "`initialize!`" : "`update!`"))" @testset verbose=true "$(name(cloud_size, seed)))" for cloud_size in cloud_sizes, seed in seeds - coords = point_cloud(cloud_size, seed = seed) NDIMS = length(cloud_size) n_points = size(coords, 2) diff --git a/test/nhs_grid.jl b/test/nhs_grid.jl index c1414497..59facc83 100644 --- a/test/nhs_grid.jl +++ b/test/nhs_grid.jl @@ -190,7 +190,8 @@ 173, 178, 179, 180, 213, 214, 215, 220, 221, 222, 227, 228, 229] update_strategies = (SerialUpdate(), ParallelUpdate()) - @testset verbose=true "eachindex_y $update_strategy" for update_strategy in update_strategies + @testset verbose=true "eachindex_y $update_strategy" for update_strategy in + update_strategies # Test that `eachindex_y` is passed correctly to the neighborhood search. # This requires `SerialUpdate` or `ParallelUpdate`. min_corner = min.(minimum(coordinates1, dims = 2),