# Any suggestions for starting doing this kind of neighbour list?

**URL:** <https://discourse.julialang.org/t/any-suggestions-for-starting-doing-this-kind-of-neighbour-list/37436>\
**Category:** Numerics\
**Created:** [April 12, 2020, 6:40am UTC](https://discourse.julialang.org/t/any-suggestions-for-starting-doing-this-kind-of-neighbour-list/37436 "2020-04-12T06:40:14Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ahmed\_Salih](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_salih/32/206579_2.png) [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Post date:** [April 12, 2020, 6:40am UTC](https://discourse.julialang.org/t/any-suggestions-for-starting-doing-this-kind-of-neighbour-list/37436/1 "2020-04-12T06:40:14Z")

</div>

Hello!

I am trying to make a neighbour list like this, where one has N number of particles. Overlayed on top of these particles is a square grid with equal sizes. The goal is then output a structure/array/something else, which would hold information like, Cell[13] = [7 2 5 8], which means that Cell number 13 would store the particle id’s 7 2 5 and 8 (i.e. the four dots in Cell 13 - see figure).

I know this kind of neighbour list is a bit primitive, but it is a step up from brute force… I am struggling to figure out where to start, since if I want to figure out where a particle is for example for this 2D case, take x as horizontal axis and y as vertical;

1. Generate grid (which I am a bit unsure about how to do efficiently)
2. Make four “if” statements to check whether a particle is inside each side of a square in the grid.
3. Then place the particle id in corresponding cell

In theory I don’t have a difficult time understanding what I want to do, the difficulty lies in how to go about it. I would like to get some feedback/suggestions, maybe I should do it in a completely other way, before I dive to deep down the rabbit hole.

Kind regards

 ![image](https://global.discourse-cdn.com/julialang/original/3X/5/3/533f6f052bf206aecea20a550bd282f18c715746.png)  
**Turbulence and Viscous Mixing using Smoothed Particle Hydrodynamics by Martin Robinson**

---

<div class="post-metadata">

**Author:** ![improbable22](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/improbable22/32/5464_2.png) [@improbable22](https://discourse.julialang.org/u/improbable22)\
**Post date:** [April 12, 2020, 6:58am UTC](https://discourse.julialang.org/t/any-suggestions-for-starting-doing-this-kind-of-neighbour-list/37436/2 "2020-04-12T06:58:05Z")

</div>

Here is one quick go at something like this…

```julia
julia> grid = [Set{Int}() for x=1:5, y=1:5];

julia> data = [(x=1+4*rand(), y=1+4*rand(), n=n) for n in 1:30];

julia> insert(grid) = p -> push!(grid[trunc(Int, p.x), trunc(Int, p.y)], p.n);

julia> foreach(insert(grid), data)

julia> grid
5×5 Array{Set{Int64},2}:
 Set([23]) Set([]) Set([26, 30]) Set([4, 27, 28, 15]) Set([])
 Set([24, 16, 20]) Set([9, 14, 10, 6, 1]) Set([3, 17, 8]) Set([5, 21]) Set([])
 Set([19, 22]) Set([7, 25]) Set([18]) Set([2, 29]) Set([])
 Set([13]) Set([11, 12]) Set([]) Set([]) Set([])
 Set([]) Set([]) Set([]) Set([]) Set([])

```

---

<div class="post-metadata">

**Author:** ![Ahmed\_Salih](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_salih/32/206579_2.png) [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Post date:** [April 12, 2020, 7:21am UTC](https://discourse.julialang.org/t/any-suggestions-for-starting-doing-this-kind-of-neighbour-list/37436/3 "2020-04-12T07:21:09Z")

</div>

> [@improbable22](#):
>
> foreach(insert(grid), data)

I will be honest, I am very impressed with this, just tested for 1 million particles and it takes 0.1 seconds. I have some further questions, but I am going to work a bit with it and see how far I can get in understanding first.

Thank you very much!

Kind regards

---

<div class="post-metadata">

**Author:** ![Ahmed\_Salih](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_salih/32/206579_2.png) [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Post date:** [April 12, 2020, 8:11am UTC](https://discourse.julialang.org/t/any-suggestions-for-starting-doing-this-kind-of-neighbour-list/37436/4 "2020-04-12T08:11:18Z")

</div>

So I have a few follow up questions now, I edited the code slightly. Note, for some reason using Atom the plots do not appear / bug occurs, but it worked for me if I just copy pasted all the code to a terminal in v1.2.

```julia
using Plots

#Define equal sized square grids
gridc = [Set{Int}() for x=1:5, y=1:5];
#Generate arbitrary 2D data
data = [(x=1+4*rand(), y=1+4*rand(), n=n) for n in 1:30];

# Define a function "insert", but I have a hard time understanding what it does
# Could you try explaining this step by step?
# I assume p is basically "data", and then we check for each square what "data"
# lies in it
insert(gridc) = p -> push!(gridc[trunc(Int, p.x), trunc(Int, p.y)], p.n);

# "foreach" syntax is a bit new for me, I get confused since we only have one
# variable "data"?
foreach(insert(gridc), data)

# We will continue onward using the GR backend
gr()
x = [x[1] for x in data]
y = [y[2] for y in data]
plot(x, y, seriestype = :scatter, title = "My Scatter Plot",aspect_ratio=:equal,series_annotations = text.(1:length(x), :bottom))

```

I hope you could explain to me, step by step, in the function, “insert(gridc)”?

Also how do I plot the “physical” limits of the squares in the grid? I can that the input in gridc right now is a 5 by 5 square, from x = 1 to x = 5, so it seems like it is a set of 4x4 squares?

Sorry for my newbie questions, I just would never have thought to do it like this.

Kind regards

---

<div class="post-metadata">

**Author:** ![improbable22](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/improbable22/32/5464_2.png) [@improbable22](https://discourse.julialang.org/u/improbable22)\
**Post date:** [April 12, 2020, 8:49am UTC](https://discourse.julialang.org/t/any-suggestions-for-starting-doing-this-kind-of-neighbour-list/37436/5 "2020-04-12T08:49:58Z")

</div>

Sure. `foreach(f, A)` is like `map(f,A)` which is a lot like `f.(A)`, except that it doesn’t save what’s returned by `f` anywhere. And the function I handed it is `f(p::NamedTuple) = push!(gridc[…], p.n)`, it doesn’t return anything useful, but it does alter the grid. In a less functional mood you could equally well write

```julia
for p in data
    i,j = trunc(Int, p.x), trunc(Int, p.y) # turn real x,y into integer i,j
    s = grid[i,j] # get the Set for this cell
    push!(s, p.n) # store the point’s number in that
end

```

You could also do `for (n,p) in enumerate(data)`, then `data[n]` is the point in question, rather than expecting `n` to be a field of each point.

Using `trunc(Int, x)` means that the all `3 <= x < 4` belong to cell `i==3`. My example is a 5x5 grid, for numbers `1<=x<=6`, although by mistake I generated data always `x<5` which is why the last column is empty (and last row). The function `searchsortedfirst(fences, x)` might be a better choice in general.

The advantage of this grid of Sets is that, from the `(i,j)` of one point, it’s easy to find the neighbouring points. In the picture you posted, figuring out all the neighbours of cell 13 is harder work. (You can still access e.g. `grid[3]` if you wish, because arrays allow linear indexing.)

---

<div class="post-metadata">

**Author:** ![Ahmed\_Salih](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_salih/32/206579_2.png) [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Post date:** [April 12, 2020, 9:24am UTC](https://discourse.julialang.org/t/any-suggestions-for-starting-doing-this-kind-of-neighbour-list/37436/6 "2020-04-12T09:24:35Z")

</div>

Thank you very much for your explanation. This makes a lot more sense now. So now I am trying to apply this functionality to an actual problem of mine seen below;

![image](https://global.discourse-cdn.com/julialang/original/3X/6/9/6924b1bee7cb52dc990277adc8fac3255b567399.png)

Note that originally this data was scattered around (0,0), so the range was from -1 to 1, but the insert function would not work, so I had to offset the data by 2 on both x and y. This is fine since I am just doing neighbour search anyways.

I have checked again and it seems like the code always produces a row and column too much of gridc, so I simply added “-1” as such;

```julia
gridc = [Set{Int}() for x=1:1:3-1, y=1:1:3-1];

```

And gridc return 2x2 now, which I expect, since I basically for the above example define squares from 1 to 2 and 2 to 3 in both axis, so 2x2. I also confirm it works by noting that each square has same number of particles, so that is great!

So now comes my final difficulty for now.

I want to use this neighbour search for an SPH (Smoothed Particle Hydrodynamics), which means that I would prefer to use much smaller squares, for example something like;

```julia
gridc = [Set{Int}() for x=1:0.05:3-1, y=1:0.05:3-1];

```

But doing this does not produce the correct results anymore, which I believe is due to “trunc” not working as expected, i.e. as you described before.

If I try to scale my data with the inverse, i.e. 1/0.05 = 20, then I start getting errors like;

```julia
BoundsError: attempt to access 21×21 Array{Set{Int64},2} at index [20, 38]
getindex at array.jl:729 [inlined]
#293 at NL.jl:16 [inlined]
foreach(::getfield(Main, Symbol("##293#294")){Array{Set{Int64},2}}, ::Array{NamedTuple{(:x, :y, :n),Tuple{Float64,Float64,Int64}},1}) at abstractarray.jl:1920
top-level scope at NL.jl:20

```

So I am a bit at a loss at what to do. Do you think it is possible to keep using your approach, since it is really neat?

The code I am using currently (without scaling by 20 on x and y) is here:

```julia
using Plots
using DelimitedFiles

#Define equal sized square grids
gridc = [Set{Int}() for x=1:0.05:3-1, y=1:0.05:3-1];

#Generate arbitrary 2D data
dataini = readdlm("CircleGridUniform.txt", '\t', Float64, '\n')
data = [(x=2+dataini[n,1], y=2+dataini[n,2], n=n) for n in 1:size(dataini)[1]];
#data = [(x=1+4*rand(), y=1+4*rand(), n=n) for n in 1:30];

# Define a function "insert", but I have a hard time understanding what it does
# Could you try explaining this step by step?
# I assume p is basically "data", and then we check for each square what "data"
# lies in it
insert(gridc) = p -> push!(gridc[trunc(Int, p.x), trunc(Int, p.y)], p.n);

# "foreach" syntax is a bit new for me, I get confused since we only have one
# variable "data"?
foreach(insert(gridc), data)

# We will continue onward using the GR backend
gr()
x = [x[1] for x in data]
y = [y[2] for y in data]
#series_annotations = text.(1:length(x), :bottom)
plot(x, y, seriestype = :scatter, title = "My Scatter Plot",aspect_ratio=:equal)

```

And the datafile is too big to copy and paste here, but using;

```julia
data = [(x=1+4*rand(), y=1+4*rand(), n=n) for n in 1:10000];

```

One can see the same problem I experience, with gridc being filled incorrectly due to the trunc I suppose.

Kind regards

---

<div class="post-metadata">

**Author:** ![Ahmed\_Salih](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_salih/32/206579_2.png) [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Post date:** [April 12, 2020, 12:00pm UTC](https://discourse.julialang.org/t/any-suggestions-for-starting-doing-this-kind-of-neighbour-list/37436/8 "2020-04-12T12:00:24Z")

</div>

So finally I got it to work using “trunc”. The full code is down below;

```julia
using Plots
using DelimitedFiles

# REMEMBER TO CHANGE DIGITS AT INSERT FUNCTION TO NUMBER OF DIGITS OF sp
sp = 0.5;
sp1 = 1/sp;

#Define equal sized square grids
gridc = [Set{Int}() for x=(1:sp:3), y=(1:sp:3)];

#Data from:
#https://www.researchgate.net/publication/256688486_Particle_packing_algorithm_for_SPH_schemes
dataini = readdlm("CircleGridUniform.txt", '\t', Float64, '\n')
#Here is the scaling of data based on sp
data = [(x=(2+dataini[n,1])*sp1, y=(2+dataini[n,2])*sp1, n=n) for n in 1:size(dataini)[1]];
#data = [(x=1+4*rand(), y=1+4*rand(), n=n) for n in 1:10000];

#insert values in gridc based on indices given as below
insert(gridc) = p -> push!(gridc[floor(Int,trunc(p.x,digits=1)-sp1+1), floor(Int,trunc(p.y,digits=1)-sp1+1)], p.n);

foreach(insert(gridc), data)

# We will continue onward using the GR backend
gr()
x = [x[1] for x in data]
y = [y[2] for y in data]
#series_annotations = text.(1:length(x), :bottom)
#p = plot(x[collect(gridc[2,2])], y[collect(gridc[2,2])], seriestype = :scatter, title = "My Scatter Plot",aspect_ratio=:equal)

for i = 1:size(gridc)[1]
    for j = 1:size(gridc)[1]
        p = plot!(x[collect(gridc[i,j])],y[collect(gridc[i,j])],seriestype = :scatter,show = true,marker=([:hex]), markersize = 2.25,aspect_ratio=:equal)
        display(p)
    end
end
#current()

```

Again this is not the most fool-proof way, some tinkering might have to be done, but for my purpose and current level this is neat. It produces this plot:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/e/a/eab639b7ec618ec8c0cf90ad57a4c645a7d75a3d.png)

Where the different colors denotes connection to different cells. One should note that in the solution I posted the last row and colum of gridc will always be irrelevant, while the rest will be the actual data. Confirm this by using sp = 0.5.

An optimization to this would therefore be to remove this (couldn’t figure out how) and another optimization would be to make a function, which could look at “0.5” and figure out by it self that it has one digit after comma. I couldn’t figure that out either.

@improbable22 I have marked your initial answer as solution - thank you for helping out!

Kind regards

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [April 12, 2020, 1:57pm UTC](https://discourse.julialang.org/t/any-suggestions-for-starting-doing-this-kind-of-neighbour-list/37436/9 "2020-04-12T13:57:18Z")

</div>

Just pointing out that equal-sized cells will cause you to have to compare (significantly) more points than needed. You could try out for example `KDTree` from [GitHub - KristofferC/NearestNeighbors.jl: High performance nearest neighbor data structures and algorithms for Julia.](https://github.com/KristofferC/NearestNeighbors.jl) which splits the cells based on the actual points like:

 ![](https://global.discourse-cdn.com/julialang/original/3X/c/8/c8a9b2b3b9790156f3d98d52771eb18dade8d7bf.png)

---

<div class="post-metadata">

**Author:** ![Ahmed\_Salih](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_salih/32/206579_2.png) [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Post date:** [April 12, 2020, 2:36pm UTC](https://discourse.julialang.org/t/any-suggestions-for-starting-doing-this-kind-of-neighbour-list/37436/10 "2020-04-12T14:36:35Z")

</div>

That is also fairly interesting! The reason that one goes for uniform cells in SPH usually afaik, is that most of the time the particles are clumped “nicely”, and not distributed more sparsely as in your picture. Sometimes on purpose the cells are also left bigger than needed, so one can “get away with” only making a neighboursearch, say every 10 time steps.

Since it is a physical simulation it is also based on the Euclidean distance between particles, and I don’t know if the KD tree follows the same metric?

Never the less, thanks for showing me 🙂

Kind regards

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [April 12, 2020, 7:34pm UTC](https://discourse.julialang.org/t/any-suggestions-for-starting-doing-this-kind-of-neighbour-list/37436/11 "2020-04-12T19:34:33Z")

</div>

> [@Ahmed\_Salih](#):
>
> Since it is a physical simulation it is also based on the Euclidean distance between particles, and I don’t know if the KD tree follows the same metric?

KD Trees work with any Minkowski distance (and thus Euclidean distance).

---

<div class="post-metadata">

**Author:** ![Ahmed\_Salih](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_salih/32/206579_2.png) [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Post date:** [April 16, 2020, 5:25pm UTC](https://discourse.julialang.org/t/any-suggestions-for-starting-doing-this-kind-of-neighbour-list/37436/12 "2020-04-16T17:25:43Z")

</div>

So now I finally got to a point, where I can try to implement some of your neighbour search techniques. Currently I am applying “rangesearch” from your NearestNeighbors package, as such:

```julia
pg_arr = reshape(collect(Iterators.flatten(pg)),(3,nTot))
    balltree = BallTree(pg_arr,reorder = false)

    for t = 1:100
        if mod(t,5) == 0
            pg_arr .= reshape(collect(Iterators.flatten(pg)),(3,nTot))
            balltree = BallTree(pg_arr,reorder = false)
        end
        PackStep(pg,pg_tmp,u,u_tmp,nCoords,nTot,balltree)
    end

```

pg is an Array of tuple, that is why I have to do that “ugly” line to turn it into a typical array which the package can work with. I can see that this does not affect performance too much, so let us avoid that for now.

Currently I am remaking the balltree every 5th time step based on updated point values as seen in the loop. Inside of “PackStep” I then do as such:

```julia
function PackStep(pg,pg_tmp,u,u_tmp,nCoords,nTot,balltree)
    #idxs = Array{Int64,1};
    @fastmath @inbounds @threads for i = 1:nCoords
        Wgx = 0.0;
        Wgz = 0.0;
        p_i = pg[i];

        idxs = inrange(balltree, collect(p_i), R, false)

        filter!(x->x≠i,idxs)

```

I remind that p\_i is a tuple, that is why I use collect. I don’t care about the ordering so that is why it is “false”. I have to use a filter operation to remove the idxs value of the current particle I am working with - this is also fine.

Now the problem lies in, I cannot figure out if I can allocate this in a smarter way, since I am experiencing imo a big number of allocations:

```julia
0.300602 seconds (1.89 M allocations: 192.338 MiB)

```

Where without using this neighbour search technique, I would be around 300k allocations and 50 MiB.

I know I am not showing much code, but maybe you could pinpoint if I am doing something terribly wrong?

I see that I actually get a 25% speed up using this, so I really want to stick with it, just afraid that memory might get out of hand quicker this way.

Kind regards

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [April 16, 2020, 7:58pm UTC](https://discourse.julialang.org/t/any-suggestions-for-starting-doing-this-kind-of-neighbour-list/37436/13 "2020-04-16T19:58:23Z")

</div>

> [@Ahmed\_Salih](#):
>
> I see that I actually get a 25% speed up using this, so I really want to stick with it, just afraid that memory might get out of hand quicker this way.

Make sure you try with `KDTree` (instead of `BallTree`) as well because it has been a bit more optimized.

> [@Ahmed\_Salih](#):
>
> I cannot figure out if I can allocate this in a smarter way, since I am experiencing imo a big number of allocations:

Allocations by themselves are really not a problem. (Note that the number presented is the accumulated allocations over the time measured, not peak memory usage). The Julia GC is fast and is optimized for small allocations. What you should be interested in in the vast majority of cases is time (of course removing allocations can also improve runtime).

I would recommend you to store your points as a `Vector` of `StaticArray`s from [GitHub - JuliaArrays/StaticArrays.jl: Statically sized arrays for Julia](https://github.com/JuliaArrays/StaticArrays.jl) (instead of tuples). You can then create a tree by just passing that vector in and also you can call `inrange` with multiple points, for example:

```julia
julia> using NearestNeighbors

julia> using StaticArrays

julia> points = [rand(SVector{3}) for i in 1:10]
10-element Array{SArray{Tuple{3},Float64,1,3},1}:
 [0.520554391484928, 0.14816182639676279, 0.7495092644743351]
 [0.4921570494179721, 0.9936676737375236, 0.876053293047401]
...

julia> tree = KDTree(points)
KDTree{SArray{Tuple{3},Float64,1,3},Euclidean,Float64}
  Number of points: 10
  Dimensions: 3
  Metric: Euclidean(0.0)
  Reordered: true

julia> query_points = [rand(SVector{3}) for i in 1:3]
3-element Array{SArray{Tuple{3},Float64,1,3},1}:
 [0.2330633929342978, 0.3178020382035527, 0.4728928083723807]
 [0.15978776776689108, 0.11545355330021723, 0.1913731910261416]
 [0.7876520459805436, 0.727945529856302, 0.6189995574567859]

julia> inrange(tree, query_points, 0.5)
3-element Array{Array{Int64,1},1}:
 [1, 5, 9, 10]
 []
 [2, 4, 6, 7, 10]

```

---

<div class="post-metadata">

**Author:** ![Ahmed\_Salih](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ahmed_salih/32/206579_2.png) [@Ahmed\_Salih](https://discourse.julialang.org/u/Ahmed_Salih)\
**Post date:** [April 16, 2020, 9:59pm UTC](https://discourse.julialang.org/t/any-suggestions-for-starting-doing-this-kind-of-neighbour-list/37436/14 "2020-04-16T21:59:37Z")

</div>

Thanks for your reply!

I tested kNN but did not really see any improvement for this case.

Thanks for your explanation about allocations. I found that if I allocated all of the treesearch initially I could cut another 0.1 seconds off, which means that using neighbour search gives about 50% performance boost now. Also I managed by preallocating, i.e. doing calculation all at once as you showed to get the following:

```julia
0.216613 seconds (405.92 k allocations: 35.082 MiB)

```

Which is much more reasonable from before. I can see that I am able to cut even one more, 0.1 seconds off, if I can remove this filtering for loop:

```julia
filter!(x->x≠i,idxs[i])

```

Where idxs[1] is Array{Int64,1}. This is my only problem with the rangesearch right now, that it also return the particle it self. If you happen to have a one liner which can take an array like:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/9/f/9fb2e43964c6efad854568149c1bdb59bb1e6101.png)

And remove the yellow elements, i.e. index of the array - then I would really appreciate it.

EDIT: I will also be looking at StaticArrays.jl it seems, since you showed that it works really well with the neighbour search algorithms.

Kind regards