# Closest distance between two Euclidian grids
**URL:** https://discourse.julialang.org/t/closest-distance-between-two-euclidian-grids/61218
**Category:** General Usage
**Tags:** question
**Created:** [May 15, 2021, 11:46pm UTC](https://discourse.julialang.org/t/closest-distance-between-two-euclidian-grids/61218 "2021-05-15T23:46:49Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@cadojo](https://discourse.julialang.org/u/cadojo)
#### Post date: [May 15, 2021, 11:46pm UTC](https://discourse.julialang.org/t/closest-distance-between-two-euclidian-grids/61218/1 "2021-05-15T23:46:49Z")
I have two large data sets: `A (size 1_500_000, 3)`, and `B (size 1_500, 3)`. I’d like to find the **closest point** between `A` and `B`. I could iterate through each possible combination, but that is `O(n^2)` and is not possible to complete in a reasonable time for datasets this size.
I asked a similar question for one dimension [here](https://discourse.julialang.org/t/approximate-set-intersection-for-floating-point-numbers/61047). I understand sorting the data would be helpful, but I’m not sure how to apply this to a 2-dimensional case.
Is there a package that can compute the **closest point** between two Euclidian grids?
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [May 16, 2021, 12:01am UTC](https://discourse.julialang.org/t/closest-distance-between-two-euclidian-grids/61218/2 "2021-05-16T00:01:05Z")
Did you check [NearestNeighbors.jl](https://www.google.com/url?sa=t&source=web&rct=j&url=https://github.com/KristofferC/NearestNeighbors.jl&ved=2ahUKEwimmOLj8czwAhWFh_0HHX68DDkQFjABegQIAxAC&usg=AOvVaw0y9f4HgSdHZVLI-l1kzYLy)?
Including herein a MWE, that hopefully is correct:
```julia
using NearestNeighbors, Random, Printf
# create input data:
Random.seed!(1); A = rand(3, 1_500_000)
Random.seed!(1234); B = rand(3, 1_500)
kdtree = KDTree(A)
idxs, dists = nn(kdtree, B) # convenience function for k=1
ixb = argmin(dists)
ixa = idxs[ixb]
@printf("Minimum distance = %.5f, at:", dists[ixb])
@printf("A[:,%i] = (%.4f, %.4f, %.4f)\n", ixa, A[:,ixa]...) # splatted values supported
@printf("B[:,%i] = (%.4f, %.4f, %.4f)\n", ixb, B[:,ixb]...) #
# Printed results:
Minimum distance = 0.00045, at:
A[:,1080399] = (0.3572, 0.1117, 0.4786)
B[:,1136] = (0.3571, 0.1119, 0.4790)
# QC solution:
using LinearAlgebra
julia> dists[ixb] ≈ norm(A[:,ixa] - B[:,ixb])
true
```
---
### Author:  [@lmiq](https://discourse.julialang.org/u/lmiq)
#### Post date: [May 16, 2021, 1:32am UTC](https://discourse.julialang.org/t/closest-distance-between-two-euclidian-grids/61218/3 "2021-05-16T01:32:55Z")
In the field of molecular simulations we always rely on [Cell lists](https://en.m.wikipedia.org/wiki/Cell_lists) to compute short range distances. That could be used to compute the nearest neighbor in roughly O(n) time.\*
I am curious about the algorithms implemented in the package above. Anyone can point a reference of the algorithms implemented there?
\*It has been done: “Near-neighbor calculations using a modified cell-linked list method - ScienceDirect” [Near-neighbor calculations using a modified cell-linked list method - ScienceDirect](https://www.sciencedirect.com/science/article/abs/pii/S0010465598002033)
Edit: I understand this is the type of thing NearestNeighbors.jl is doing: [k-d tree - Wikipedia](https://en.m.wikipedia.org/wiki/K-d_tree)
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [May 16, 2021, 1:47pm UTC](https://discourse.julialang.org/t/closest-distance-between-two-euclidian-grids/61218/4 "2021-05-16T13:47:24Z")
Do not know if the [nearest neighbor algorithm](https://en.wikipedia.org/wiki/K-nearest_neighbors_algorithm) is among the most efficient to answer OP’s question.
In the MWE example above, it seems that one may need two steps: (i) find in A the 1500 nearest neighbors of B, and (ii) find the nearest among the nearest.
There may be a better solution while using NearestNeighbors.jl.
---
### Author:  [@cadojo](https://discourse.julialang.org/u/cadojo)
#### Post date: [May 16, 2021, 4:16pm UTC](https://discourse.julialang.org/t/closest-distance-between-two-euclidian-grids/61218/5 "2021-05-16T16:16:33Z")
Love the Julia community, this is perfect. Nearest neighbor is found in just a few seconds, where my previous O(n^2) implementation never returned. Thanks for pointing this out!
---
### Author:  [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)
#### Post date: [May 16, 2021, 7:24pm UTC](https://discourse.julialang.org/t/closest-distance-between-two-euclidian-grids/61218/6 "2021-05-16T19:24:21Z")
Isn’t this two one-dimensional searches really?
---
### Author:  [@cadojo](https://discourse.julialang.org/u/cadojo)
#### Post date: [May 16, 2021, 7:34pm UTC](https://discourse.julialang.org/t/closest-distance-between-two-euclidian-grids/61218/7 "2021-05-16T19:34:15Z")
Might be missing something, but I don’t see how.
I have a collection of departure points that I want to travel _from_: `A`. And I have a collection of arrival points that I want to travel _to_: `B`.
I want to find the pair of points (one from `A`, one from `B`) that are the **closest** together. The naive solution for this is…
```julia
minimum = Inf
optimal = ([NaN, NaN, NaN], [NaN, NaN, NaN])
for departure in A
for arrival in B
difference = norm(departure .- arrival)
if difference < minimum
minimum = difference
optimal = (departure, arrival)
end
end
end
```
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [May 16, 2021, 7:39pm UTC](https://discourse.julialang.org/t/closest-distance-between-two-euclidian-grids/61218/8 "2021-05-16T19:39:32Z")
Please check [this article](https://en.wikipedia.org/wiki/Closest_pair_of_points_problem) out, that points to a Julia solution [here](http://rosettacode.org/wiki/Closest-pair_problem#Julia).
_ **NB:** unfortunately, the code posted there is labelled “Brute-force algorithm”_
---
### Author:  [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)
#### Post date: [May 16, 2021, 7:40pm UTC](https://discourse.julialang.org/t/closest-distance-between-two-euclidian-grids/61218/9 "2021-05-16T19:40:28Z")
Sort in X, Y, Z using the Manhattan norm. Then sequential search.
---
### Author:  [@lmiq](https://discourse.julialang.org/u/lmiq)
#### Post date: [May 16, 2021, 7:54pm UTC](https://discourse.julialang.org/t/closest-distance-between-two-euclidian-grids/61218/10 "2021-05-16T19:54:12Z")
If I understand what the Manhattan norm is, computing it alone would already be a O(n\*m) computation.
(The algorithms must avoid looping over all pairs of points, and for that some sort of spacial classification is needed)
@cadojo one curiosity: since your points are three dimensional, I thought you were dealing with atoms. What kind of travel that is though? Between galaxies?
---
### Author:  [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)
#### Post date: [May 16, 2021, 8:33pm UTC](https://discourse.julialang.org/t/closest-distance-between-two-euclidian-grids/61218/11 "2021-05-16T20:33:00Z")
I misunderstood the OP.
---
### Author:  [@lmiq](https://discourse.julialang.org/u/lmiq)
#### Post date: [May 16, 2021, 9:52pm UTC](https://discourse.julialang.org/t/closest-distance-between-two-euclidian-grids/61218/12 "2021-05-16T21:52:34Z")
Just curiosity, I implemented a (non-general, not very optimized) form of the cell lists above. It resulted to be be ~8 times slower than the method implemented in `NearestNeighbors.jl`. I will have to check if the method there applied to my actual use case 🙂 . I hope it does, having it 8 times faster won’t be bad at all.
```julia
NN(A, B) = (0.00044736197216615207, 1080399, 1136)
0.918342 seconds (7.54 k allocations: 96.740 MiB)
MD(A, B) = (0.000447361972166152, 1080399, 1136)
7.259169 seconds (4 allocations: 11.452 MiB, 0.11% gc time)
(0.000447361972166152, 1080399, 1136)
```
```julia
import Pkg
Pkg.activate("nn")
using BenchmarkTools
using NearestNeighbors, Random, Printf
N1 = 15 * 10^5
N2 = 15 * 10^2
# create input data:
Random.seed!(1); A = rand(3, N1)
Random.seed!(1234); B = rand(3, N2)
# NearestNeighbors algorithm:
function NN(A,B)
kdtree = KDTree(A)
idxs, dists = nn(kdtree, B) # convenience function for k=1
ixb = argmin(dists)
ixa = idxs[ixb]
return dists[ixb], ixa, ixb
end
@show NN(A,B)
@time NN(A,B)
function icell3D(x, n)
i = trunc(Int,(x[1]/n))+1
j = trunc(Int,(x[2]/n))+1
k = trunc(Int,(x[3]/n))+1
return i, j, k
end
function MD(A,B)
d = 0.1 # assume that we know that the shortest distance is smaller than this
n = floor(Int,1/d)
# cell lists
fatomA = zeros(Int,n,n,n) # first atom of A in the cell
natomA = zeros(Int,size(A,2)) # next atom of A in the cell
# fill cells
for i in axes(A,2)
ic, jc, kc = icell3D(@view(A[:,i]),n)
natomA[i] = fatomA[ic,jc,kc]
fatomA[ic,jc,kc] = i
end
iamin = 0
ibmin = 0
dmin = +Inf
# loop over atoms of B
@inbounds for ib in axes(B,2)
xb = @view(B[:,ib])
i, j, k = icell3D(xb,n)
# Loop over vicinal cells to compute distances to solvent atoms, and
# add data to dc structure (includes current cell)
for ic in i-1:i+1
(ic == 0 || ic > n) && continue
for jc in j-1:j+1
(jc == 0 || jc > n) && continue
for kc in k-1:k+1
(kc == 0 || kc > n) && continue
ia = fatomA[ic,jc,kc]
while ia > 0
dij = (A[1,ia]-B[1,ib])^2 +
(A[2,ia]-B[2,ib])^2 +
(A[3,ia]-B[3,ib])^2
if dij < dmin
dmin = dij
iamin = ia
ibmin = ib
end
ia = natomA[ia]
end
end
end
end
end
return sqrt(dmin), iamin, ibmin
end
@show MD(A,B)
@time MD(A,B)
```
---
### Author:  [@cadojo](https://discourse.julialang.org/u/cadojo)
#### Post date: [May 17, 2021, 7:58pm UTC](https://discourse.julialang.org/t/closest-distance-between-two-euclidian-grids/61218/13 "2021-05-17T19:58:07Z")
> What kind of travel that is though? Between galaxies?
Travel in our solar system! I’m using invariant manifolds within the Circular Restricted Three-body Problem to find low-cost transfers from Earth to Jupiter. [Repo here!](https://github.com/cadojo/Halo-Orbit-Explorations)
---
### Author:  [@lmiq](https://discourse.julialang.org/u/lmiq)
#### Post date: [May 17, 2021, 9:04pm UTC](https://discourse.julialang.org/t/closest-distance-between-two-euclidian-grids/61218/14 "2021-05-17T21:04:17Z")
Cool! And why there are so many objects?
---
### Author:  [@cadojo](https://discourse.julialang.org/u/cadojo)
#### Post date: [May 17, 2021, 9:31pm UTC](https://discourse.julialang.org/t/closest-distance-between-two-euclidian-grids/61218/15 "2021-05-17T21:31:06Z")
I’m trying to find the intersection of a manifold that **departs** the Sun-Earth system, and a manifold that **converges to** a Halo orbit in the Sun-Jupiter system. I’d ideally like to find a machine-level precise intersection between the manifold positions. Once I find two positions (one in each manifold) that are nearly identical, I can treat the difference in **velocity** between the two overlapping points in the manifolds as an impulsive maneuver.
I don’t have any analytical bounds for the manifolds, I just have the raw manifold positions and velocities that I can compute. So, to find a precise intersection, I’m collection **a lot** of data points, with the hope that there’s a better chance there is an intersection.
Unfortunately, it’s seeming like the closest approach between the two manifolds I can find (using a numerical optimizer with Sun-Jupiter Halo orbit amplitude as the single parameter) is around 3000 km. This is **not** machine level precision! I could use interpolation between data points, but I’ve received guidance that interpolation isn’t ideal here.
So my current idea is to build a better optimization. I’m thinking of collection **a lot** of orbital states within the Sun-Earth Halo, and the Sun-Jupiter Halo. I’ll then use a box optimization with the following parameters:
1. Phase angle of the Sun-Earth Halo departure (onto the departure manifold)
2. Phase angle of the Sun-Jupiter Halo departure (onto the arrival manifold)
3. Sun-Earth departure perturbation amplitude
4. Sun-Jupiter arrival perturbation amplitude
So rather than compute several hundred trajectories off of each Halo orbit, and use nearest neighbors to find the closest approach between the two manifolds, I’ll have the optimizer…
1. Pick parameters
2. Propagate the departure and arrival trajectories for some pre-set duration
3. Find the closest approach between the two trajectories
4. Return the distance between the two closest orbital states as the _cost_ for the optimizer to minimize
If anyone has a better idea for how to do this, it would be much appreciated!
As a brief self-plug, all of this functionality is available in `GeneralAstrodynamics` on [GitHub](https://github.com/cadojo/GeneralAstrodynamics.jl) and Julia’s [general registry](https://juliahub.com/ui/Packages/GeneralAstrodynamics/o9DTc). I’m planning on polishing this up and doing a package announcement for a 1.0 release this summer!
---
### Author:  [@cortner](https://discourse.julialang.org/u/cortner)
#### Post date: [May 17, 2021, 10:08pm UTC](https://discourse.julialang.org/t/closest-distance-between-two-euclidian-grids/61218/16 "2021-05-17T22:08:34Z")
FWIW cell lists are implemented in
> **[GitHub - JuliaMolSim/NeighbourLists.jl: neighbour list for particle...](https://github.com/JuliaMolSim/NeighbourLists.jl)**
>
> neighbour list for particle simulations based on matscipy - GitHub - JuliaMolSim/NeighbourLists.jl: neighbour list for particle simulations based on matscipy
---
### Author:  [@lmiq](https://discourse.julialang.org/u/lmiq)
#### Post date: [May 17, 2021, 11:23pm UTC](https://discourse.julialang.org/t/closest-distance-between-two-euclidian-grids/61218/17 "2021-05-17T23:23:52Z")
I implemented cell lists in the [ComplexMixtures](https://github.com/m3g/ComplexMixtures.jl) package, with support to periodic boundaries (with straight angles only). Do you have some easy benchmark so I can compare?
I should have made it more shearable, but when I started that I didn’t understand the modularity that Julia offers. Now I have to do some rewriting to make that useful to others if that is the case.
---
### Author:  [@cortner](https://discourse.julialang.org/u/cortner)
#### Post date: [May 17, 2021, 11:32pm UTC](https://discourse.julialang.org/t/closest-distance-between-two-euclidian-grids/61218/18 "2021-05-17T23:32:52Z")
I used to have some benchmarks to compare against matscipy. I might be able to revive them. Would be great to combine codes (or get rid of one) since mine is still missing some functionality
---
### Author:  [@cortner](https://discourse.julialang.org/u/cortner)
#### Post date: [May 17, 2021, 11:33pm UTC](https://discourse.julialang.org/t/closest-distance-between-two-euclidian-grids/61218/19 "2021-05-17T23:33:32Z")
PS - it’s fairly well debugged and should Desk with arbitrary cell shapes. Not sure whether it is robust when the cell collapses though
---
### Author:  [@lmiq](https://discourse.julialang.org/u/lmiq)
#### Post date: [June 5, 2021, 3:20am UTC](https://discourse.julialang.org/t/closest-distance-between-two-euclidian-grids/61218/20 "2021-06-05T03:20:28Z")
Just if you are interested. I have put up a package to compute function mappings dependent on short-ranged pairwise particle distances, and it can be used to compute the minimum distance between two sets of points. Depending on what you know from your point distribution, it can be faster than the nearest neighbour search above. It is here: [https://github.com/m3g/CellListMap.jl](https://github.com/m3g/CellListMap.jl)
I added how to compute the same thing above in the example, and got (with 4 threads):
```julia
NN(A, B) = (0.00044736197216615207, 1080399, 1136)
987.730 ms (7541 allocations: 96.74 MiB)
getmind(x, y) = (0.000447361972166152, 1136, 1080399)
71.849 ms (141 allocations: 22.90 MiB)
```
single threaded I get 266.265 ms for getmind.
The “box side” and the “cutoff” must be tunned with the knowledge of your system. In particular, if you do not have periodic boundary conditions, as was the case here, I had to set the box side to be larger than the maximum dimensions plus the cutoff. Of course the time is very much dependent on the cutoff used, and there is some optimal value.
> **Code**
>
> ```julia
> import Pkg
> Pkg.activate("nn")
>
> #Pkg.add(url="https://github.com/m3g/CellListMap.jl")
> using CellListMap, StaticArrays
> using BenchmarkTools, Test
> using NearestNeighbors, Random, Printf
>
> N1 = 15 * 10^5
> N2 = 15 * 10^2
>
> # create input data:
> Random.seed!(1); A = rand(3, N1)
> Random.seed!(1234); B = rand(3, N2)
>
> function NN(A,B)
> kdtree = KDTree(A)
> idxs, dists = nn(kdtree, B) # convenience function for k=1
> ixb = argmin(dists)
> ixa = idxs[ixb]
> return dists[ixb], ixa, ixb
> end
>
> x = [SVector{3,Float64}(B[:,i]) for i in 1:N2 ]
> y = [SVector{3,Float64}(A[:,i]) for i in 1:N1 ]
>
> function getmind(x,y)
>
> # Number of particles, sides and cutoff
> sides = [1.2,1.2,1.2]
> cutoff = 0.03
> box = Box(sides,cutoff)
>
> # Initialize auxiliary linked lists (largest set!)
> lc = LinkedLists(length(y))
>
> # Initializing linked cells with these positions (largest set!)
> initcells!(y,box,lc)
>
> # Function that keeps the minimum distance
> f(i,j,d2,mind) = d2 < mind[1] ? (d2,i,j) : mind
>
> # We have to define our own reduce function here for the tuple
> function reduce_mind(output_threaded)
> mind = output_threaded[1]
> for i in 2:Threads.nthreads()
> if output_threaded[i][1] < mind[1]
> mind = output_threaded[i]
> end
> end
> return mind
> end
>
> # Result
> mind = ( +Inf, 0, 0 )
>
> # Run pairwise computation
> mind = map_pairwise!(
> (x,y,i,j,d2,out) -> f(i,j,d2,out),
> mind,x,y,box,lc;reduce=reduce_mind
> )
> return (sqrt(mind[1]),mind[2],mind[3])
>
> end
>
> @show NN(A,B)
> @show getmind(x,y)
>
> @btime NN(A,B)
> @btime getmind(x,y)
>
> nothing
>
> ```