Efficient way of computing norm between all the rows of a matrix

New to Julia
,
6

This kind of data is not the one where cell lists shine the most, because it is very sparse (cell lists are ideal for more or less homogeneous distributions of points). Still, there it goes (I have adapted the naive version to use static arrays):

julia> include("./vinicius.jl")
cell list - serial: 
  5.258 s (13576 allocations: 68.23 MiB)
cell list - parallel: 
  2.096 s (77501 allocations: 131.05 MiB)
naive - serial: 
# still running, will update when finished...

Here is the code:

using CellListMap
using StaticArrays
using DelimitedFiles
using LinearAlgebra: norm
using BenchmarkTools

function count_close(p::AbstractVector{SVector{N,T}},r;parallel=true) where {N,T}
    p_length = limits(p) # length of set on each dimension
    # To prevent too many cells, set the cell side such that the maximum number of cells 
    # is at most equal to the number of points
    max_length = maximum(p_length.limits) 
    max_bins = max_length/r
    n = ceil(Int,min(max_bins,length(p)^(1/N))) 
    cell_side = max(r,max_length/n)
    # Prepare cell lists
    box = Box(p_length,cell_side)
    cl = CellList(p,box,parallel=parallel)
    r2 = r^2
    n_close = 0
    n_close = map_pairwise!(
        (x,y,i,j,d2,n_close) -> d2 < r2 ? n_close + 1 : n_close,
        n_close,
        box, cl, parallel=parallel
    )
    return n_close
end

data = readdlm("pontos.txt")
pontos = [ SVector{2,Float64}(data[i,1],data[i,2]) for i in axes(data,1) ]
pontos = pontos[1:100000]
r = 2e-2

println("cell list - serial: ")
@btime count_close($pontos,$r,parallel=false)

println("cell list - parallel: ")
@btime count_close($pontos,$r,parallel=true)

function C(pontos,r)
    contar = 0
    for r1 in 1:length(pontos)
        for r2 in r1+1:length(pontos)
            if norm(pontos[r2]-pontos[r1]) < r
                contar += 1 
            end
        end
    end
    return contar
end

println("naive - serial: ")
@btime C(pontos,r)

up: For 100_000 points this is the difference (for one million I think it will never end, had to kill it):

julia> include("./vinicius.jl")
cell list - serial: 
  74.476 ms (8778 allocations: 7.61 MiB)
cell list - parallel: 
  27.496 ms (40732 allocations: 17.19 MiB)
naive - serial: 
  8.912 s (1 allocation: 16 bytes)
15579644
6 Likes