Fix bug due to rounding errors in periodicity computation

src/cell_lists/full_grid.jl

@@ -81,33 +81,11 @@ 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
+@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)

src/nhs_grid.jl

@@ -586,37 +586,44 @@ 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 `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 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
 
 @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_list, 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
+    return nonperiodic_cell_coords(coords, cell_list, cell_size)
+end
 
-    # 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
+@inline function nonperiodic_cell_coords(coords, cell_list, cell_size)
+    return Tuple(floor_to_int.(coords ./ cell_size))
 end
 
 function copy_neighborhood_search(nhs::GridNeighborhoodSearch, search_radius, n_points;

test/nhs_grid.jl

@@ -62,25 +62,24 @@
         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
         # bounds for the finite domain of the full grid cell list.
         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
@@ -239,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]
@@ -249,28 +248,44 @@
             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]),
+            # 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])
         ]
 
         @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