From b7a96e9a578911f9748b0478357afa215448f08d Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Wed, 14 May 2025 18:12:50 +0200 Subject: [PATCH 1/6] Add option to skip inactive points in `foreach_point_neighbor` --- src/neighborhood_search.jl | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/src/neighborhood_search.jl b/src/neighborhood_search.jl index 9fb6dcbb..c19b256d 100644 --- a/src/neighborhood_search.jl +++ b/src/neighborhood_search.jl @@ -171,7 +171,8 @@ See also [`initialize!`](@ref), [`update!`](@ref). """ function foreach_point_neighbor(f::T, system_coords, neighbor_coords, neighborhood_search; parallelization_backend::ParallelizationBackend = default_backend(system_coords), - points = axes(system_coords, 2)) where {T} + points = axes(system_coords, 2), + points_active = nothing) where {T} # The type annotation above is to make Julia specialize on the type of the function. # Otherwise, unspecialized code will cause a lot of allocations # and heavily impact performance. @@ -181,14 +182,20 @@ function foreach_point_neighbor(f::T, system_coords, neighbor_coords, neighborho @boundscheck checkbounds(system_coords, ndims(neighborhood_search), points) @threaded parallelization_backend for point in points - # Now we can safely assume that `point` is inbounds - @inbounds foreach_neighbor(f, system_coords, neighbor_coords, - neighborhood_search, point) + # This check is optimized away when `points_active` is `nothing` + if point_active(points_active, point) + # After the explicit boundscheck, we can safely assume that `point` is inbounds + @inbounds foreach_neighbor(f, system_coords, neighbor_coords, + neighborhood_search, point) + end end return nothing end +@inline point_active(::Nothing, _) = true +@inline point_active(points_active, point) = points_active[point] + @propagate_inbounds function foreach_neighbor(f, system_coords, neighbor_system_coords, neighborhood_search::AbstractNeighborhoodSearch, point; From 5c8b76cd88e6010ba1b2027febe8ffb99222b4a1 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Wed, 14 May 2025 18:17:29 +0200 Subject: [PATCH 2/6] Add kwarg `points_active` to `update!` and `initialize!` --- src/neighborhood_search.jl | 4 +-- src/nhs_grid.jl | 52 ++++++++++++++++++++++---------------- src/nhs_precomputed.jl | 14 +++++----- src/nhs_trivial.jl | 4 +-- 4 files changed, 42 insertions(+), 32 deletions(-) diff --git a/src/neighborhood_search.jl b/src/neighborhood_search.jl index c19b256d..8d431ed4 100644 --- a/src/neighborhood_search.jl +++ b/src/neighborhood_search.jl @@ -35,7 +35,7 @@ See also [`update!`](@ref). """ @inline function initialize!(search::AbstractNeighborhoodSearch, x, y; parallelization_backend = default_backend(x), - eachindex_y = axes(y, 2)) + eachindex_y = axes(y, 2), points_active = nothing) return search end @@ -73,7 +73,7 @@ See also [`initialize!`](@ref). @inline function update!(search::AbstractNeighborhoodSearch, x, y; points_moving = (true, true), parallelization_backend = default_backend(x), - eachindex_y = axes(y, 2)) + eachindex_y = axes(y, 2), points_active = nothing) return search end diff --git a/src/nhs_grid.jl b/src/nhs_grid.jl index 08139ba9..418e7cd7 100644 --- a/src/nhs_grid.jl +++ b/src/nhs_grid.jl @@ -201,7 +201,7 @@ end function initialize_grid!(neighborhood_search::GridNeighborhoodSearch, y::AbstractMatrix; parallelization_backend = default_backend(y), - eachindex_y = axes(y, 2)) + eachindex_y = axes(y, 2), points_active = nothing) (; cell_list) = neighborhood_search empty!(cell_list) @@ -216,12 +216,15 @@ function initialize_grid!(neighborhood_search::GridNeighborhoodSearch, y::Abstra # Ignore the parallelization backend here. This cannot be parallelized. for point in eachindex_y - # Get cell index of the point's cell - point_coords = @inbounds extract_svector(y, Val(ndims(neighborhood_search)), point) - cell = cell_coords(point_coords, neighborhood_search) - - # Add point to corresponding cell - push_cell!(cell_list, cell, point) + # This check is optimized away when `points_active` is `nothing` + if point_active(points_active, point) + # Get cell index of the point's cell + point_coords = @inbounds extract_svector(y, Val(ndims(neighborhood_search)), point) + cell = cell_coords(point_coords, neighborhood_search) + + # Add point to corresponding cell + push_cell!(cell_list, cell, point) + end end return neighborhood_search @@ -230,7 +233,7 @@ end function initialize_grid!(neighborhood_search::GridNeighborhoodSearch{<:Any, ParallelUpdate}, y::AbstractMatrix; parallelization_backend = default_backend(y), - eachindex_y = axes(y, 2)) + eachindex_y = axes(y, 2), points_active = nothing) (; cell_list) = neighborhood_search empty!(cell_list) @@ -244,12 +247,15 @@ function initialize_grid!(neighborhood_search::GridNeighborhoodSearch{<:Any, @boundscheck checkbounds(y, eachindex_y) @threaded parallelization_backend for point in eachindex_y - # Get cell index of the point's cell - point_coords = @inbounds extract_svector(y, Val(ndims(neighborhood_search)), point) - cell = cell_coords(point_coords, neighborhood_search) - - # Add point to corresponding cell - push_cell_atomic!(cell_list, cell, point) + # This check is optimized away when `points_active` is `nothing` + if point_active(points_active, point) + # Get cell index of the point's cell + point_coords = @inbounds extract_svector(y, Val(ndims(neighborhood_search)), point) + cell = cell_coords(point_coords, neighborhood_search) + + # Add point to corresponding cell + push_cell_atomic!(cell_list, cell, point) + end end return neighborhood_search @@ -258,12 +264,13 @@ end function update!(neighborhood_search::GridNeighborhoodSearch, x::AbstractMatrix, y::AbstractMatrix; points_moving = (true, true), parallelization_backend = default_backend(x), - eachindex_y = axes(y, 2)) + eachindex_y = axes(y, 2), points_active = nothing) # 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 - update_grid!(neighborhood_search, y; eachindex_y, parallelization_backend) + update_grid!(neighborhood_search, y; eachindex_y, points_active, + parallelization_backend) end # Update only with neighbor coordinates @@ -273,10 +280,10 @@ function update_grid!(neighborhood_search::Union{GridNeighborhoodSearch{<:Any, SemiParallelUpdate}}, y::AbstractMatrix; parallelization_backend = default_backend(y), - eachindex_y = axes(y, 2)) + eachindex_y = axes(y, 2), points_active = nothing) (; cell_list, update_buffer) = neighborhood_search - if eachindex_y != axes(y, 2) + if eachindex_y != axes(y, 2) || points_active !== nothing # Incremental update doesn't support inactive points error("this neighborhood search/update strategy does not support inactive points") end @@ -381,10 +388,10 @@ end function update_grid!(neighborhood_search::GridNeighborhoodSearch{<:Any, ParallelIncrementalUpdate}, y::AbstractMatrix; parallelization_backend = default_backend(y), - eachindex_y = axes(y, 2)) + eachindex_y = axes(y, 2), points_active = nothing) (; cell_list, update_buffer) = neighborhood_search - if eachindex_y != axes(y, 2) + if eachindex_y != axes(y, 2) || points_active !== nothing # Incremental update doesn't support inactive points error("this neighborhood search/update strategy does not support inactive points") end @@ -445,8 +452,9 @@ function update_grid!(neighborhood_search::Union{GridNeighborhoodSearch{<:Any, GridNeighborhoodSearch{<:Any, SerialUpdate}}, y::AbstractMatrix; parallelization_backend = default_backend(y), - eachindex_y = axes(y, 2)) - initialize_grid!(neighborhood_search, y; parallelization_backend, eachindex_y) + eachindex_y = axes(y, 2), points_active = nothing) + initialize_grid!(neighborhood_search, y; parallelization_backend, + eachindex_y, points_active) end # Specialized version of the function in `neighborhood_search.jl`, which is faster diff --git a/src/nhs_precomputed.jl b/src/nhs_precomputed.jl index 69749f1c..98417642 100644 --- a/src/nhs_precomputed.jl +++ b/src/nhs_precomputed.jl @@ -72,21 +72,22 @@ end function update!(search::PrecomputedNeighborhoodSearch, x::AbstractMatrix, y::AbstractMatrix; points_moving = (true, true), parallelization_backend = default_backend(x), - eachindex_y = axes(y, 2)) + eachindex_y = axes(y, 2), points_active = nothing) (; neighborhood_search, neighbor_lists) = search # Update grid NHS - update!(neighborhood_search, x, y; eachindex_y, points_moving, parallelization_backend) + update!(neighborhood_search, x, y; eachindex_y, points_moving, points_active, + parallelization_backend) # Skip update if both point sets are static if any(points_moving) initialize_neighbor_lists!(neighbor_lists, neighborhood_search, x, y, - parallelization_backend, eachindex_y) + parallelization_backend, eachindex_y, points_active) end end function initialize_neighbor_lists!(neighbor_lists, neighborhood_search, x, y, - parallelization_backend, eachindex_y) + parallelization_backend, eachindex_y, points_active) # Initialize neighbor lists empty!(neighbor_lists) resize!(neighbor_lists, size(x, 2)) @@ -95,8 +96,9 @@ function initialize_neighbor_lists!(neighbor_lists, neighborhood_search, x, y, end # Fill neighbor lists - foreach_point_neighbor(x, y, neighborhood_search; parallelization_backend, - points = eachindex_y) do point, neighbor, _, _ + foreach_point_neighbor(x, y, neighborhood_search; + parallelization_backend, points = eachindex_y, + points_active) do point, neighbor, _, _ push!(neighbor_lists[point], neighbor) end end diff --git a/src/nhs_trivial.jl b/src/nhs_trivial.jl index 814ac03a..584564ad 100644 --- a/src/nhs_trivial.jl +++ b/src/nhs_trivial.jl @@ -34,14 +34,14 @@ end @inline function initialize!(search::TrivialNeighborhoodSearch, x, y; parallelization_backend = default_backend(x), - eachindex_y = axes(y, 2)) + eachindex_y = axes(y, 2), points_active = nothing) return search end @inline function update!(search::TrivialNeighborhoodSearch, x, y; points_moving = (true, true), parallelization_backend = default_backend(x), - eachindex_y = axes(y, 2)) + eachindex_y = axes(y, 2), points_active = nothing) return search end From 45c5731d4eb7348d128b05f5e90acbbf84ce59ea Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Wed, 14 May 2025 18:54:57 +0200 Subject: [PATCH 3/6] Update src/neighborhood_search.jl Co-authored-by: Niklas Neher <73897120+LasNikas@users.noreply.github.com> --- src/neighborhood_search.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/neighborhood_search.jl b/src/neighborhood_search.jl index 8d431ed4..7fc9678f 100644 --- a/src/neighborhood_search.jl +++ b/src/neighborhood_search.jl @@ -194,7 +194,7 @@ function foreach_point_neighbor(f::T, system_coords, neighbor_coords, neighborho end @inline point_active(::Nothing, _) = true -@inline point_active(points_active, point) = points_active[point] +@inline point_active(points_active, point) = points_active[point] == true @propagate_inbounds function foreach_neighbor(f, system_coords, neighbor_system_coords, neighborhood_search::AbstractNeighborhoodSearch, From 6dab97f46c9a01abf57459bf261f96fcba5a7c17 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Wed, 14 May 2025 18:56:17 +0200 Subject: [PATCH 4/6] Fix tests --- src/nhs_precomputed.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/nhs_precomputed.jl b/src/nhs_precomputed.jl index 98417642..1f3acf73 100644 --- a/src/nhs_precomputed.jl +++ b/src/nhs_precomputed.jl @@ -54,14 +54,14 @@ end function initialize!(search::PrecomputedNeighborhoodSearch, x::AbstractMatrix, y::AbstractMatrix; parallelization_backend = default_backend(x), - eachindex_y = axes(y, 2)) + eachindex_y = axes(y, 2), points_active = nothing) (; neighborhood_search, neighbor_lists) = search # Initialize grid NHS initialize!(neighborhood_search, x, y; eachindex_y, parallelization_backend) initialize_neighbor_lists!(neighbor_lists, neighborhood_search, x, y, - parallelization_backend, eachindex_y) + parallelization_backend, eachindex_y, points_active) end # WARNING! Experimental feature: From 546637a4111f7c69312e94ee89ac4e2e3041f0c6 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Wed, 14 May 2025 18:56:40 +0200 Subject: [PATCH 5/6] Reformat --- src/neighborhood_search.jl | 2 +- src/nhs_grid.jl | 6 ++++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/neighborhood_search.jl b/src/neighborhood_search.jl index 7fc9678f..36427896 100644 --- a/src/neighborhood_search.jl +++ b/src/neighborhood_search.jl @@ -186,7 +186,7 @@ function foreach_point_neighbor(f::T, system_coords, neighbor_coords, neighborho if point_active(points_active, point) # After the explicit boundscheck, we can safely assume that `point` is inbounds @inbounds foreach_neighbor(f, system_coords, neighbor_coords, - neighborhood_search, point) + neighborhood_search, point) end end diff --git a/src/nhs_grid.jl b/src/nhs_grid.jl index 418e7cd7..f33bccab 100644 --- a/src/nhs_grid.jl +++ b/src/nhs_grid.jl @@ -219,7 +219,8 @@ function initialize_grid!(neighborhood_search::GridNeighborhoodSearch, y::Abstra # This check is optimized away when `points_active` is `nothing` if point_active(points_active, point) # Get cell index of the point's cell - point_coords = @inbounds extract_svector(y, Val(ndims(neighborhood_search)), point) + point_coords = @inbounds extract_svector(y, Val(ndims(neighborhood_search)), + point) cell = cell_coords(point_coords, neighborhood_search) # Add point to corresponding cell @@ -250,7 +251,8 @@ function initialize_grid!(neighborhood_search::GridNeighborhoodSearch{<:Any, # This check is optimized away when `points_active` is `nothing` if point_active(points_active, point) # Get cell index of the point's cell - point_coords = @inbounds extract_svector(y, Val(ndims(neighborhood_search)), point) + point_coords = @inbounds extract_svector(y, Val(ndims(neighborhood_search)), + point) cell = cell_coords(point_coords, neighborhood_search) # Add point to corresponding cell From 22e543dfe49dd441f5ac98a95543c16c43051431 Mon Sep 17 00:00:00 2001 From: Erik Faulhaber <44124897+efaulhaber@users.noreply.github.com> Date: Thu, 15 May 2025 10:11:34 +0200 Subject: [PATCH 6/6] Fix --- src/nhs_grid.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/nhs_grid.jl b/src/nhs_grid.jl index f33bccab..ed878830 100644 --- a/src/nhs_grid.jl +++ b/src/nhs_grid.jl @@ -195,8 +195,9 @@ end function initialize!(neighborhood_search::GridNeighborhoodSearch, x::AbstractMatrix, y::AbstractMatrix; parallelization_backend = default_backend(x), - eachindex_y = axes(y, 2)) - initialize_grid!(neighborhood_search, y; parallelization_backend, eachindex_y) + eachindex_y = axes(y, 2), points_active = nothing) + initialize_grid!(neighborhood_search, y; parallelization_backend, + eachindex_y, points_active) end function initialize_grid!(neighborhood_search::GridNeighborhoodSearch, y::AbstractMatrix;