# 求助！Cell List算法写出来没有两重for循环快

1. 两重for循环 for i = 1:NumberNod-1 for j = i+1:NumerNod
2. Cell List算法： 先将点放进cell中，然后只需要搜索一个邻近的8个cell即可

=============================================================================
#cell list algorithm

function FRNNMain(NumNod::Int64,Dim::Int64,A::Array{Float64,2},s::Float64)
r = s/sqrt(Dim)
PoiPair = Array{Tuple{Int64,Int64},1}(undef,1)
Cell = PoiCell(A,NumNod,r)
for key in keys(Cell)
vaule = Cell[key]
InnerCellPair!(vaule,PoiPair)
LoopNeighbor!(key,Cell,PoiPair,A,s)
end
return PoiPair
end
function PoiCell(A::Matrix{Float64},NumNod::Int64,r::Float64)
Cell = Dict()
for i = 1:1:NumNod
x = Int(floor(A[i,1]/r))
y = Int(floor(A[i,2]/r))
z = (x,y)
if z ∉ keys(Cell)
Cell[z] = [i]
else
push!(Cell[z],i)
end
end
return Cell
end
function GetNeighborCell(key::Tuple{Int64,Int64})
NeiCell = Array{Tuple}(undef,8,1)
NeiCell[1] = (key[1]+1,key[2]+1)
NeiCell[2] = (key[1]+1,key[2])
NeiCell[3] = (key[1]+1,key[2]-1)
NeiCell[4] = (key[1],key[2]+1)
NeiCell[5] = (key[1],key[2]-1)
NeiCell[6] = (key[1]-1,key[2])
NeiCell[7] = (key[1]-1,key[2]+1)
NeiCell[8] = (key[1]-1,key[2]-1)
return NeiCell
end
function InnerCellPair!(value::Array{Int64},pairs::Array{Tuple{Int64,Int64},1})
Len = length(value)
if Len >= 2
for i = 1:Len-1
for j = (i+1):Len
z = (value[i],value[j])
push!(pairs,z)
end
end
end
return pairs
end
function LoopNeighbor!(key::Tuple{Int64,Int64},Cell::Dict{Any,Any},PoiPair::Array{Tuple{Int64,Int64},1},A::Array{Float64,2},s::Float64)
value = Cell[key]
NeiCell = GetNeighborCell(key)
for i = 1:8
key1 = NeiCell[i]
if key1 ∈ keys(Cell)
value1 = Cell[NeiCell[i]]
FindPair!(value,value1,PoiPair,A,s)
end
end
return PoiPair
end
function FindPair!(value,value1,PoiPair::Array{Tuple{Int64,Int64},1},A::Array{Float64,2},s::Float64)
Len1 = length(value)
Len2 = length(value1)
for i = 1:Len1
for j = 1:Len2
P1 = value[i]
P2 = value1[j]
Point1 = A[P1,:]
Point2 = A[P2,:]
Dist = Distance(Point1,Point2)
if Dist <= s^2
end
end
end
return PoiPair
end
function Distance(Point1::Array{Float64,1},Point2::Array{Float64,1})
Dist = (Point1[1]-Point2[1])^2+(Point1[2]-Point2[2])^2
return Dist
end
pair1 = (P1,P2)
pair2 = (P2,P1)
if pair1 ∉ PoiPair
if pair2 ∉ PoiPair
push!(PoiPair,pair1)
end
end
return PoiPair
end
NumNod = 1000
Dim = 2
A = rand(NumNod,Dim)
s = 0.1;
a = FRNNMain(NumNod,Dim,A,s)


=============================================================================#forloop

NumNod = 1000
Dim = 2
A = rand(NumNod,Dim)
s = 0.1;
function forloop(A::Matrix{Float64},NumNod::Int64,s::Float64)
PoiPair = []
for i = 1:NumNod-1
for j = i+1:NumNod
Point1 = A[i,:]
Point2 = A[j,:]
Dist = (Point1[1]-Point2[1])^2+(Point1[2]-Point2[2])^2
if Dist < s
push!(PoiPair,(i,j))
end
end
end
return PoiPair
end
b = forloop(A,NumNod,s)


