From 82350ddcc007e2b157c399e1f0454229d986e4d9 Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Wed, 8 May 2024 15:41:25 +0200
Subject: [PATCH 01/17] Implement neighborhood search based on CellListMap.jl

---
 Project.toml           |  2 ++
 src/PointNeighbors.jl  |  1 +
 src/celllistmap_nhs.jl | 63 ++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 66 insertions(+)
 create mode 100644 src/celllistmap_nhs.jl

diff --git a/Project.toml b/Project.toml
index fc362260..b1cee123 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,12 +4,14 @@ authors = ["Erik Faulhaber <erik.faulhaber@uni-koeln.de>"]
 version = "0.2.2-dev"
 
 [deps]
+CellListMap = "69e1c6dd-3888-40e6-b3c8-31ac5f578864"
 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]
+CellListMap = "0.9"
 LinearAlgebra = "1"
 Polyester = "0.7.5"
 Reexport = "1"
diff --git a/src/PointNeighbors.jl b/src/PointNeighbors.jl
index 8c184006..0b7d896c 100644
--- a/src/PointNeighbors.jl
+++ b/src/PointNeighbors.jl
@@ -10,6 +10,7 @@ include("util.jl")
 include("neighborhood_search.jl")
 include("trivial_nhs.jl")
 include("grid_nhs.jl")
+include("celllistmap_nhs.jl")
 
 export for_particle_neighbor
 export TrivialNeighborhoodSearch, GridNeighborhoodSearch
diff --git a/src/celllistmap_nhs.jl b/src/celllistmap_nhs.jl
new file mode 100644
index 00000000..3e927263
--- /dev/null
+++ b/src/celllistmap_nhs.jl
@@ -0,0 +1,63 @@
+mutable struct CellListMapNeighborhoodSearch{CL, B}
+    cell_list :: CL
+    box       :: B
+
+    function CellListMapNeighborhoodSearch{NDIMS}(search_radius) where {NDIMS}
+        # Create a cell list with only one particle and resize it later
+        x = zeros(NDIMS, 1)
+        box = CellListMap.Box(CellListMap.limits(x, x), search_radius)
+        cell_list = CellListMap.CellList(x, x, box)
+
+        return new{typeof(cell_list), typeof(box)}(cell_list, box)
+    end
+end
+
+function initialize!(neighborhood_search::CellListMapNeighborhoodSearch,
+                     x::AbstractMatrix, y::AbstractMatrix)
+    update!(neighborhood_search, x, y)
+end
+
+function update!(neighborhood_search::CellListMapNeighborhoodSearch,
+                 x::AbstractMatrix, y::AbstractMatrix;
+                 particles_moving = (true, true))
+    (; cell_list) = neighborhood_search
+
+    # Resize box
+    box = CellListMap.Box(CellListMap.limits(x, y), neighborhood_search.box.cutoff)
+    neighborhood_search.box = box
+
+    # Resize and update cell list
+    CellListMap.UpdateCellList!(x, y, box, cell_list)
+
+    # Recalculate number of batches for multithreading
+    CellListMap.set_number_of_batches!(cell_list)
+
+    return neighborhood_search
+end
+
+# The type annotation is to make Julia specialize on the type of the function.
+# Otherwise, unspecialized code will cause a lot of allocations
+# and heavily impact performance.
+# See https://docs.julialang.org/en/v1/manual/performance-tips/#Be-aware-of-when-Julia-avoids-specializing
+function for_particle_neighbor(f::T, system_coords, neighbor_coords,
+                               neighborhood_search::CellListMapNeighborhoodSearch;
+                               particles = axes(system_coords, 2),
+                               parallel = true) where {T}
+    (; cell_list, box) = neighborhood_search
+
+    # `0` is the returned output, which we don't use
+    CellListMap.map_pairwise!(0, box, cell_list,
+                              parallel = parallel) do x, y, i, j, d2, output
+        # Skip all indices not in `particles`
+        i in particles || return output
+
+        pos_diff = x - y
+        distance = sqrt(d2)
+
+        @inline f(i, j, pos_diff, distance)
+
+        return output
+    end
+
+    return nothing
+end

From f7f48e8c8e6175c7e221d550da375f2d95e1659f Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Fri, 10 May 2024 09:05:12 +0200
Subject: [PATCH 02/17] using CellListMap

---
 src/PointNeighbors.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/PointNeighbors.jl b/src/PointNeighbors.jl
index 0b7d896c..43d20d2d 100644
--- a/src/PointNeighbors.jl
+++ b/src/PointNeighbors.jl
@@ -2,6 +2,7 @@ module PointNeighbors
 
 using Reexport: @reexport
 
+using CellListMap: CellListMap
 using LinearAlgebra: dot
 using Polyester: @batch
 @reexport using StaticArrays: SVector

From 357ebd252b7682d214e5431f8586d80b59dcab01 Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Thu, 4 Jul 2024 15:44:00 +0200
Subject: [PATCH 03/17] Add the new NHS as a package extension

---
 Project.toml                |  7 ++++-
 src/PointNeighbors.jl       |  5 ++-
 src/celllistmap_nhs.jl      | 63 -------------------------------------
 src/neighborhood_search.jl  |  2 ++
 test/Project.toml           |  2 ++
 test/neighborhood_search.jl |  3 ++
 test/test_util.jl           |  3 ++
 7 files changed, 18 insertions(+), 67 deletions(-)
 delete mode 100644 src/celllistmap_nhs.jl

