Skip to content

Commit 85be8bc

Browse files
efaulhaberefaulhaber
authored
Implement dynamic vector of vectors data structure (#34)
* Implement dynamic vector of vectors data structure * Add comments * Test for internal size
1 parent 927120c commit 85be8bc

7 files changed

Lines changed: 213 additions & 49 deletions

File tree

‎src/PointNeighbors.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@ using Polyester: @batch
77
@reexport using StaticArrays: SVector
88

99
include("util.jl")
10+
include("vector_of_vectors.jl")
1011
include("neighborhood_search.jl")
1112
include("nhs_trivial.jl")
1213
include("cell_lists/cell_lists.jl")

‎src/nhs_grid.jl‎

Lines changed: 32 additions & 45 deletions
Original file line numberDiff line numberDiff line change
@@ -51,24 +51,25 @@ since not sorting makes our implementation a lot faster (although less paralleli
5151
In: Computer Graphics Forum 30.1 (2011), pages 99–112.
5252
[doi: 10.1111/J.1467-8659.2010.01832.X](https://doi.org/10.1111/J.1467-8659.2010.01832.X)
5353
"""
54-
struct GridNeighborhoodSearch{NDIMS, ELTYPE, CL, CB, PB} <: AbstractNeighborhoodSearch
55-
cell_list :: CL
56-
search_radius :: ELTYPE
57-
periodic_box :: PB
58-
n_cells :: NTuple{NDIMS, Int} # Required to calculate periodic cell index
59-
cell_size :: NTuple{NDIMS, ELTYPE} # Required to calculate cell index
60-
cell_buffer :: CB # Multithreaded buffer for `update!`
61-
cell_buffer_indices :: Vector{Int} # Store which entries of `cell_buffer` are initialized
62-
threaded_update :: Bool
54+
struct GridNeighborhoodSearch{NDIMS, ELTYPE, CL, PB, UB} <: AbstractNeighborhoodSearch
55+
cell_list :: CL
56+
search_radius :: ELTYPE
57+
periodic_box :: PB
58+
n_cells :: NTuple{NDIMS, Int} # Required to calculate periodic cell index
59+
cell_size :: NTuple{NDIMS, ELTYPE} # Required to calculate cell index
60+
update_buffer :: UB # Multithreaded buffer for `update!`
61+
threaded_update :: Bool
6362

6463
function GridNeighborhoodSearch{NDIMS}(; search_radius = 0.0, n_points = 0,
6564
periodic_box = nothing,
6665
cell_list = DictionaryCellList{NDIMS}(),
6766
threaded_update = true) where {NDIMS}
6867
ELTYPE = typeof(search_radius)
6968

70-
cell_buffer = Array{index_type(cell_list), 2}(undef, n_points, Threads.nthreads())
71-
cell_buffer_indices = zeros(Int, Threads.nthreads())
69+
# Create update buffer and initialize it with empty vectors
70+
update_buffer = DynamicVectorOfVectors{index_type(cell_list)}(max_outer_length = Threads.nthreads(),
71+
max_inner_length = n_points)
72+
push!(update_buffer, (NTuple{NDIMS, Int}[] for _ in 1:Threads.nthreads())...)
7273

7374
if search_radius < eps() || isnothing(periodic_box)
7475
# No periodicity
@@ -90,37 +91,28 @@ struct GridNeighborhoodSearch{NDIMS, ELTYPE, CL, CB, PB} <: AbstractNeighborhood
9091
end
9192
end
9293

93-
new{NDIMS, ELTYPE, typeof(cell_list), typeof(cell_buffer),
94-
typeof(periodic_box)}(cell_list, search_radius, periodic_box, n_cells,
95-
cell_size, cell_buffer, cell_buffer_indices,
96-
threaded_update)
94+
new{NDIMS, ELTYPE, typeof(cell_list), typeof(periodic_box),
95+
typeof(update_buffer)}(cell_list, search_radius, periodic_box, n_cells,
96+
cell_size, update_buffer, threaded_update)
9797
end
9898
end
9999

100100
@inline Base.ndims(::GridNeighborhoodSearch{NDIMS}) where {NDIMS} = NDIMS
101101

102-
@inline function npoints(neighborhood_search::GridNeighborhoodSearch)
103-
return size(neighborhood_search.cell_buffer, 1)
104-
end
105-
106102
function initialize!(neighborhood_search::GridNeighborhoodSearch,
107103
x::AbstractMatrix, y::AbstractMatrix)
108104
initialize_grid!(neighborhood_search, y)
109105
end
110106

111-
function initialize_grid!(neighborhood_search::GridNeighborhoodSearch{NDIMS},
112-
y::AbstractMatrix) where {NDIMS}
113-
initialize_grid!(neighborhood_search, i -> extract_svector(y, Val(NDIMS), i))
114-
end
115-
116-
function initialize_grid!(neighborhood_search::GridNeighborhoodSearch, coords_fun)
107+
function initialize_grid!(neighborhood_search::GridNeighborhoodSearch, y::AbstractMatrix)
117108
(; cell_list) = neighborhood_search
118109

119110
empty!(cell_list)
120111

121-
for point in 1:npoints(neighborhood_search)
112+
for point in axes(y, 2)
122113
# Get cell index of the point's cell
123-
cell = cell_coords(coords_fun(point), neighborhood_search)
114+
point_coords = extract_svector(y, Val(ndims(neighborhood_search)), point)
115+
cell = cell_coords(point_coords, neighborhood_search)
124116

125117
# Add point to corresponding cell
126118
push_cell!(cell_list, cell, point)
@@ -147,20 +139,19 @@ end
147139

148140
# Modify the existing hash table by moving points into their new cells
149141
function update_grid!(neighborhood_search::GridNeighborhoodSearch, coords_fun)
150-
(; cell_list, cell_buffer, cell_buffer_indices, threaded_update) = neighborhood_search
142+
(; cell_list, update_buffer, threaded_update) = neighborhood_search
151143

152-
# Reset `cell_buffer` by moving all pointers to the beginning
153-
cell_buffer_indices .= 0
144+
# Empty each thread's list
145+
for i in eachindex(update_buffer)
146+
emptyat!(update_buffer, i)
147+
end
154148

155149
# Find all cells containing points that now belong to another cell
156-
mark_changed_cell!(neighborhood_search, cell_list, coords_fun,
157-
Val(threaded_update))
158-
159-
# Iterate over all marked cells and move the points into their new cells.
160-
for thread in 1:Threads.nthreads()
161-
# Only the entries `1:cell_buffer_indices[thread]` are initialized for `thread`.
162-
for i in 1:cell_buffer_indices[thread]
163-
cell_index = cell_buffer[i, thread]
150+
mark_changed_cell!(neighborhood_search, cell_list, coords_fun, Val(threaded_update))
151+
152+
# Iterate over all marked cells and move the points into their new cells
153+
for j in eachindex(update_buffer)
154+
for cell_index in update_buffer[j]
164155
points = cell_list[cell_index]
165156

166157
# Find all points whose coordinates do not match this cell
@@ -207,7 +198,7 @@ end
207198
# Otherwise, `@threaded` does not work here with Julia ARM on macOS.
208199
# See https://github.com/JuliaSIMD/Polyester.jl/issues/88.
209200
@inline function mark_changed_cell!(neighborhood_search, cell_index, coords_fun)
210-
(; cell_list, cell_buffer, cell_buffer_indices) = neighborhood_search
201+
(; cell_list, update_buffer) = neighborhood_search
211202

212203
for point in cell_list[cell_index]
213204
cell = cell_coords(coords_fun(point), neighborhood_search)
@@ -216,12 +207,8 @@ end
216207
# cell list to store cells inside `cell`.
217208
# These can be identical (see `DictionaryCellList`).
218209
if !is_correct_cell(cell_list, cell, cell_index)
219-
# Mark this cell and continue with the next one.
220-
#
221-
# `cell_buffer` is preallocated,
222-
# but only the entries 1:i are used for this thread.
223-
i = cell_buffer_indices[Threads.threadid()] += 1
224-
cell_buffer[i, Threads.threadid()] = cell_index
210+
# Mark this cell and continue with the next one
211+
pushat!(update_buffer, Threads.threadid(), cell_index)
225212
break
226213
end
227214
end

‎src/vector_of_vectors.jl‎

Lines changed: 96 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,96 @@
1+
# Data structure that behaves like a `Vector{Vector}`, but uses a contiguous memory layout.
2+
# Similar to `VectorOfVectors` of ArraysOfArrays.jl, but allows to resize the inner vectors.
3+
struct DynamicVectorOfVectors{T, ARRAY2D, ARRAY1D} <: AbstractVector{Array{T, 1}}
4+
backend::ARRAY2D # Array{T, 2}, where each column represents a vector
5+
length_::Base.RefValue{Int32} # Number of vectors
6+
lengths::ARRAY1D # Array{Int32, 1} storing the lengths of the vectors
7+
end
8+
9+
function DynamicVectorOfVectors{T}(; max_outer_length, max_inner_length) where {T}
10+
backend = Array{T, 2}(undef, max_inner_length, max_outer_length)
11+
length_ = Ref(zero(Int32))
12+
lengths = zeros(Int32, max_outer_length)
13+
14+
return DynamicVectorOfVectors{T, typeof(backend), typeof(lengths)}(backend, length_,
15+
lengths)
16+
end
17+
18+
@inline Base.size(vov::DynamicVectorOfVectors) = (vov.length_[],)
19+
20+
@inline function Base.getindex(vov::DynamicVectorOfVectors, i)
21+
(; backend, lengths) = vov
22+
23+
@boundscheck checkbounds(vov, i)
24+
25+
return view(backend, 1:lengths[i], i)
26+
end
27+
28+
@inline function Base.push!(vov::DynamicVectorOfVectors, vector::AbstractVector)
29+
(; backend, length_, lengths) = vov
30+
31+
# This data structure only supports one-based indexing
32+
Base.require_one_based_indexing(vector)
33+
34+
# Activate a new column of `backend`
35+
j = length_[] += 1
36+
lengths[j] = length(vector)
37+
38+
# Fill the new column
39+
for i in eachindex(vector)
40+
backend[i, j] = vector[i]
41+
end
42+
43+
return vov
44+
end
45+
46+
@inline function Base.push!(vov::DynamicVectorOfVectors, vector::AbstractVector, vectors...)
47+
push!(vov, vector)
48+
push!(vov, vectors...)
49+
end
50+
51+
# `push!(vov[i], value)`
52+
@inline function pushat!(vov::DynamicVectorOfVectors, i, value)
53+
(; backend, lengths) = vov
54+
55+
@boundscheck checkbounds(vov, i)
56+
57+
# Activate new entry in column `i`
58+
backend[lengths[i] += 1, i] = value
59+
60+
return vov
61+
end
62+
63+
# `deleteat!(vov[i], j)`
64+
@inline function deleteatat!(vov::DynamicVectorOfVectors, i, j)
65+
(; backend, lengths) = vov
66+
67+
# Outer bounds check
68+
@boundscheck checkbounds(vov, i)
69+
# Inner bounds check
70+
@boundscheck checkbounds(1:lengths[i], j)
71+
72+
# Replace value to delete by the last value in this column
73+
last_value = backend[lengths[i], i]
74+
backend[j, i] = last_value
75+
76+
# Remove the last value in this column
77+
lengths[i] -= 1
78+
79+
return vov
80+
end
81+
82+
@inline function Base.empty!(vov::DynamicVectorOfVectors)
83+
# Move all pointers to the beginning
84+
vov.lengths .= zero(Int32)
85+
vov.length_[] = zero(Int32)
86+
87+
return vov
88+
end
89+
90+
# `empty!(vov[i])`
91+
@inline function emptyat!(vov::DynamicVectorOfVectors, i)
92+
# Move length pointer to the beginning
93+
vov.lengths[i] = zero(Int32)
94+
95+
return vov
96+
end

‎test/nhs_grid.jl‎

Lines changed: 1 addition & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -91,9 +91,7 @@
9191

9292
# Create neighborhood search
9393
nhs1 = GridNeighborhoodSearch{3}(; search_radius, n_points)
94-
95-
coords_fun(i) = coordinates1[:, i]
96-
initialize_grid!(nhs1, coords_fun)
94+
initialize_grid!(nhs1, coordinates1)
9795

9896
# Get each neighbor for `point_position1`
9997
neighbors1 = sort(collect(PointNeighbors.eachneighbor(point_position1, nhs1)))

‎test/test_util.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
# All `using` calls are in this file, so that one can run any test file
22
# after running only this file.
3-
using Test: @test, @testset
3+
using Test: @test, @testset, @test_throws
44
using PointNeighbors
55

66
"""

‎test/unittest.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,7 @@
11
# Separate file that can be executed to only run unit tests.
22
# Include `test_util.jl` first.
33
@testset verbose=true "Unit Tests" begin
4+
include("vector_of_vectors.jl")
45
include("nhs_trivial.jl")
56
include("nhs_grid.jl")
67
include("neighborhood_search.jl")

‎test/vector_of_vectors.jl‎

Lines changed: 81 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,81 @@
1+
@testset verbose=true "`DynamicVectorOfVectors`" begin
2+
# Test different types by defining a function to convert to this type
3+
types = [Int32, Float64, i -> (i, i)]
4+
5+
@testset verbose=true "Eltype $(eltype(type(1)))" for type in types
6+
ELTYPE = typeof(type(1))
7+
vov_ref = Vector{Vector{ELTYPE}}()
8+
vov = PointNeighbors.DynamicVectorOfVectors{ELTYPE}(max_outer_length = 20,
9+
max_inner_length = 10)
10+
11+
# Test internal size
12+
@test size(vov.backend) == (10, 20)
13+
14+
function verify(vov, vov_ref)
15+
@test length(vov) == length(vov_ref)
16+
@test eachindex(vov) == eachindex(vov_ref)
17+
@test axes(vov) == axes(vov_ref)
18+
19+
@test_throws BoundsError vov[0]
20+
@test_throws BoundsError vov[length(vov) + 1]
21+
22+
for i in eachindex(vov_ref)
23+
@test vov[i] == vov_ref[i]
24+
end
25+
end
26+
27+
# Initial check
28+
verify(vov, vov_ref)
29+
30+
# First `push!`
31+
push!(vov_ref, type.([1, 2, 3]))
32+
push!(vov, type.([1, 2, 3]))
33+
34+
verify(vov, vov_ref)
35+
36+
# `push!` multiple items
37+
push!(vov_ref, type.([4]), type.([5, 6, 7, 8]))
38+
push!(vov, type.([4]), type.([5, 6, 7, 8]))
39+
40+
verify(vov, vov_ref)
41+
42+
# `push!` to an inner vector
43+
push!(vov_ref[1], type(12))
44+
PointNeighbors.pushat!(vov, 1, type(12))
45+
46+
verify(vov, vov_ref)
47+
48+
# Delete entry of inner vector. Note that this changes the order of the elements.
49+
deleteat!(vov_ref[3], 2)
50+
PointNeighbors.deleteatat!(vov, 3, 2)
51+
52+
@test vov_ref[3] == type.([5, 7, 8])
53+
@test vov[3] == type.([5, 8, 7])
54+
55+
# Delete second to last entry
56+
deleteat!(vov_ref[3], 2)
57+
PointNeighbors.deleteatat!(vov, 3, 2)
58+
59+
@test vov_ref[3] == type.([5, 8])
60+
@test vov[3] == type.([5, 7])
61+
62+
# Delete last entry
63+
deleteat!(vov_ref[3], 2)
64+
PointNeighbors.deleteatat!(vov, 3, 2)
65+
66+
# Now they are identical again
67+
verify(vov, vov_ref)
68+
69+
# Delete the remaining entry of this vector
70+
deleteat!(vov_ref[3], 1)
71+
PointNeighbors.deleteatat!(vov, 3, 1)
72+
73+
verify(vov, vov_ref)
74+
75+
# `empty!`
76+
empty!(vov_ref)
77+
empty!(vov)
78+
79+
verify(vov, vov_ref)
80+
end
81+
end

0 commit comments

Comments
 (0)