=============================================================================

@btime FRNNMain(NumNod,Dim,A,s)
142.194 ms (79608 allocations: 7.33 MiB)
@btime forloop(A,NumNod,s)
46.899 ms (1117806 allocations: 97.09 MiB)

julia
println(“hello world”)


r = s/sqrt(2)是cell的尺寸，这样cell里的点的距离就都小于s了

cell是什么？

@btime FRNNMain(NumNod,Dim,A,s)
142.194 ms (79608 allocations: 7.33 MiB)
@btime forloop(A,NumNod,s)
46.899 ms (1117806 allocations: 97.09 MiB)


cell就是小格子，相当于把二维平面划分成一个个的小格子，把点放到这些小格子里面。

Point1 = A[i,:]
Point2 = A[j,:]
Dist = (Point1[1]-Point2[1])^2+(Point1[2]-Point2[2])^2


    Len1 = length(value)
Len2 = length(value1)
for i = 1:Len1
for j = 1:Len2
P1 = value[i]
P2 = value1[j]


for p1 in value
for p2 in value1


@btime FRNNMain(NumNod,Dim,A,s)
201.197 ms (7702 allocations: 761.42 KiB)


# fied-radius nearest neighbor
#cell list algorithm
function FRNNMain(NumNod::Int64,Dim::Int64,A::Array{Float64,2},s::Float64)
r = s/sqrt(Dim)
PoiPair = Array{Tuple{Int64,Int64},1}(undef,1)
popfirst!(PoiPair)
Cell = PoiCell(A,NumNod,r)
for key in keys(Cell)
vaule = Cell[key]
InnerCellPair!(vaule,PoiPair)
LoopNeighbor!(key,Cell,PoiPair,A,s)
end
return PoiPair
end
function PoiCell(A::Matrix{Float64},NumNod::Int64,r::Float64)
Cell = Dict{Tuple{Int64,Int64},Array{Int64,1}}()
for i = 1:1:NumNod
x = fld(A[1,i],r)
y = fld(A[2,i],r)
z = (x,y)
if z ∉ keys(Cell)
Cell[z] = [i]
else
push!(Cell[z],i)
end
end
return Cell
end
function GetNeighborCell(key::Tuple{Int64,Int64})
NeiCell = Array{Tuple{Int64,Int64}}(undef,8,1)
NeiCell[1] = (key[1]+1,key[2]+1)
NeiCell[2] = (key[1]+1,key[2])
NeiCell[3] = (key[1]+1,key[2]-1)
NeiCell[4] = (key[1],key[2]+1)
NeiCell[5] = (key[1],key[2]-1)
NeiCell[6] = (key[1]-1,key[2])
NeiCell[7] = (key[1]-1,key[2]+1)
NeiCell[8] = (key[1]-1,key[2]-1)
return NeiCell
end
function InnerCellPair!(value::Array{Int64},pairs::Array{Tuple{Int64,Int64},1})
Len = length(value)::Int64
if Len >= 2
for i = 1:Len-1
for j = (i+1):Len
z = (value[i],value[j])
push!(pairs,z)
end
end
end
return pairs
end
function LoopNeighbor!(key::Tuple{Int64,Int64},Cell::Dict{Tuple{Int64,Int64},Array{Int64,1}},PoiPair::Array{Tuple{Int64,Int64},1},A::Array{Float64,2},s::Float64)
NeiCell = Array{Tuple{Int64,Int64}}(undef,8,1)
NeiCell = GetNeighborCell(key)
for i = 1:8
if NeiCell[i] ∈ keys(Cell)
FindPair!(Cell[key],Cell[NeiCell[i]],PoiPair,A,s)
end
end
return PoiPair
end
function FindPair!(value::Array{Int64},value1::Array{Int64},PoiPair::Array{Tuple{Int64,Int64},1},A::Array{Float64,2},s::Float64)
for P1 in value
for P2 in value1
#Dist = Distance(A[P1,:],A[P2,:])
Dist::Float64 = (A[1,P1]-A[1,P2])^2+(A[2,P1]-A[2,P2])^2
if Dist <= s^2
end
end
end
return PoiPair
end
function Distance(Point1::Array{Float64,1},Point2::Array{Float64,1})
Dist = (Point1[1]-Point2[1])^2+(Point1[2]-Point2[2])^2::Float64
return Dist
end
pair1 = (P1,P2)
pair2 = (P2,P1)
if pair1 ∉ PoiPair
if pair2 ∉ PoiPair
push!(PoiPair,pair1)
end
end
return PoiPair
end
NumNod = 1000
Dim = 2
A = rand(Dim,NumNod)
s = 0.1;
a = FRNNMain(NumNod,Dim,A,s)