diff --git a/Project.toml b/Project.toml
index 9f486e05..8a41d674 100644
--- a/Project.toml
+++ b/Project.toml
@@ -6,7 +6,6 @@ version = "0.4.1-dev"
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 Atomix = "a9b6321e-bd34-4604-b9c9-b65b8de01458"
-CellListMap = "69e1c6dd-3888-40e6-b3c8-31ac5f578864"
 GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -14,6 +13,12 @@ Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
+[weakdeps]
+CellListMap = "69e1c6dd-3888-40e6-b3c8-31ac5f578864"
+
+[extensions]
+PointNeighborsCellListMapExt = "CellListMap"
+
 [compat]
 Adapt = "4"
 Atomix = "0.1"
diff --git a/src/PointNeighbors.jl b/src/PointNeighbors.jl
index a6b30757..8314d158 100644
--- a/src/PointNeighbors.jl
+++ b/src/PointNeighbors.jl
@@ -4,7 +4,6 @@ using Reexport: @reexport
 
 using Adapt: Adapt
 using Atomix: Atomix
-using CellListMap: CellListMap
 using GPUArraysCore: AbstractGPUArray
 using KernelAbstractions: KernelAbstractions, @kernel, @index
 using LinearAlgebra: dot
@@ -18,11 +17,11 @@ include("nhs_trivial.jl")
 include("cell_lists/cell_lists.jl")
 include("nhs_grid.jl")
 include("nhs_precomputed.jl")
-include("celllistmap_nhs.jl")
 include("gpu.jl")
 
 export foreach_point_neighbor, foreach_neighbor
-export TrivialNeighborhoodSearch, GridNeighborhoodSearch, PrecomputedNeighborhoodSearch
+export TrivialNeighborhoodSearch, GridNeighborhoodSearch, PrecomputedNeighborhoodSearch,
+       CellListMapNeighborhoodSearch
 export DictionaryCellList, FullGridCellList
 export ParallelUpdate, SemiParallelUpdate, SerialUpdate
 export initialize!, update!, initialize_grid!, update_grid!
diff --git a/src/celllistmap_nhs.jl b/src/celllistmap_nhs.jl
deleted file mode 100644
index 3e927263..00000000
--- a/src/celllistmap_nhs.jl
+++ /dev/null
@@ -1,63 +0,0 @@
-mutable struct CellListMapNeighborhoodSearch{CL, B}
-    cell_list :: CL
-    box       :: B
-
-    function CellListMapNeighborhoodSearch{NDIMS}(search_radius) where {NDIMS}
-        # Create a cell list with only one particle and resize it later
-        x = zeros(NDIMS, 1)
-        box = CellListMap.Box(CellListMap.limits(x, x), search_radius)
-        cell_list = CellListMap.CellList(x, x, box)
-
-        return new{typeof(cell_list), typeof(box)}(cell_list, box)
-    end
-end
-
-function initialize!(neighborhood_search::CellListMapNeighborhoodSearch,
-                     x::AbstractMatrix, y::AbstractMatrix)
-    update!(neighborhood_search, x, y)
-end
-
-function update!(neighborhood_search::CellListMapNeighborhoodSearch,
-                 x::AbstractMatrix, y::AbstractMatrix;
-                 particles_moving = (true, true))
-    (; cell_list) = neighborhood_search
-
-    # Resize box
-    box = CellListMap.Box(CellListMap.limits(x, y), neighborhood_search.box.cutoff)
-    neighborhood_search.box = box
-
-    # Resize and update cell list
-    CellListMap.UpdateCellList!(x, y, box, cell_list)
-
-    # Recalculate number of batches for multithreading
-    CellListMap.set_number_of_batches!(cell_list)
-
-    return neighborhood_search
-end
-
-# The type annotation is to make Julia specialize on the type of the function.
-# Otherwise, unspecialized code will cause a lot of allocations
-# and heavily impact performance.
-# See https://docs.julialang.org/en/v1/manual/performance-tips/#Be-aware-of-when-Julia-avoids-specializing
-function for_particle_neighbor(f::T, system_coords, neighbor_coords,
-                               neighborhood_search::CellListMapNeighborhoodSearch;
-                               particles = axes(system_coords, 2),
-                               parallel = true) where {T}
-    (; cell_list, box) = neighborhood_search
-
-    # `0` is the returned output, which we don't use
-    CellListMap.map_pairwise!(0, box, cell_list,
-                              parallel = parallel) do x, y, i, j, d2, output
-        # Skip all indices not in `particles`
-        i in particles || return output
-
-        pos_diff = x - y
-        distance = sqrt(d2)
-
-        @inline f(i, j, pos_diff, distance)
-
-        return output
-    end
-
-    return nothing
-end
diff --git a/src/neighborhood_search.jl b/src/neighborhood_search.jl
index dacbc103..3299b983 100644
--- a/src/neighborhood_search.jl
+++ b/src/neighborhood_search.jl
@@ -236,3 +236,5 @@ end
 @inline function periodic_coords(coords, periodic_box::Nothing)
     return coords
 end
+
+CellListMapNeighborhoodSearch() = error("CellListMap.jl has to be imported to use this")
diff --git a/test/Project.toml b/test/Project.toml
index f57fdf2e..3e2bd507 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,10 +1,12 @@
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+CellListMap = "69e1c6dd-3888-40e6-b3c8-31ac5f578864"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 BenchmarkTools = "1"
+CellListMap = "0.9"
 Plots = "1"
 Test = "1"
diff --git a/test/neighborhood_search.jl b/test/neighborhood_search.jl
index a49ad2b0..c54d62c1 100644
--- a/test/neighborhood_search.jl
+++ b/test/neighborhood_search.jl
@@ -178,6 +178,7 @@
                                                                            search_radius,
                                                                            backend = Vector{Vector{Int}})),
                 PrecomputedNeighborhoodSearch{NDIMS}(; search_radius, n_points),
