From 5d367b3214adc04bce458f050ed43f0a533ff028 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Tue, 17 Mar 2026 12:30:25 +0100 Subject: [PATCH 1/4] Fix bug due to rounding errors in periodicity computation --- src/cell_lists/full_grid.jl | 29 ------------------------- src/nhs_grid.jl | 42 ++++++++++++++++++------------------- 2 files changed, 21 insertions(+), 50 deletions(-) diff --git a/src/cell_lists/full_grid.jl b/src/cell_lists/full_grid.jl index 1edd165e..1fd4d0a6 100644 --- a/src/cell_lists/full_grid.jl +++ b/src/cell_lists/full_grid.jl @@ -81,35 +81,6 @@ function FullGridCellList(; min_corner, max_corner, return FullGridCellList(cells, linear_indices, min_corner, max_corner) end -@inline function cell_coords(coords, periodic_box::Nothing, cell_list::FullGridCellList, - cell_size) - (; min_corner) = cell_list - - # Subtract `min_corner` to offset coordinates so that the min corner of the grid - # corresponds to the (1, 1, 1) cell. - return Tuple(floor_to_int.((coords .- min_corner) ./ cell_size)) .+ 1 -end - -@inline function cell_coords(coords, periodic_box::PeriodicBox, cell_list::FullGridCellList, - cell_size) - # Subtract `periodic_box.min_corner` to offset coordinates so that the min corner - # of the grid corresponds to the (0, 0, 0) cell. - offset_coords = periodic_coords(coords, periodic_box) .- periodic_box.min_corner - - # Add 2, so that the min corner will be the (2, 2, 2)-cell. - # With this, we still have one padding layer in each direction around the periodic box, - # just like without using a periodic box. - # This is not needed for finding neighbor cells, but to make the bounds check - # work the same way as without a periodic box. - return Tuple(floor_to_int.(offset_coords ./ cell_size)) .+ 2 -end - -@inline function periodic_cell_index(cell_index, ::PeriodicBox, n_cells, - cell_list::FullGridCellList) - # 2-based modulo to match the indexing of the periodic box explained above. - return mod.(cell_index .- 2, n_cells) .+ 2 -end - function Base.empty!(cell_list::FullGridCellList) (; cells) = cell_list diff --git a/src/nhs_grid.jl b/src/nhs_grid.jl index 21759ab2..6b95b109 100644 --- a/src/nhs_grid.jl +++ b/src/nhs_grid.jl @@ -579,37 +579,37 @@ end end @inline function periodic_cell_index(cell_index, neighborhood_search) - (; n_cells, periodic_box, cell_list) = neighborhood_search + (; n_cells, periodic_box) = neighborhood_search - periodic_cell_index(cell_index, periodic_box, n_cells, cell_list) + periodic_cell_index(cell_index, periodic_box, n_cells) end -@inline periodic_cell_index(cell_index, ::Nothing, n_cells, cell_list) = cell_index - -@inline function periodic_cell_index(cell_index, ::PeriodicBox, n_cells, cell_list) - # 1-based modulo - return mod1.(cell_index, n_cells) +@inline periodic_cell_index(cell_index, ::Nothing, n_cells) = cell_index + +@inline function periodic_cell_index(cell_index, ::PeriodicBox, n_cells) + # With the `DictionaryCellList`, the grid is conceptually infinite, + # so we can use any offset for the modulo operation. + # With the `FullGridCellList`, we need 2-based modulo, so that the min corner + # will be the (2, 2, 2)-cell. + # With this, we still have one padding layer in each direction around the periodic box, + # just like without using a periodic box. + # This is not needed for finding neighbor cells, but to make the bounds check + # work the same way as without a periodic box (cells in bounds are 2:end-1). + return mod.(cell_index .- 2, n_cells) .+ 2 end @inline function cell_coords(coords, neighborhood_search) - (; periodic_box, cell_list, cell_size) = neighborhood_search + # Compute cell coordinates ignoring periodicity at first + nonperiodic_cell_coords_ = nonperiodic_cell_coords(coords, neighborhood_search) - return cell_coords(coords, periodic_box, cell_list, cell_size) + # Now return (without periodicity) or convert to periodic cell index (with periodicity) + return periodic_cell_index(nonperiodic_cell_coords_, neighborhood_search) end -@inline function cell_coords(coords, periodic_box::Nothing, cell_list, cell_size) - return Tuple(floor_to_int.(coords ./ cell_size)) -end +@inline function nonperiodic_cell_coords(coords, neighborhood_search) + (; cell_size) = neighborhood_search -@inline function cell_coords(coords, periodic_box::PeriodicBox, cell_list, cell_size) - # Subtract `min_corner` to offset coordinates so that the min corner of the periodic - # box corresponds to the (0, 0, 0) cell of the NHS. - # This way, there are no partial cells in the domain if the domain size is an integer - # multiple of the cell size (which is required, see the constructor). - offset_coords = periodic_coords(coords, periodic_box) .- periodic_box.min_corner - - # Add one for 1-based indexing. The min corner will be the (1, 1, 1)-cell. - return Tuple(floor_to_int.(offset_coords ./ cell_size)) .+ 1 + return Tuple(floor_to_int.(coords ./ cell_size)) end function copy_neighborhood_search(nhs::GridNeighborhoodSearch, search_radius, n_points; From 52e799b588804589af8f21546599075489729af2 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Tue, 23 Jun 2026 12:58:56 +0200 Subject: [PATCH 2/4] Fix --- src/cell_lists/full_grid.jl | 7 ++++ src/nhs_grid.jl | 19 +++++++---- test/nhs_grid.jl | 67 +++++++++++++++++++++---------------- 3 files changed, 59 insertions(+), 34 deletions(-) diff --git a/src/cell_lists/full_grid.jl b/src/cell_lists/full_grid.jl index 1fd4d0a6..c1702df6 100644 --- a/src/cell_lists/full_grid.jl +++ b/src/cell_lists/full_grid.jl @@ -81,6 +81,13 @@ function FullGridCellList(; min_corner, max_corner, return FullGridCellList(cells, linear_indices, min_corner, max_corner) end +@inline function nonperiodic_cell_coords(coords, cell_list::FullGridCellList, cell_size) + # Subtract `min_corner` to offset coordinates so that the (padded) min corner + # of the grid corresponds to the (1, 1, 1) cell. + # The unpadded min corner (that is passed by the user) corresponds to the (2, 2, 2) cell. + return Tuple(floor_to_int.((coords .- cell_list.min_corner) ./ cell_size) .+ 1) +end + function Base.empty!(cell_list::FullGridCellList) (; cells) = cell_list diff --git a/src/nhs_grid.jl b/src/nhs_grid.jl index 6b95b109..c3365fff 100644 --- a/src/nhs_grid.jl +++ b/src/nhs_grid.jl @@ -587,14 +587,17 @@ end @inline periodic_cell_index(cell_index, ::Nothing, n_cells) = cell_index @inline function periodic_cell_index(cell_index, ::PeriodicBox, n_cells) - # With the `DictionaryCellList`, the grid is conceptually infinite, - # so we can use any offset for the modulo operation. - # With the `FullGridCellList`, we need 2-based modulo, so that the min corner - # will be the (2, 2, 2)-cell. + # With the `FullGridCellList`, nonperiodic cell coords are in the range 2:n_cells+1, + # so we need an offset of 2 to preserve this, so that the min corner will still + # be the (2, 2, 2)-cell. # With this, we still have one padding layer in each direction around the periodic box, # just like without using a periodic box. - # This is not needed for finding neighbor cells, but to make the bounds check + # This padding is not needed for finding neighbor cells, but to make the bounds check # work the same way as without a periodic box (cells in bounds are 2:end-1). + # + # With the `DictionaryCellList`, the grid is conceptually infinite, + # so we can use any offset for the modulo operation. + # The offset only determines to which part of R^n the periodic box is mapped. return mod.(cell_index .- 2, n_cells) .+ 2 end @@ -607,8 +610,12 @@ end end @inline function nonperiodic_cell_coords(coords, neighborhood_search) - (; cell_size) = neighborhood_search + (; cell_list, cell_size) = neighborhood_search + + return nonperiodic_cell_coords(coords, cell_list, cell_size) +end +@inline function nonperiodic_cell_coords(coords, cell_list, cell_size) return Tuple(floor_to_int.(coords ./ cell_size)) end diff --git a/test/nhs_grid.jl b/test/nhs_grid.jl index 17ff2fa9..3f4cbb29 100644 --- a/test/nhs_grid.jl +++ b/test/nhs_grid.jl @@ -62,12 +62,11 @@ coords2 = [NaN, 0] coords3 = [typemax(Int) + 1.0, -typemax(Int) - 1.0] - @test PointNeighbors.cell_coords(coords1, nothing, nothing, (1.0, 1.0)) == - (typemax(Int), typemin(Int)) - @test PointNeighbors.cell_coords(coords2, nothing, nothing, (1.0, 1.0)) == - (typemax(Int), 0) - @test PointNeighbors.cell_coords(coords3, nothing, nothing, (1.0, 1.0)) == - (typemax(Int), typemin(Int)) + nhs = GridNeighborhoodSearch{2}(search_radius = 1.0, n_points = 1) + + @test PointNeighbors.cell_coords(coords1, nhs) == (typemax(Int), typemin(Int)) + @test PointNeighbors.cell_coords(coords2, nhs) == (typemax(Int), 0) + @test PointNeighbors.cell_coords(coords3, nhs) == (typemax(Int), typemin(Int)) # The full grid cell list adds one to the coordinates to avoid zero-indexing. # This corner case is not relevant, as `typemax` coordinates will always be out of @@ -75,12 +74,12 @@ cell_list = FullGridCellList(min_corner = (0.0, 0.0), max_corner = (1.0, 1.0), search_radius = 1.0) - @test PointNeighbors.cell_coords(coords1, nothing, cell_list, (1.0, 1.0)) == - (typemax(Int), typemin(Int)) .+ 1 - @test PointNeighbors.cell_coords(coords2, nothing, cell_list, (1.0, 1.0)) == - (typemax(Int), 1) .+ 1 - @test PointNeighbors.cell_coords(coords3, nothing, cell_list, (1.0, 1.0)) == - (typemax(Int), typemin(Int)) .+ 1 + nhs = GridNeighborhoodSearch{2}(search_radius = 1.0, n_points = 1, + cell_list = cell_list) + + @test PointNeighbors.cell_coords(coords1, nhs) == (typemax(Int), typemin(Int)) .+ 1 + @test PointNeighbors.cell_coords(coords2, nhs) == (typemax(Int), 1) .+ 1 + @test PointNeighbors.cell_coords(coords3, nhs) == (typemax(Int), typemin(Int)) .+ 1 end @testset "Rectangular Point Cloud 2D" begin @@ -249,28 +248,40 @@ PeriodicBox(min_corner = [-0.1, -0.2], max_corner = [0.2, 0.4]), # The `GridNeighborhoodSearch` is forced to round up the cell sizes in this test # to avoid split cells. - PeriodicBox(min_corner = [-0.1, -0.2], max_corner = [0.205, 0.43]), + PeriodicBox(min_corner = [-0.1, -0.205], max_corner = [0.205, 0.435]), PeriodicBox(min_corner = [-0.1, -0.2, 0.05], max_corner = [0.2, 0.4, 0.35]) ] @testset verbose=true "$(names[i])" for i in eachindex(names) coords = coordinates[i] - nhs = GridNeighborhoodSearch{size(coords, 1)}(search_radius = 0.1, - n_points = size(coords, 2), - periodic_box = periodic_boxes[i]) - - initialize_grid!(nhs, coords) - - neighbors = [sort(collect(PointNeighbors.eachneighbor(coords[:, i], nhs))) - for i in 1:5] - - # Note that (1, 2) and (2, 3) are not neighbors, but they are in neighboring cells - @test neighbors[1] == [1, 2, 3, 5] - @test neighbors[2] == [1, 2, 3] - @test neighbors[3] == [1, 2, 3] - @test neighbors[4] == [4] - @test neighbors[5] == [1, 5] + # `DictionaryCellList` + nhs1 = GridNeighborhoodSearch{size(coords, 1)}(search_radius = 0.1, + n_points = size(coords, 2), + periodic_box = periodic_boxes[i]) + + # `FullGridCellList` + cell_list = FullGridCellList(; min_corner = periodic_boxes[i].min_corner, + max_corner = periodic_boxes[i].max_corner, + search_radius = 0.1) + nhs2 = GridNeighborhoodSearch{size(coords, 1)}(search_radius = 0.1, + n_points = size(coords, 2), + periodic_box = periodic_boxes[i], + cell_list = cell_list) + + for nhs in (nhs1, nhs2) + initialize_grid!(nhs, coords) + + neighbors = [sort(collect(PointNeighbors.eachneighbor(coords[:, i], nhs))) + for i in 1:5] + + # Note that (1, 2) and (2, 3) are not neighbors, but they are in neighboring cells + @test neighbors[1] == [1, 2, 3, 5] + @test neighbors[2] == [1, 2, 3] + @test neighbors[3] == [1, 2, 3] + @test neighbors[4] == [4] + @test neighbors[5] == [1, 5] + end end @testset "Offset Domain Triggering Split Cells" begin From fa3e5dc42dc60c34c26153c7ca52437d14055558 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Tue, 23 Jun 2026 13:13:12 +0200 Subject: [PATCH 3/4] Fix --- test/nhs_grid.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/test/nhs_grid.jl b/test/nhs_grid.jl index 3f4cbb29..7b891634 100644 --- a/test/nhs_grid.jl +++ b/test/nhs_grid.jl @@ -238,7 +238,7 @@ [-0.08 0.0 0.18 0.1 -0.08 -0.12 -0.05 -0.09 0.15 0.39], [-0.08 0.0 0.18 0.1 -0.08 - -0.12 -0.05 -0.09 0.15 0.42], + -0.12 -0.05 -0.09 0.15 0.42] .- [0.0017, 0.005], [-0.08 0.0 0.18 0.1 -0.08 -0.12 -0.05 -0.09 0.15 0.39 0.14 0.34 0.12 0.06 0.13] @@ -248,7 +248,11 @@ PeriodicBox(min_corner = [-0.1, -0.2], max_corner = [0.2, 0.4]), # The `GridNeighborhoodSearch` is forced to round up the cell sizes in this test # to avoid split cells. - PeriodicBox(min_corner = [-0.1, -0.205], max_corner = [0.205, 0.435]), + # In order for the points to have the same relative positions inside their cells + # with the two cell lists, we need to offset the coordinates and periodic box + # to make the grid aligned with the origin. + PeriodicBox(min_corner = [-0.1, -0.2] .- [0.0017, 0.005], + max_corner = [0.205, 0.43] .- [0.0017, 0.005]), PeriodicBox(min_corner = [-0.1, -0.2, 0.05], max_corner = [0.2, 0.4, 0.35]) ] From a55b4718145008d796039c5d4a84b7288ef2c93b Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Wed, 24 Jun 2026 12:41:34 +0200 Subject: [PATCH 4/4] Reformat --- test/nhs_grid.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/nhs_grid.jl b/test/nhs_grid.jl index 7b891634..9caabf8d 100644 --- a/test/nhs_grid.jl +++ b/test/nhs_grid.jl @@ -277,7 +277,7 @@ initialize_grid!(nhs, coords) neighbors = [sort(collect(PointNeighbors.eachneighbor(coords[:, i], nhs))) - for i in 1:5] + for i in 1:5] # Note that (1, 2) and (2, 3) are not neighbors, but they are in neighboring cells @test neighbors[1] == [1, 2, 3, 5]