@btime FRNNMain(NumNod,Dim,A,s)
209.789 ms (961 allocations: 670.45 KiB)


NumNod = 1000
Dim = 2
A = rand(Dim,NumNod)
s = 0.1;
function forloop(A::Matrix{Float64},NumNod::Int64,s::Float64)
PoiPair = Array{Tuple{Int64,Int64},1}(undef,1)
popfirst!(PoiPair)
for i = 1:NumNod-1
for j = i+1:NumNod
Dist = (A[1,i]-A[1,j])^2+(A[2,i]-A[2,j])^2
if Dist < s^2
push!(PoiPair,(i,j))
end
end
end
return PoiPair
end
b = forloop(A,NumNod,s)


@btime forloop(A,NumNod,s)
972.200 μs (14 allocations: 512.67 KiB)


@code_warntype FRNNMain(NumNod,Dim,A,s)
Variables
#self#::Core.Compiler.Const(FRNNMain, false)
NumNod::Int64
Dim::Int64
A::Array{Float64,2}
s::Float64
r::Float64
PoiPair::Array{Tuple{Int64,Int64},1}
Cell::Dict{Tuple{Int64,Int64},Array{Int64,1}}
@_9::Union{Nothing, Tuple{Tuple{Int64,Int64},Int64}}
key::Tuple{Int64,Int64}
vaule::Array{Int64,1}

Body::Array{Tuple{Int64,Int64},1}
1 ─ %1  = Main.sqrt(Dim)::Float64
│         (r = s / %1)
│   %3  = Core.apply_type(Main.Tuple, Main.Int64, Main.Int64)::Core.Compiler.Const(Tuple{Int64,Int64}, false)
│   %4  = Core.apply_type(Main.Array, %3, 1)::Core.Compiler.Const(Array{Tuple{Int64,Int64},1}, false)
│         (PoiPair = (%4)(Main.undef, 1))
│         Main.popfirst!(PoiPair)
│         (Cell = Main.PoiCell(A, NumNod, r))
│   %8  = Main.keys(Cell)::Base.KeySet{Tuple{Int64,Int64},Dict{Tuple{Int64,Int64},Array{Int64,1}}}
│         (@_9 = Base.iterate(%8))
│   %10 = (@_9 === nothing)::Bool
│   %11 = Base.not_int(%10)::Bool
└──       goto #4 if not %11
2 ┄ %13 = @_9::Tuple{Tuple{Int64,Int64},Int64}::Tuple{Tuple{Int64,Int64},Int64}
│         (key = Core.getfield(%13, 1))
│   %15 = Core.getfield(%13, 2)::Int64
│         (vaule = Base.getindex(Cell, key))
│         Main.InnerCellPair!(vaule, PoiPair)
│         Main.LoopNeighbor!(key, Cell, PoiPair, A, s)
│         (@_9 = Base.iterate(%8, %15))
│   %20 = (@_9 === nothing)::Bool
│   %21 = Base.not_int(%20)::Bool
└──       goto #4 if not %21
3 ─       goto #2
4 ┄       return PoiPair