+                CellListMapNeighborhoodSearch(NDIMS, search_radius),
             ]
 
             names = [
@@ -187,6 +188,7 @@
                 "`GridNeighborhoodSearch` with `FullGridCellList` with `DynamicVectorOfVectors` and `SemiParallelUpdate`",
                 "`GridNeighborhoodSearch` with `FullGridCellList` with `Vector{Vector}`",
                 "`PrecomputedNeighborhoodSearch`",
+                "`CellListMapNeighborhoodSearch`",
             ]
 
             # Also test copied templates
@@ -202,6 +204,7 @@
                                                                            max_corner,
                                                                            backend = Vector{Vector{Int32}})),
                 PrecomputedNeighborhoodSearch{NDIMS}(),
+                CellListMapNeighborhoodSearch(NDIMS, 1.0),
             ]
             copied_nhs = copy_neighborhood_search.(template_nhs, search_radius, n_points)
             append!(neighborhood_searches, copied_nhs)
diff --git a/test/test_util.jl b/test/test_util.jl
index 07a3ac98..a5c19110 100644
--- a/test/test_util.jl
+++ b/test/test_util.jl
@@ -3,6 +3,9 @@
 using Test: @test, @testset, @test_throws
 using PointNeighbors
 
+# Load `PointNeighborsCellListMapExt`
+import CellListMap
+
 """
     @trixi_testset "name of the testset" #= code to test #=
 

From c8f1ce2d3e4163faab49e78a69ec57cfa93cb9e1 Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Thu, 4 Jul 2024 15:57:04 +0200
Subject: [PATCH 04/17] Add missing extension file

---
 ext/PointNeighborsCellListMapExt.jl | 78 +++++++++++++++++++++++++++++
 1 file changed, 78 insertions(+)
 create mode 100644 ext/PointNeighborsCellListMapExt.jl

diff --git a/ext/PointNeighborsCellListMapExt.jl b/ext/PointNeighborsCellListMapExt.jl
new file mode 100644
index 00000000..617a5a23
--- /dev/null
+++ b/ext/PointNeighborsCellListMapExt.jl
@@ -0,0 +1,78 @@
+module PointNeighborsCellListMapExt
+
+using PointNeighbors
+using CellListMap: CellListMap
+
+mutable struct CellListMapNeighborhoodSearch{CL, B}
+    cell_list :: CL
+    box       :: B
+
+    function PointNeighbors.CellListMapNeighborhoodSearch(NDIMS, search_radius)
+        # Create a cell list with only one point and resize it later
+        x = zeros(NDIMS, 1)
+        box = CellListMap.Box(CellListMap.limits(x, x), search_radius)
+        cell_list = CellListMap.CellList(x, x, box)
+
+        return new{typeof(cell_list), typeof(box)}(cell_list, box)
+    end
+end
+
+function PointNeighbors.initialize!(neighborhood_search::CellListMapNeighborhoodSearch,
+                                    x::AbstractMatrix, y::AbstractMatrix)
+    PointNeighbors.update!(neighborhood_search, x, y)
+end
+
+function PointNeighbors.update!(neighborhood_search::CellListMapNeighborhoodSearch,
+                                x::AbstractMatrix, y::AbstractMatrix;
+                                points_moving = (true, true))
+    (; cell_list) = neighborhood_search
+
+    # Resize box
+    box = CellListMap.Box(CellListMap.limits(x, y), neighborhood_search.box.cutoff)
+    neighborhood_search.box = box
+
+    # Resize and update cell list
+    CellListMap.UpdateCellList!(x, y, box, cell_list)
+
+    # Recalculate number of batches for multithreading
+    CellListMap.set_number_of_batches!(cell_list)
+
+    return neighborhood_search
+end
+
+# The type annotation is to make Julia specialize on the type of the function.
+# Otherwise, unspecialized code will cause a lot of allocations
+# and heavily impact performance.
+# See https://docs.julialang.org/en/v1/manual/performance-tips/#Be-aware-of-when-Julia-avoids-specializing
+function PointNeighbors.foreach_point_neighbor(f::T, system_coords, neighbor_coords,
+                                               neighborhood_search::CellListMapNeighborhoodSearch;
+                                               points = axes(system_coords, 2),
+                                               parallel = true) where {T}
+    (; cell_list, box) = neighborhood_search
+
+    # `0` is the returned output, which we don't use
+    CellListMap.map_pairwise!(0, box, cell_list,
+                              parallel = parallel) do x, y, i, j, d2, output
+        # Skip all indices not in `points`
+        i in points || return output
+
+        pos_diff = x - y
+        distance = sqrt(d2)
+
+        @inline f(i, j, pos_diff, distance)
+
+        return output
+    end
+
+    return nothing
+end
+
+function PointNeighbors.copy_neighborhood_search(nhs::CellListMapNeighborhoodSearch,
+                                                 search_radius, n_points;
+                                                 eachpoint = 1:n_points)
+    NDIMS = length(nhs.box.cell_size)
+
+    return PointNeighbors.CellListMapNeighborhoodSearch(NDIMS, search_radius)
+end
+
+end

From 1d1c078c34326c8f46e1e8489241c6aab130a279 Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Fri, 5 Jul 2024 09:41:01 +0200
Subject: [PATCH 05/17] Add functions `search_radius` and `ndims`

---
 ext/PointNeighborsCellListMapExt.jl | 12 +++++++++---
 1 file changed, 9 insertions(+), 3 deletions(-)

diff --git a/ext/PointNeighborsCellListMapExt.jl b/ext/PointNeighborsCellListMapExt.jl
index 617a5a23..5ccb7bc0 100644
--- a/ext/PointNeighborsCellListMapExt.jl
+++ b/ext/PointNeighborsCellListMapExt.jl
@@ -17,6 +17,14 @@ mutable struct CellListMapNeighborhoodSearch{CL, B}
     end
 end
 
+function PointNeighbors.search_radius(neighborhood_search::CellListMapNeighborhoodSearch)
+    return neighborhood_search.box.cutoff
+end
+
+function Base.ndims(neighborhood_search::CellListMapNeighborhoodSearch)
+    return length(neighborhood_search.box.cell_size)
+end
+
 function PointNeighbors.initialize!(neighborhood_search::CellListMapNeighborhoodSearch,
                                     x::AbstractMatrix, y::AbstractMatrix)
     PointNeighbors.update!(neighborhood_search, x, y)
@@ -70,9 +78,7 @@ end
 function PointNeighbors.copy_neighborhood_search(nhs::CellListMapNeighborhoodSearch,
                                                  search_radius, n_points;
                                                  eachpoint = 1:n_points)
-    NDIMS = length(nhs.box.cell_size)
-
-    return PointNeighbors.CellListMapNeighborhoodSearch(NDIMS, search_radius)
+    return PointNeighbors.CellListMapNeighborhoodSearch(ndims(nhs), search_radius)
 end
 
 end

From e3a96372aaebdbfdfdf7b5fc26cbe26c3f96f2f3 Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Tue, 26 Nov 2024 15:59:33 +0100
Subject: [PATCH 06/17] Fix `parallel` and error message

---
 ext/PointNeighborsCellListMapExt.jl | 5 +++--
 src/neighborhood_search.jl          | 4 +++-
 2 files changed, 6 insertions(+), 3 deletions(-)

diff --git a/ext/PointNeighborsCellListMapExt.jl b/ext/PointNeighborsCellListMapExt.jl
index 5ccb7bc0..99cedc88 100644
--- a/ext/PointNeighborsCellListMapExt.jl
+++ b/ext/PointNeighborsCellListMapExt.jl
@@ -58,9 +58,10 @@ function PointNeighbors.foreach_point_neighbor(f::T, system_coords, neighbor_coo
                                                parallel = true) where {T}
     (; cell_list, box) = neighborhood_search
 
-    # `0` is the returned output, which we don't use
+    # `0` is the returned output, which we don't use.
+    # Note that `parallel !== false` is `true` when `parallel` is a PointNeighbors backend.
     CellListMap.map_pairwise!(0, box, cell_list,
-                              parallel = parallel) do x, y, i, j, d2, output
+                              parallel = parallel !== false) do x, y, i, j, d2, output
         # Skip all indices not in `points`
         i in points || return output
 
diff --git a/src/neighborhood_search.jl b/src/neighborhood_search.jl
index 2ef11800..0dce3cf9 100644
--- a/src/neighborhood_search.jl
+++ b/src/neighborhood_search.jl
@@ -261,4 +261,6 @@ end
     return coords
 end
 
-CellListMapNeighborhoodSearch() = error("CellListMap.jl has to be imported to use this")
+function CellListMapNeighborhoodSearch(_...)
+    error("CellListMap.jl has to be imported to use `CellListMapNeighborhoodSearch`")
+end

From b874b5fd002add1251792a8cf5fc98148ac53cfa Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Mon, 2 Dec 2024 15:46:00 +0100
Subject: [PATCH 07/17] Use different CellListMap implementation when `points
 === neighbors`

---
 ext/PointNeighborsCellListMapExt.jl | 39 ++++++++++++++++++++++++++---
 1 file changed, 35 insertions(+), 4 deletions(-)

diff --git a/ext/PointNeighborsCellListMapExt.jl b/ext/PointNeighborsCellListMapExt.jl
index 99cedc88..ceab7b3d 100644
--- a/ext/PointNeighborsCellListMapExt.jl
+++ b/ext/PointNeighborsCellListMapExt.jl
@@ -1,17 +1,23 @@
 module PointNeighborsCellListMapExt
 
 using PointNeighbors
-using CellListMap: CellListMap
+using CellListMap: CellListMap, CellList, CellListPair
 
 mutable struct CellListMapNeighborhoodSearch{CL, B}
     cell_list :: CL
     box       :: B
 
-    function PointNeighbors.CellListMapNeighborhoodSearch(NDIMS, search_radius)
+    function PointNeighbors.CellListMapNeighborhoodSearch(NDIMS, search_radius,
+                                                          points_equal_neighbors = true)
         # Create a cell list with only one point and resize it later
         x = zeros(NDIMS, 1)
         box = CellListMap.Box(CellListMap.limits(x, x), search_radius)
-        cell_list = CellListMap.CellList(x, x, box)
+
+        if points_equal_neighbors
+            cell_list = CellListMap.CellList(x, box)
+        else
+            cell_list = CellListMap.CellList(x, x, box)
+        end
 
         return new{typeof(cell_list), typeof(box)}(cell_list, box)
     end
@@ -30,7 +36,8 @@ function PointNeighbors.initialize!(neighborhood_search::CellListMapNeighborhood
     PointNeighbors.update!(neighborhood_search, x, y)
 end
 
-function PointNeighbors.update!(neighborhood_search::CellListMapNeighborhoodSearch,
+# When `x !== y`, a `CellListPair` must be used
+function PointNeighbors.update!(neighborhood_search::CellListMapNeighborhoodSearch{<:CellListPair},
                                 x::AbstractMatrix, y::AbstractMatrix;
                                 points_moving = (true, true))
     (; cell_list) = neighborhood_search
@@ -48,6 +55,30 @@ function PointNeighbors.update!(neighborhood_search::CellListMapNeighborhoodSear
     return neighborhood_search
 end
 
+# When `points_equal_neighbors == true`, a `CellList` is used and `x` must be equal to `y`
+function PointNeighbors.update!(neighborhood_search::CellListMapNeighborhoodSearch{<:CellList},
+                                x::AbstractMatrix, y::AbstractMatrix;
+                                points_moving = (true, true))
+    (; cell_list) = neighborhood_search
+
+    @assert x === y "When `points_equal_neighbors == true`, `x` must be equal to `y`"
+
+    # Resize box
+    box = CellListMap.Box(CellListMap.limits(x), neighborhood_search.box.cutoff)
+    neighborhood_search.box = box
+
+    # Resize and update cell list
+    CellListMap.UpdateCellList!(x, box, cell_list)
+
+    # Recalculate number of batches for multithreading
+    CellListMap.set_number_of_batches!(cell_list)
+
+    # Due to https://github.com/m3g/CellListMap.jl/issues/106, we have to update again
+    CellListMap.UpdateCellList!(x, box, cell_list)
+
+    return neighborhood_search
+end
+
 # The type annotation is to make Julia specialize on the type of the function.
 # Otherwise, unspecialized code will cause a lot of allocations
 # and heavily impact performance.

From 0ba9dad009791cbf94354b379e12c56ff4d84305 Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Mon, 2 Dec 2024 17:18:46 +0100
Subject: [PATCH 08/17] Add docs, tests and fixes

---
 ext/PointNeighborsCellListMapExt.jl | 79 +++++++++++++++++++++++++++--
 src/neighborhood_search.jl          |  2 +-
 test/neighborhood_search.jl         | 12 +++--
 3 files changed, 84 insertions(+), 9 deletions(-)

diff --git a/ext/PointNeighborsCellListMapExt.jl b/ext/PointNeighborsCellListMapExt.jl
index ceab7b3d..2db7231e 100644
--- a/ext/PointNeighborsCellListMapExt.jl
+++ b/ext/PointNeighborsCellListMapExt.jl
@@ -3,12 +3,39 @@ module PointNeighborsCellListMapExt
 using PointNeighbors
 using CellListMap: CellListMap, CellList, CellListPair
 
+"""
+    CellListMapNeighborhoodSearch(NDIMS; search_radius = 1.0, points_equal_neighbors = false)
+
+Neighborhood search based on the package [CellListMap.jl](https://github.com/m3g/CellListMap.jl).
+This package provides a similar implementation to the [`GridNeighborhoodSearch`](@ref)
+with [`FullGridCellList`](@ref), but with better support for periodic boundaries.
+This is just a wrapper to use CellListMap.jl with the PointNeighbors.jl API.
+Note that periodic boundaries are not yet supported.
+
+# Arguments
+- `NDIMS`: Number of dimensions.
+
+# Keywords
+- `search_radius = 1.0`:    The fixed search radius. The default of `1.0` is useful together
+                            with [`copy_neighborhood_search`](@ref).
+- `points_equal_neighbors = false`: If `true`, a `CellListMap.CellList` is used instead of
+                                    a `CellListMap.CellListPair`. This requires that `x === y`
+                                    in [`initialize!`](@ref) and [`update!`](@ref).
+                                    This option exists only for benchmarking purposes.
+                                    It makes the main loop awkward because CellListMap.jl
+                                    only computes pairs with `i < j` and PointNeighbors.jl
+                                    computes all pairs, so we have to manually use symmetry
+                                    to add the missing pairs.
+
+!!! warning "Experimental implementation"
+    This is an experimental feature and may change in future releases.
+"""
 mutable struct CellListMapNeighborhoodSearch{CL, B}
     cell_list :: CL
     box       :: B
 
-    function PointNeighbors.CellListMapNeighborhoodSearch(NDIMS, search_radius,
-                                                          points_equal_neighbors = true)
+    function PointNeighbors.CellListMapNeighborhoodSearch(NDIMS; search_radius = 1.0,
+                                                          points_equal_neighbors = false)
         # Create a cell list with only one point and resize it later
         x = zeros(NDIMS, 1)
         box = CellListMap.Box(CellListMap.limits(x, x), search_radius)
@@ -84,7 +111,7 @@ end
 # and heavily impact performance.
 # See https://docs.julialang.org/en/v1/manual/performance-tips/#Be-aware-of-when-Julia-avoids-specializing
 function PointNeighbors.foreach_point_neighbor(f::T, system_coords, neighbor_coords,
-                                               neighborhood_search::CellListMapNeighborhoodSearch;
+                                               neighborhood_search::CellListMapNeighborhoodSearch{<:CellListPair};
                                                points = axes(system_coords, 2),
                                                parallel = true) where {T}
     (; cell_list, box) = neighborhood_search
@@ -107,10 +134,52 @@ function PointNeighbors.foreach_point_neighbor(f::T, system_coords, neighbor_coo
     return nothing
 end
 
-function PointNeighbors.copy_neighborhood_search(nhs::CellListMapNeighborhoodSearch,
+function PointNeighbors.foreach_point_neighbor(f::T, system_coords, neighbor_coords,
+                                               neighborhood_search::CellListMapNeighborhoodSearch{<:CellList};
+                                               points = axes(system_coords, 2),
+                                               parallel = true) where {T}
+    (; cell_list, box) = neighborhood_search
+
+    # `0` is the returned output, which we don't use.
+    # Note that `parallel !== false` is `true` when `parallel` is a PointNeighbors backend.
+    CellListMap.map_pairwise!(0, box, cell_list,
+                              parallel = parallel !== false) do x, y, i, j, d2, output
+        # Skip all indices not in `points`
+        i in points || return output
+
+        pos_diff = x - y
+        distance = sqrt(d2)
+
+        # When `points_equal_neighbors == true`, a `CellList` is used.
+        # With a `CellList`, we only see each pair once and have to use symmetry manually.
+        @inline f(i, j, pos_diff, distance)
+        @inline f(j, i, -pos_diff, distance)
+
+        return output
+    end
+
+    # With a `CellList`, interaction of a particle with itself is not included
+    PointNeighbors.@threaded system_coords for point in points
+        zero_pos_diff = zero(PointNeighbors.SVector{ndims(neighborhood_search),
+                                                    eltype(system_coords)})
+        @inline f(point, point, zero_pos_diff, zero(eltype(system_coords)))
+    end
+
+    return nothing
+end
+
+function PointNeighbors.copy_neighborhood_search(nhs::CellListMapNeighborhoodSearch{<:CellListPair},
+                                                 search_radius, n_points;
+                                                 eachpoint = 1:n_points)
+    return PointNeighbors.CellListMapNeighborhoodSearch(ndims(nhs); search_radius,
+                                                        points_equal_neighbors = false)
+end
+
+function PointNeighbors.copy_neighborhood_search(nhs::CellListMapNeighborhoodSearch{<:CellList},
                                                  search_radius, n_points;
                                                  eachpoint = 1:n_points)
-    return PointNeighbors.CellListMapNeighborhoodSearch(ndims(nhs), search_radius)
+    return PointNeighbors.CellListMapNeighborhoodSearch(ndims(nhs); search_radius,
+                                                        points_equal_neighbors = true)
 end
 
 end
diff --git a/src/neighborhood_search.jl b/src/neighborhood_search.jl
index 0dce3cf9..3551c4b9 100644
--- a/src/neighborhood_search.jl
+++ b/src/neighborhood_search.jl
@@ -261,6 +261,6 @@ end
     return coords
 end
 
-function CellListMapNeighborhoodSearch(_...)
+function CellListMapNeighborhoodSearch(args...; kwargs...)
     error("CellListMap.jl has to be imported to use `CellListMapNeighborhoodSearch`")
 end
diff --git a/test/neighborhood_search.jl b/test/neighborhood_search.jl
index a8d862a2..47d84b2a 100644
--- a/test/neighborhood_search.jl
+++ b/test/neighborhood_search.jl
@@ -91,6 +91,7 @@
             # Run this for every neighborhood search
             @testset "$(names[j])" for j in eachindex(names)
                 nhs = neighborhood_searches[j]
+                @test PointNeighbors.search_radius(nhs) == search_radius
 
                 initialize!(nhs, coords, coords)
 
@@ -176,7 +177,9 @@
                                                                            search_radius,
                                                                            backend = Vector{Vector{Int}})),
                 PrecomputedNeighborhoodSearch{NDIMS}(; search_radius, n_points),
-                CellListMapNeighborhoodSearch(NDIMS, search_radius)
+                CellListMapNeighborhoodSearch(NDIMS; search_radius),
+                CellListMapNeighborhoodSearch(NDIMS; search_radius,
+                                              points_equal_neighbors = true)
             ]
 
             names = [
@@ -186,7 +189,8 @@
                 "`GridNeighborhoodSearch` with `FullGridCellList` with `DynamicVectorOfVectors` and `SemiParallelUpdate`",
                 "`GridNeighborhoodSearch` with `FullGridCellList` with `Vector{Vector}`",
                 "`PrecomputedNeighborhoodSearch`",
-                "`CellListMapNeighborhoodSearch`"
+                "`CellListMapNeighborhoodSearch`",
+                "`CellListMapNeighborhoodSearch` with `points_equal_neighbors = true`"
             ]
 
             # Also test copied templates
@@ -202,7 +206,8 @@
                                                                            max_corner,
                                                                            backend = Vector{Vector{Int32}})),
                 PrecomputedNeighborhoodSearch{NDIMS}(),
-                CellListMapNeighborhoodSearch(NDIMS, 1.0)
+                CellListMapNeighborhoodSearch(NDIMS),
+                CellListMapNeighborhoodSearch(NDIMS, points_equal_neighbors = true)
             ]
             copied_nhs = copy_neighborhood_search.(template_nhs, search_radius, n_points)
             append!(neighborhood_searches, copied_nhs)
@@ -212,6 +217,7 @@
 
             @testset "$(names[i])" for i in eachindex(names)
                 nhs = neighborhood_searches[i]
+                @test PointNeighbors.search_radius(nhs) == search_radius
 
                 # Initialize with `seed = 1`
                 initialize!(nhs, coords_initialize, coords_initialize)

From aad3e13c8269b05ac1b6cff83d8ea31c2277aded Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Mon, 2 Dec 2024 17:50:17 +0100
Subject: [PATCH 09/17] Improve usage error messages

---
 src/neighborhood_search.jl | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/neighborhood_search.jl b/src/neighborhood_search.jl
index 3551c4b9..1adaf8e6 100644
--- a/src/neighborhood_search.jl
+++ b/src/neighborhood_search.jl
@@ -261,6 +261,8 @@ end
     return coords
 end
 
-function CellListMapNeighborhoodSearch(args...; kwargs...)
-    error("CellListMap.jl has to be imported to use `CellListMapNeighborhoodSearch`")
+function CellListMapNeighborhoodSearch(NDIMS; search_radius = 1.0,
+                                       points_equal_neighbors = false)
+    error("CellListMap.jl has to be imported to use `CellListMapNeighborhoodSearch`. " *
+          "If you did import CellListMap.jl, please try `]instantiate`.")
 end

From cbf25c599f3f241147b91cad10ccd7938e1545a2 Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Mon, 2 Dec 2024 18:21:00 +0100
Subject: [PATCH 10/17] Reformat code

---
 ext/PointNeighborsCellListMapExt.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/ext/PointNeighborsCellListMapExt.jl b/ext/PointNeighborsCellListMapExt.jl
index 2db7231e..e905f08a 100644
--- a/ext/PointNeighborsCellListMapExt.jl
+++ b/ext/PointNeighborsCellListMapExt.jl
@@ -88,7 +88,7 @@ function PointNeighbors.update!(neighborhood_search::CellListMapNeighborhoodSear
                                 points_moving = (true, true))
     (; cell_list) = neighborhood_search
 
-    @assert x === y "When `points_equal_neighbors == true`, `x` must be equal to `y`"
+    @assert x===y "When `points_equal_neighbors == true`, `x` must be equal to `y`"
 
     # Resize box
     box = CellListMap.Box(CellListMap.limits(x), neighborhood_search.box.cutoff)

From 7779a7e341db437b1da8038787d67a57131b50bb Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Tue, 3 Dec 2024 11:23:19 +0100
Subject: [PATCH 11/17] Fix method overwrite warning

---
 ext/PointNeighborsCellListMapExt.jl | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/ext/PointNeighborsCellListMapExt.jl b/ext/PointNeighborsCellListMapExt.jl
index e905f08a..fdadea2b 100644
--- a/ext/PointNeighborsCellListMapExt.jl
+++ b/ext/PointNeighborsCellListMapExt.jl
@@ -34,7 +34,9 @@ mutable struct CellListMapNeighborhoodSearch{CL, B}
     cell_list :: CL
     box       :: B
 
-    function PointNeighbors.CellListMapNeighborhoodSearch(NDIMS; search_radius = 1.0,
+    # Add dispatch on `NDIMS` to avoid method overwriting of the function in PointNeighbors.jl
+    function PointNeighbors.CellListMapNeighborhoodSearch(NDIMS::Integer;
+                                                          search_radius = 1.0,
                                                           points_equal_neighbors = false)
         # Create a cell list with only one point and resize it later
         x = zeros(NDIMS, 1)

From 4e4ec1d46ed47eb8fe1ccc3bb25962b946779e7f Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Wed, 4 Dec 2024 16:53:58 +0100
Subject: [PATCH 12/17] Add CellListMap.jl NHS to benchmark code

---
 benchmarks/plot.jl | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/benchmarks/plot.jl b/benchmarks/plot.jl
index 976cf045..e56585d0 100644
--- a/benchmarks/plot.jl
+++ b/benchmarks/plot.jl
@@ -1,5 +1,6 @@
 using Plots
 using BenchmarkTools
+import CellListMap
 
 # Generate a rectangular point cloud
 include("../test/point_cloud.jl")
@@ -40,7 +41,8 @@ function plot_benchmarks(benchmark, n_points_per_dimension, iterations;
                          seed = 1, perturbation_factor_position = 1.0)
     neighborhood_searches_names = ["TrivialNeighborhoodSearch";;
                                    "GridNeighborhoodSearch";;
-                                   "PrecomputedNeighborhoodSearch"]
+                                   "PrecomputedNeighborhoodSearch";;
+                                   "CellListMapNeighborhoodSearch"]
 
     # Multiply number of points in each iteration (roughly) by this factor
     scaling_factor = 4
@@ -62,7 +64,8 @@ function plot_benchmarks(benchmark, n_points_per_dimension, iterations;
         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)
+            PrecomputedNeighborhoodSearch{NDIMS}(; search_radius, n_points = n_particles),
+            CellListMapNeighborhoodSearch(NDIMS; search_radius)
         ]
 
         for i in eachindex(neighborhood_searches)

From 39ac331d3797cfeef2754180b583ccbc9f2eea30 Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Fri, 6 Dec 2024 12:53:36 +0100
Subject: [PATCH 13/17] Add comments

---
 ext/PointNeighborsCellListMapExt.jl | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/ext/PointNeighborsCellListMapExt.jl b/ext/PointNeighborsCellListMapExt.jl
index fdadea2b..94c934ce 100644
--- a/ext/PointNeighborsCellListMapExt.jl
+++ b/ext/PointNeighborsCellListMapExt.jl
@@ -32,6 +32,7 @@ Note that periodic boundaries are not yet supported.
 """
 mutable struct CellListMapNeighborhoodSearch{CL, B}
     cell_list :: CL
+    # Note that we need this struct to be mutable to replace the box in `update!`
     box       :: B
 
     # Add dispatch on `NDIMS` to avoid method overwriting of the function in PointNeighbors.jl
@@ -160,7 +161,8 @@ function PointNeighbors.foreach_point_neighbor(f::T, system_coords, neighbor_coo
         return output
     end
 
-    # With a `CellList`, interaction of a particle with itself is not included
+    # With a `CellList`, only pairs with `i < j` are considered.
+    # We can cover `i > j` with symmetry above, but `i = j` has to be computed separately.
     PointNeighbors.@threaded system_coords for point in points
         zero_pos_diff = zero(PointNeighbors.SVector{ndims(neighborhood_search),
                                                     eltype(system_coords)})

From fbf2e993006ef836982d8fc199cff98a81e5a72d Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Fri, 6 Dec 2024 12:57:20 +0100
Subject: [PATCH 14/17] Fix benchmark tests

---
 benchmarks/plot.jl | 16 ++++++++++++----
 1 file changed, 12 insertions(+), 4 deletions(-)

diff --git a/benchmarks/plot.jl b/benchmarks/plot.jl
index e56585d0..c5e1323c 100644
--- a/benchmarks/plot.jl
+++ b/benchmarks/plot.jl
@@ -41,8 +41,12 @@ function plot_benchmarks(benchmark, n_points_per_dimension, iterations;
                          seed = 1, perturbation_factor_position = 1.0)
     neighborhood_searches_names = ["TrivialNeighborhoodSearch";;
                                    "GridNeighborhoodSearch";;
-                                   "PrecomputedNeighborhoodSearch";;
-                                   "CellListMapNeighborhoodSearch"]
+                                   "PrecomputedNeighborhoodSearch"]
+
+    if length(n_points_per_dimension) > 1
+        # Not implemented for 1D
+        push!(neighborhood_searches_names, "CellListMapNeighborhoodSearch")
+    end
 
     # Multiply number of points in each iteration (roughly) by this factor
     scaling_factor = 4
@@ -64,10 +68,14 @@ function plot_benchmarks(benchmark, n_points_per_dimension, iterations;
         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),
-            CellListMapNeighborhoodSearch(NDIMS; search_radius)
+            PrecomputedNeighborhoodSearch{NDIMS}(; search_radius, n_points = n_particles)
         ]
 
