Just for the sake of didactics:
The “natural” way to solve that problem is probably:
julia> function meandist(n)
d = 0.
for i in 1:n
x = rand(2)
y = rand(2)
d += sqrt((x[1]-y[1])^2+(x[2]-y[2])^2)
end
return d/n
end
julia> @btime meandist(1000)
88.329 μs (2000 allocations: 187.50 KiB)
0.5351316591605562
This is quite slow and allocates a lot, because of the allocations of x and y random vectors at each iteration of the loop. These allocations can be solved by using static allocations, with tuples:
julia> function meandist(n)
d = 0.
for i in 1:n
x = (rand(),rand())
y = (rand(),rand())
d += sqrt((x[1]-y[1])^2+(x[2]-y[2])^2)
end
return d/n
end
meandist (generic function with 1 method)
julia> @btime meandist(1000)
13.012 μs (0 allocations: 0 bytes)
0.5122565299210201
or with the StaticArrays package:
julia> using StaticArrays
julia> function meandist(n)
d = 0.
for i in 1:n
x = SVector{2,Float64}(rand(),rand())
y = SVector{2,Float64}(rand(),rand())
d += sqrt((x[1]-y[1])^2+(x[2]-y[2])^2)
end
return d/n
end
meandist (generic function with 1 method)
julia> @btime meandist(1000)
13.027 μs (0 allocations: 0 bytes)
0.5096925569126717
These options are not as fast as the one that uses the Distances package because there the computation of the distances is loop-vectorized. We can do that “by hand” (edit: without manually using @avx the performance is the same - I guess loop vectorization occurs anyway, automatically):
julia> using LoopVectorization
julia> function meandist(n)
d = 0.
x = rand(2,n)
y = rand(2,n)
@avx for i in 1:n
d += sqrt((x[1,i]-y[1,i])^2+(x[2,i]-y[2,i])^2)
end
return d/n
end
meandist (generic function with 1 method)
julia> @btime meandist(1000)
5.006 μs (2 allocations: 31.50 KiB)
0.5336431684976665
and, of course, lot of that time is just the generation of the two matrices:
julia> function meandist(n,x,y)
d = 0.
@avx for i in 1:n
d += sqrt((x[1,i]-y[1,i])^2+(x[2,i]-y[2,i])^2)
end
return d/n
end
meandist (generic function with 2 methods)
julia> x = rand(2,1000); y = rand(2,1000);
julia> @btime meandist(1000,$x,$y)
1.530 μs (0 allocations: 0 bytes)
0.5345441925830112
Edit: Without generating the list of random points inside the function, the option with static arrays is as fast as the loop-vectorized (perhaps it gets loop-vectorized by the compiler automatically):
julia> function meandist(n,x,y)
d = 0.
for i in 1:n
d += sqrt((x[i][1]-y[i][1])^2+(x[i][2]-y[i][2])^2)
end
return d/n
end
meandist (generic function with 2 methods)
julia> x = rand(SVector{2,Float64},1000)
julia> y = rand(SVector{2,Float64},1000)
julia> @btime meandist(1000,$x,$y)
1.508 μs (0 allocations: 0 bytes)
0.5239587680180073
The equivalent version without static arrays is slower, but not bad:
julia> x = [ rand(2) for i in 1:1000 ];
julia> y = [ rand(2) for i in 1:1000 ];
julia> @btime meandist(1000,$x,$y)
3.292 μs (0 allocations: 0 bytes)
0.5311847109064113
Edit 2:
Given that allocating the vectors outside the function gives similar results for static arrays or matrices, the difference is in the time required to generate those. There are important differences there, indeed:
julia> @btime [ SVector{2,Float64}(rand(),rand()) for i in 1:1000 ];
12.237 μs (1 allocation: 15.75 KiB)
julia> @btime rand(SVector{2,Float64},1000);
3.195 μs (1 allocation: 15.75 KiB)
julia> @btime rand(2,1000);
1.605 μs (1 allocation: 15.75 KiB)