julia 的 for 很给力。小样本量基本优势很大，大样本才能看出算法上的差距

### 改变数组大小 NumNod

julia> for NumNod in collect(1:10) * 1000
A = rand(NumNod,Dim)*100;
@time forloop(A,NumNod,s);
end
0.196111 seconds (999.02 k allocations: 91.462 MiB, 51.76% gc time)
0.249301 seconds (4.00 M allocations: 366.031 MiB, 19.79% gc time)
0.397443 seconds (9.00 M allocations: 823.709 MiB, 3.72% gc time)
1.482190 seconds (16.00 M allocations: 1.430 GiB, 2.91% gc time)
1.201030 seconds (25.00 M allocations: 2.235 GiB, 3.60% gc time)
1.675108 seconds (35.99 M allocations: 3.218 GiB, 3.59% gc time)
2.311076 seconds (48.99 M allocations: 4.380 GiB, 3.65% gc time)
4.310118 seconds (63.99 M allocations: 5.721 GiB, 3.03% gc time)

julia> for NumNod in collect(1:10) * 10_000
A = rand(NumNod,Dim)*100;
@time forloop(A,NumNod,s);
end
4.261473 seconds (99.99 M allocations: 8.940 GiB, 4.27% gc time)
18.061562 seconds (399.99 M allocations: 35.761 GiB, 4.63% gc time)
37.787276 seconds (899.98 M allocations: 80.464 GiB, 4.09% gc time)
114.836089 seconds (1.60 G allocations: 143.049 GiB, 3.38% gc time)
^C

julia> for NumNod in collect(1:10) * 1000
A = rand(NumNod,Dim)*100;
@time FRNNMain(NumNod,Dim,A,s);
end
0.037698 seconds (55.18 k allocations: 2.433 MiB)
0.054817 seconds (92.95 k allocations: 3.535 MiB)
0.005024 seconds (76.64 k allocations: 2.737 MiB)
0.010408 seconds (103.70 k allocations: 3.610 MiB)
0.022965 seconds (134.55 k allocations: 4.594 MiB)
0.022106 seconds (167.60 k allocations: 5.660 MiB)
0.040841 seconds (213.71 k allocations: 7.119 MiB, 33.68% gc time)
0.025413 seconds (270.34 k allocations: 8.902 MiB)
0.042115 seconds (348.90 k allocations: 11.345 MiB)
0.038399 seconds (475.34 k allocations: 15.266 MiB)

julia> for NumNod in collect(1:10) * 10_000
A = rand(NumNod,Dim)*100;
@time FRNNMain(NumNod,Dim,A,s);
end
0.039367 seconds (441.21 k allocations: 14.224 MiB, 25.38% gc time)
0.050425 seconds (539.30 k allocations: 18.547 MiB)
0.089731 seconds (961.16 k allocations: 32.120 MiB, 13.99% gc time)
0.184787 seconds (1.80 M allocations: 58.545 MiB, 7.14% gc time)
0.177084 seconds (1.32 M allocations: 47.563 MiB, 20.97% gc time)
0.241691 seconds (1.57 M allocations: 56.436 MiB, 25.85% gc time)
0.224645 seconds (1.87 M allocations: 66.499 MiB, 5.91% gc time)
0.279534 seconds (2.18 M allocations: 77.317 MiB, 5.94% gc time)
0.348275 seconds (2.52 M allocations: 88.803 MiB, 6.01% gc time)
0.521723 seconds (2.91 M allocations: 101.983 MiB, 7.30% gc time)

julia> for NumNod in collect(1:10) * 100_000
A = rand(NumNod,Dim)*100;
@time FRNNMain(NumNod,Dim,A,s);
end
0.535916 seconds (2.91 M allocations: 102.119 MiB, 9.12% gc time)
5.450482 seconds (6.25 M allocations: 225.615 MiB, 4.79% gc time)
28.951260 seconds (11.60 M allocations: 416.061 MiB, 1.47% gc time)
87.724544 seconds (13.12 M allocations: 504.621 MiB, 0.63% gc time)
336.967927 seconds (17.41 M allocations: 680.904 MiB, 0.22% gc time)
^C


