Implementation of closest-points D&C algorithm

First, we generate the points.

import numpy as np
import plotly.graph_objects as go
import plotly.io as pio
pio.renderers.default = 'notebook'

def generate_2d_points(N):
    return 200*np.random.rand(N,2)-100

def plot_2d_points(points):
    fig = go.Figure(data=go.Scatter(
        x=points[:, 0],
        y=points[:, 1],
        mode='markers',
        marker=dict(
            size=10,
            opacity=0.7
        )
    ))

    fig.update_layout(width=600, height=600)
    return fig

N = 50
points = generate_2d_points(N)
fig = plot_2d_points(points)
fig.show()

Brute Force Implementation

The brute force implementation has no surprises. It returns a tuple of the two points, and the distance between them.

def brute_force(points):
    bestpair = (points[0], points[1], np.linalg.norm(points[0]-points[1]))
    N, _ = points.shape

    for i in range(N-1):
        for j in range(i+1,N):
            dist = np.linalg.norm( points[i,:]-points[j,:] )
            if dist<bestpair[2]:
                bestpair = (points[i], points[j], dist.item())
    return bestpair

first, second, distance = brute_force(points)
print(first,second,distance)

[ 43.40312713 -22.57030696] [ 41.21756384 -26.76626592] 4.731042007145598

def highlight_closest(fig, point1, point2):
    fig.add_trace(go.Scatter(
        x=[point1[0], point2[0]],
        y=[point1[1], point2[1]],
        mode='markers',
        marker=dict(
            color='red',
            size=12
        )
    ))
    fig.update_layout(showlegend=False)
    return fig

fig = highlight_closest(fig, first, second )
fig.show()

Divide and Conquer Implementation

Again, the algorithm is:

• Presort the points by the x-axis.

• Recursively divide by the median x-value. When the number of points gets small, solve using the brute force approach.

• When combining, the solution is either the solution from the left half, the solution from the right half, or a point pair that straddles the boundary.

  • To find pairs that straddle the boundary, we first select only the points closer to the boundary than the closest pair from the left or right halves

  • We sort these points in that strip by the y-value

  • Each point only needs to be compared to a small, constant number of nearby points

def strip_closest(strip_points, min_dist_so_far):
    min_dist = min_dist_so_far
    closest_pair = (None, None)
    
    for i in range(strip_points.shape[0]):
        for j in range(i + 1, min(i + 8, strip_points.shape[0])):
            p1 = strip_points[i]
            p2 = strip_points[j]
            dist = np.linalg.norm(p1-p2)
            
            if dist < min_dist:
                min_dist = dist
                closest_pair = (p1, p2)
                
    return closest_pair[0], closest_pair[1], min_dist

def find_closest_pair_dc(points_sorted_x):
    N = points_sorted_x.shape[0]
    
    # Base Case: If N <= 3, use brute force
    if N <= 3:
        return brute_force(points_sorted_x)

    # Divide Step
    mid = N // 2
    mid_x = points_sorted_x[mid, 0]
    
    p1_l, p2_l, dist_l = find_closest_pair_dc(points_sorted_x[:mid])
    p1_r, p2_r, dist_r = find_closest_pair_dc(points_sorted_x[mid:])

    if dist_l < dist_r:
        min_dist, closest_p1, closest_p2 = dist_l, p1_l, p2_l
    else:
        min_dist, closest_p1, closest_p2 = dist_r, p1_r, p2_r
    
    strip = points_sorted_x[np.abs(points_sorted_x[:, 0] - mid_x) < min_dist]
    
    strip_sorted_y = strip[strip[:, 1].argsort()]

    p1_s, p2_s, dist_strip = strip_closest(strip_sorted_y, min_dist)
    
    if dist_strip < min_dist:
        return p1_s, p2_s, dist_strip
    else:
        return closest_p1, closest_p2, min_dist

def closest_pair_starter(points):
    if points.shape[0] < 2:
        return np.inf, None, None
    
    points_sorted_x = points[points[:, 0].argsort()]
    return find_closest_pair_dc(points_sorted_x)

first, second, dist = closest_pair_starter(points)
print(first, second, dist)

[ 41.21756384 -26.76626592] [ 43.40312713 -22.57030696] 4.731042007145598

fig = plot_2d_points(points)
fig = highlight_closest(fig, first, second)
fig.show()

Comparison

With only 5000 points, the O(n^2) and O(n log n) algorithms have a significant real difference in speed.

import time

N = 5000
points = generate_2d_points(N)

start = time.time()
first, second, distance = brute_force(points)
end = time.time()
print(f'Best points are {distance} apart, and brute force took {end-start} seconds.')

start = time.time()
first, second, distance = closest_pair_starter(points)
end = time.time()
print(f'Best points are {distance} apart, and D&C took {end-start} seconds.')

Best points are 0.0196075187992648 apart, and brute force took 16.488074779510498 seconds.
Best points are 0.0196075187992648 apart, and D&C took 0.11477851867675781 seconds.