+        if NDIMS > 1
+            # Not implemented for 1D
+            push!(neighborhood_searches, CellListMapNeighborhoodSearch(NDIMS; search_radius))
+        end
+
         for i in eachindex(neighborhood_searches)
             neighborhood_search = neighborhood_searches[i]
             initialize!(neighborhood_search, coordinates, coordinates)

From d7000c8867f4ff8bd75e540537c8be4ae76bdee7 Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Mon, 9 Dec 2024 08:19:48 +0100
Subject: [PATCH 15/17] Fix plotting code

---
 benchmarks/plot.jl | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/benchmarks/plot.jl b/benchmarks/plot.jl
index c5e1323c..c64ea56d 100644
--- a/benchmarks/plot.jl
+++ b/benchmarks/plot.jl
@@ -39,8 +39,8 @@ 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";;
+    neighborhood_searches_names = ["TrivialNeighborhoodSearch",
+                                   "GridNeighborhoodSearch",
                                    "PrecomputedNeighborhoodSearch"]
 
     if length(n_points_per_dimension) > 1
@@ -88,10 +88,12 @@ function plot_benchmarks(benchmark, n_points_per_dimension, iterations;
         end
     end
 
+    labels = reshape(neighborhood_searches_names, (1, length(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)
+         label = labels, title = title)
 end

From 1326b5238b4606a7ec9ea94e2df868f3c0e475b7 Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Mon, 9 Dec 2024 08:21:03 +0100
Subject: [PATCH 16/17] Reformat

---
 benchmarks/plot.jl                  | 10 ++++++----
 ext/PointNeighborsCellListMapExt.jl |  4 ++--
 2 files changed, 8 insertions(+), 6 deletions(-)

diff --git a/benchmarks/plot.jl b/benchmarks/plot.jl
index c64ea56d..5832984c 100644
--- a/benchmarks/plot.jl
+++ b/benchmarks/plot.jl
@@ -39,9 +39,10 @@ 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 = [
+        "TrivialNeighborhoodSearch",
+        "GridNeighborhoodSearch",
+        "PrecomputedNeighborhoodSearch"]
 
     if length(n_points_per_dimension) > 1
         # Not implemented for 1D
@@ -73,7 +74,8 @@ function plot_benchmarks(benchmark, n_points_per_dimension, iterations;
 
         if NDIMS > 1
             # Not implemented for 1D
-            push!(neighborhood_searches, CellListMapNeighborhoodSearch(NDIMS; search_radius))
+            push!(neighborhood_searches,
+                  CellListMapNeighborhoodSearch(NDIMS; search_radius))
         end
 
         for i in eachindex(neighborhood_searches)
diff --git a/ext/PointNeighborsCellListMapExt.jl b/ext/PointNeighborsCellListMapExt.jl
index 94c934ce..7b290aa7 100644
--- a/ext/PointNeighborsCellListMapExt.jl
+++ b/ext/PointNeighborsCellListMapExt.jl
@@ -31,9 +31,9 @@ Note that periodic boundaries are not yet supported.
     This is an experimental feature and may change in future releases.
 """
 mutable struct CellListMapNeighborhoodSearch{CL, B}
-    cell_list :: CL
+    cell_list::CL
     # Note that we need this struct to be mutable to replace the box in `update!`
-    box       :: B
+    box::B
 
     # Add dispatch on `NDIMS` to avoid method overwriting of the function in PointNeighbors.jl
     function PointNeighbors.CellListMapNeighborhoodSearch(NDIMS::Integer;

From d183c1429ff9bb4846ce20cc4ed625c32444d014 Mon Sep 17 00:00:00 2001
From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com>
Date: Mon, 9 Dec 2024 09:45:22 +0100
Subject: [PATCH 17/17] Compute coverage for `ext` directory

---
 .github/workflows/ci.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 00906f7d..ebc6dff6 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -80,7 +80,7 @@ jobs:
         if: matrix.os == 'ubuntu-latest' && matrix.version == '1'
         uses: julia-actions/julia-processcoverage@v1
         with:
-          directories: src,test
+          directories: src,ext,test
       - name: Upload coverage report to Codecov
         # Only run coverage in one Job (Ubuntu and latest Julia version)
         if: matrix.os == 'ubuntu-latest' && matrix.version == '1'