### 改变分散程度

julia> NumNod = 5000
5000

julia> for scale in 1:10
A = rand(NumNod,Dim)*scale*100;
@time forloop(A,NumNod,s);
end
1.398829 seconds (25.00 M allocations: 2.235 GiB, 3.46% gc time)
1.336941 seconds (25.00 M allocations: 2.235 GiB, 3.55% gc time)
1.528647 seconds (25.00 M allocations: 2.235 GiB, 3.32% gc time)
1.312679 seconds (25.00 M allocations: 2.235 GiB, 3.64% gc time)
1.321496 seconds (25.00 M allocations: 2.235 GiB, 3.57% gc time)
1.362860 seconds (25.00 M allocations: 2.235 GiB, 3.53% gc time)
1.601800 seconds (25.00 M allocations: 2.235 GiB, 3.53% gc time)
1.664082 seconds (25.00 M allocations: 2.235 GiB, 3.60% gc time)
2.062444 seconds (25.00 M allocations: 2.235 GiB, 2.87% gc time)
2.135823 seconds (25.00 M allocations: 2.235 GiB, 2.97% gc time)


julia> NumNod = 200_000
200000

julia> for scale in 1:10
A = rand(NumNod,Dim)*scale*100;
@time FRNNMain(NumNod,Dim,A,s);
end
7.165528 seconds (6.25 M allocations: 225.755 MiB, 4.89% gc time)
1.157459 seconds (6.20 M allocations: 210.068 MiB, 12.50% gc time)
1.126173 seconds (6.16 M allocations: 206.322 MiB, 31.87% gc time)
0.874867 seconds (6.18 M allocations: 205.859 MiB, 16.87% gc time)
1.068134 seconds (6.16 M allocations: 205.089 MiB, 25.65% gc time)
1.082324 seconds (6.17 M allocations: 205.182 MiB, 23.96% gc time)
1.161492 seconds (6.18 M allocations: 205.316 MiB, 24.34% gc time)
1.142114 seconds (6.18 M allocations: 205.177 MiB, 23.17% gc time)
0.886259 seconds (6.17 M allocations: 204.759 MiB, 20.16% gc time)
1.234273 seconds (6.18 M allocations: 204.897 MiB, 28.85% gc time)

julia> NumNod = 300_000
300000

julia> for scale in 1:10
A = rand(NumNod,Dim)*scale*100;
@time FRNNMain(NumNod,Dim,A,s);
end
29.940059 seconds (11.67 M allocations: 418.380 MiB, 1.40% gc time)
3.063644 seconds (12.31 M allocations: 407.787 MiB, 7.55% gc time)
2.303256 seconds (12.58 M allocations: 409.451 MiB, 16.61% gc time)
1.653237 seconds (12.58 M allocations: 407.519 MiB, 18.47% gc time)
1.580190 seconds (12.66 M allocations: 409.065 MiB, 13.93% gc time)
4.677850 seconds (12.64 M allocations: 407.666 MiB, 18.80% gc time)
2.241075 seconds (12.70 M allocations: 409.340 MiB, 19.28% gc time)
3.397634 seconds (12.69 M allocations: 408.701 MiB, 8.65% gc time)
2.453740 seconds (12.70 M allocations: 408.790 MiB, 21.53% gc time)
1.601407 seconds (12.69 M allocations: 408.511 MiB, 21.20% gc time)

julia>


   NNN    index deltas:    (-1,  2) ( 0,  2) ( 1,  2)
NNNNN           (-2,  1) (-1,  1) ( 0,  1) ( 1,  1) ( 2,  1)
NNANN           (-2,  0) (-1,  0) ( 0,  0) ( 1,  0) ( 2,  0)
NNNNN           (-2, -1) (-1, -1) ( 0, -1) ( 1, -1) ( 2, -1)
NNN                     (-1, -2) ( 0, -2) ( 1, -2)


