Any suggestions for starting doing this kind of neighbour list?

[Ahmed_Salih]

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

Turbulence and Viscous Mixing using Smoothed Particle Hydrodynamics by Martin Robinson

[improbable22]

Here is one quick go at something like this…

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([])

[Ahmed_Salih]

@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

[Ahmed_Salih]

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.

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

[improbable22]

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

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.)

[Ahmed_Salih]

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;

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;

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;

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;

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:

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;

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

[Ahmed_Salih]

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

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:

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

[kristoffer.carlsson]

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. which splits the cells based on the actual points like:

[Ahmed_Salih]

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

[kristoffer.carlsson]

@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).

[Ahmed_Salih]

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:

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:

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:

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

[kristoffer.carlsson]

@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 StaticArrays from GitHub - JuliaArrays/StaticArrays.jl: Statically sized arrays for Julia (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> 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]

[Ahmed_Salih]

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:

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:

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:

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