sf 的评论里提到过 r 的选取问题。

NumNod = 1000
Dim = 2
s = 0.1
println("Cell List")
for i = collect(1:10)*NumNod
A = rand(Dim,i)*100
@time FRNNMain(NumNod,Dim,A,s)
end
println("For Loop")
for i = collect(1:10)*NumNod
A = rand(Dim,i)*100
@time forloop(A,NumNod,s)
end


Cell List
0.001097 seconds (3.02 k allocations: 1010.063 KiB)
0.000870 seconds (3.02 k allocations: 1010.141 KiB)
0.000847 seconds (3.02 k allocations: 1010.063 KiB)
0.000886 seconds (3.02 k allocations: 1010.141 KiB)
0.000884 seconds (3.02 k allocations: 1010.141 KiB)
0.001272 seconds (3.02 k allocations: 1010.281 KiB)
0.000889 seconds (3.02 k allocations: 1010.063 KiB)
0.001526 seconds (3.02 k allocations: 1010.063 KiB)
0.000702 seconds (3.02 k allocations: 1009.313 KiB)
0.000887 seconds (3.02 k allocations: 1010.063 KiB)
For Loop
0.001759 seconds (1.00 k allocations: 3.961 MiB)
0.005201 seconds (1.00 k allocations: 3.961 MiB)
0.003741 seconds (1.00 k allocations: 3.961 MiB)
0.003309 seconds (1.00 k allocations: 3.961 MiB)
0.004977 seconds (1.00 k allocations: 3.961 MiB)
0.003765 seconds (1.00 k allocations: 3.961 MiB)
0.004728 seconds (1.00 k allocations: 3.961 MiB)
0.003427 seconds (1.00 k allocations: 3.961 MiB)
0.003933 seconds (1.00 k allocations: 3.961 MiB)
0.016121 seconds (1.00 k allocations: 3.961 MiB, 75.12% gc time)


NumNod = 1000
Dim = 2
s = 0.1*100
println("Cell List")
for i = collect(1:10)*NumNod
A = rand(Dim,i)*100
@time FRNNMain(NumNod,Dim,A,s)
end
println("For Loop")
for i = collect(1:10)*NumNod
A = rand(Dim,i)*100
@time forloop(A,NumNod,s)
end


Cell List
0.150363 seconds (978 allocations: 751.547 KiB)
0.142117 seconds (965 allocations: 748.063 KiB)
0.139737 seconds (963 allocations: 747.750 KiB)
0.150158 seconds (956 allocations: 746.719 KiB)
0.146481 seconds (979 allocations: 752.359 KiB)
0.142613 seconds (958 allocations: 747.078 KiB)
0.138542 seconds (974 allocations: 752.422 KiB)
0.137803 seconds (966 allocations: 747.016 KiB)
0.141860 seconds (977 allocations: 750.859 KiB)
0.149367 seconds (972 allocations: 750.219 KiB)
For Loop
0.002087 seconds (1.01 k allocations: 4.461 MiB)
0.002893 seconds (1.01 k allocations: 4.461 MiB)
0.004925 seconds (1.01 k allocations: 4.461 MiB)
0.003822 seconds (1.01 k allocations: 4.461 MiB)
0.015138 seconds (1.01 k allocations: 4.461 MiB, 76.00% gc time)
0.001978 seconds (1.01 k allocations: 4.461 MiB)
0.004069 seconds (1.01 k allocations: 4.461 MiB)
0.003221 seconds (1.01 k allocations: 4.461 MiB)
0.003352 seconds (1.01 k allocations: 4.461 MiB)
0.003903 seconds (1.01 k allocations: 4.461 MiB)


FRNNMain的速度还是远慢于forloop

forloop里面不管s是多少，每个点都要搜一遍，所以肯定比不搜来得慢鸭。

（困在家里的俺只有轻薄本，i5-8250U的CPU10000个点完全跑不动）