Closest-pair problem

From Rosetta Code
(Redirected from Closest pair problem)
Jump to: navigation, search
Task
Closest-pair problem
You are encouraged to solve this task according to the task description, using any language you may know.
This page uses content from Wikipedia. The original article was at Closest pair of points problem. The list of authors can be seen in the page history. As with Rosetta Code, the text of Wikipedia is available under the GNU FDL. (See links for details on variance)

The aim of this task is to provide a function to find the closest two points among a set of given points in two dimensions, i.e. to solve the Closest pair of points problem in the planar case.

The straightforward solution is a O(n2) algorithm (which we can call brute-force algorithm); the pseudocode (using indexes) could be simply:

bruteForceClosestPair of P(1), P(2), ... P(N)
if N < 2 then
  returnelse
  minDistance ← |P(1) - P(2)|
  minPoints ← { P(1), P(2) }
  foreach i ∈ [1, N-1]
    foreach j ∈ [i+1, N]
      if |P(i) - P(j)| < minDistance then
        minDistance ← |P(i) - P(j)|
        minPoints ← { P(i), P(j) } 
      endif
    endfor
  endfor
  return minDistance, minPoints
 endif

A better algorithm is based on the recursive divide&conquer approach, as explained also at Wikipedia, which is O(n log n); a pseudocode could be:

closestPair of (xP, yP)
               where xP is P(1) .. P(N) sorted by x coordinate, and
                     yP is P(1) .. P(N) sorted by y coordinate (ascending order)
if N ≤ 3 then
  return closest points of xP using brute-force algorithm
else
  xL ← points of xP from 1 to ⌈N/2⌉
  xR ← points of xP from ⌈N/2⌉+1 to N
  xm ← xP(⌈N/2⌉)x
  yL ← { p ∈ yP : px ≤ xm }
  yR ← { p ∈ yP : px > xm }
  (dL, pairL) ← closestPair of (xL, yL)
  (dR, pairR) ← closestPair of (xR, yR)
  (dmin, pairMin) ← (dR, pairR)
  if dL < dR then
    (dmin, pairMin) ← (dL, pairL)
  endif
  yS ← { p ∈ yP : |xm - px| < dmin }
  nS ← number of points in yS
  (closest, closestPair) ← (dmin, pairMin)
  for i from 1 to nS - 1
    k ← i + 1
    while k ≤ nS and yS(k)y - yS(i)y < dmin
      if |yS(k) - yS(i)| < closest then
        (closest, closestPair) ← (|yS(k) - yS(i)|, {yS(k), yS(i)})
      endif
      k ← k + 1
    endwhile
  endfor
  return closest, closestPair
endif


References and further readings

Contents

[edit] Ada

Dimension independant, but has to be defined at procedure call time (could be a parameter). Output is simple, can be formatted using Float_IO.

closest.adb: (uses brute force algorithm)

with Ada.Numerics.Generic_Elementary_Functions;
with Ada.Text_IO;
 
procedure Closest is
package Math is new Ada.Numerics.Generic_Elementary_Functions (Float);
 
Dimension : constant := 2;
type Vector is array (1 .. Dimension) of Float;
type Matrix is array (Positive range <>) of Vector;
 
-- calculate the distance of two points
function Distance (Left, Right : Vector) return Float is
Result : Float := 0.0;
Offset : Natural := 0;
begin
loop
Result := Result + (Left(Left'First + Offset) - Right(Right'First + Offset))**2;
Offset := Offset + 1;
exit when Offset >= Left'Length;
end loop;
return Math.Sqrt (Result);
end Distance;
 
-- determine the two closest points inside a cloud of vectors
function Get_Closest_Points (Cloud : Matrix) return Matrix is
Result : Matrix (1..2);
Min_Distance : Float;
begin
if Cloud'Length(1) < 2 then
raise Constraint_Error;
end if;
Result := (Cloud (Cloud'First), Cloud (Cloud'First + 1));
Min_Distance := Distance (Cloud (Cloud'First), Cloud (Cloud'First + 1));
for I in Cloud'First (1) .. Cloud'Last(1) - 1 loop
for J in I + 1 .. Cloud'Last(1) loop
if Distance (Cloud (I), Cloud (J)) < Min_Distance then
Min_Distance := Distance (Cloud (I), Cloud (J));
Result := (Cloud (I), Cloud (J));
end if;
end loop;
end loop;
return Result;
end Get_Closest_Points;
 
Test_Cloud : constant Matrix (1 .. 10) := ( (5.0, 9.0), (9.0, 3.0),
(2.0, 0.0), (8.0, 4.0),
(7.0, 4.0), (9.0, 10.0),
(1.0, 9.0), (8.0, 2.0),
(0.0, 10.0), (9.0, 6.0));
Closest_Points : Matrix := Get_Closest_Points (Test_Cloud);
 
Second_Test : constant Matrix (1 .. 10) := ( (0.654682, 0.925557), (0.409382, 0.619391),
(0.891663, 0.888594), (0.716629, 0.9962),
(0.477721, 0.946355), (0.925092, 0.81822),
(0.624291, 0.142924), (0.211332, 0.221507),
(0.293786, 0.691701), (0.839186, 0.72826));
Second_Points : Matrix := Get_Closest_Points (Second_Test);
begin
Ada.Text_IO.Put_Line ("Closest Points:");
Ada.Text_IO.Put_Line ("P1: " & Float'Image (Closest_Points (1) (1)) & " " & Float'Image (Closest_Points (1) (2)));
Ada.Text_IO.Put_Line ("P2: " & Float'Image (Closest_Points (2) (1)) & " " & Float'Image (Closest_Points (2) (2)));
Ada.Text_IO.Put_Line ("Distance: " & Float'Image (Distance (Closest_Points (1), Closest_Points (2))));
Ada.Text_IO.Put_Line ("Closest Points 2:");
Ada.Text_IO.Put_Line ("P1: " & Float'Image (Second_Points (1) (1)) & " " & Float'Image (Second_Points (1) (2)));
Ada.Text_IO.Put_Line ("P2: " & Float'Image (Second_Points (2) (1)) & " " & Float'Image (Second_Points (2) (2)));
Ada.Text_IO.Put_Line ("Distance: " & Float'Image (Distance (Second_Points (1), Second_Points (2))));
end Closest;
Output:
Closest Points:
P1:  8.00000E+00  4.00000E+00
P2:  7.00000E+00  4.00000E+00
Distance:  1.00000E+00
Closest Points 2:
P1:  8.91663E-01  8.88594E-01
P2:  9.25092E-01  8.18220E-01
Distance:  7.79101E-02

[edit] AWK

 
# syntax: GAWK -f CLOSEST-PAIR_PROBLEM.AWK
BEGIN {
x[++n] = 0.654682 ; y[n] = 0.925557
x[++n] = 0.409382 ; y[n] = 0.619391
x[++n] = 0.891663 ; y[n] = 0.888594
x[++n] = 0.716629 ; y[n] = 0.996200
x[++n] = 0.477721 ; y[n] = 0.946355
x[++n] = 0.925092 ; y[n] = 0.818220
x[++n] = 0.624291 ; y[n] = 0.142924
x[++n] = 0.211332 ; y[n] = 0.221507
x[++n] = 0.293786 ; y[n] = 0.691701
x[++n] = 0.839186 ; y[n] = 0.728260
min = 1E20
for (i=1; i<=n-1; i++) {
for (j=i+1; j<=n; j++) {
dsq = (x[i]-x[j])^2 + (y[i]-y[j])^2
if (dsq < min) {
min = dsq
mini = i
minj = j
}
}
}
printf("distance between (%.6f,%.6f) and (%.6f,%.6f) is %g\n",x[mini],y[mini],x[minj],y[minj],sqrt(min))
exit(0)
}
 
Output:
distance between (0.891663,0.888594) and (0.925092,0.818220) is 0.0779102

[edit] BBC BASIC

To find the closest pair it is sufficient to compare the squared-distances, it is not necessary to perform the square root for each pair!

      DIM x(9), y(9)
 
FOR I% = 0 TO 9
READ x(I%), y(I%)
NEXT
 
min = 1E30
FOR I% = 0 TO 8
FOR J% = I%+1 TO 9
dsq = (x(I%) - x(J%))^2 + (y(I%) - y(J%))^2
IF dsq < min min = dsq : mini% = I% : minj% = J%
NEXT
NEXT I%
PRINT "Closest pair is ";mini% " and ";minj% " at distance "; SQR(min)
END
 
DATA 0.654682, 0.925557
DATA 0.409382, 0.619391
DATA 0.891663, 0.888594
DATA 0.716629, 0.996200
DATA 0.477721, 0.946355
DATA 0.925092, 0.818220
DATA 0.624291, 0.142924
DATA 0.211332, 0.221507
DATA 0.293786, 0.691701
DATA 0.839186, 0.728260
 
Output:
Closest pair is 2 and 5 at distance 0.0779101913

[edit] C

See Closest-pair problem/C

[edit] C++

/*
Author: Kevin Bacon
Date: 04/03/2014
Task: Closest-pair problem
*/

 
#include <iostream>
#include <vector>
#include <utility>
#include <cmath>
#include <random>
#include <chrono>
#include <algorithm>
#include <iterator>
 
typedef std::pair<double, double> point_t;
typedef std::pair<point_t, point_t> points_t;
 
double distance_between(const point_t& a, const point_t& b) {
return std::sqrt(std::pow(b.first - a.first, 2)
+ std::pow(b.second - a.second, 2));
}
 
std::pair<double, points_t> find_closest_brute(const std::vector<point_t>& points) {
if (points.size() < 2) {
return { -1, { { 0, 0 }, { 0, 0 } } };
}
auto minDistance = std::abs(distance_between(points.at(0), points.at(1)));
points_t minPoints = { points.at(0), points.at(1) };
for (auto i = std::begin(points); i != (std::end(points) - 1); ++i) {
for (auto j = i + 1; j < std::end(points); ++j) {
auto newDistance = std::abs(distance_between(*i, *j));
if (newDistance < minDistance) {
minDistance = newDistance;
minPoints.first = *i;
minPoints.second = *j;
}
}
}
return { minDistance, minPoints };
}
 
std::pair<double, points_t> find_closest_optimized(const std::vector<point_t>& xP,
const std::vector<point_t>& yP) {
if (xP.size() <= 3) {
return find_closest_brute(xP);
}
auto N = xP.size();
auto xL = std::vector<point_t>();
auto xR = std::vector<point_t>();
std::copy(std::begin(xP), std::begin(xP) + (N / 2), std::back_inserter(xL));
std::copy(std::begin(xP) + (N / 2), std::end(xP), std::back_inserter(xR));
auto xM = xP.at(N / 2).first;
auto yL = std::vector<point_t>();
auto yR = std::vector<point_t>();
std::copy_if(std::begin(yP), std::end(yP), std::back_inserter(yL), [&xM](const point_t& p) {
return p.first <= xM;
});
std::copy_if(std::begin(yP), std::end(yP), std::back_inserter(yR), [&xM](const point_t& p) {
return p.first > xM;
});
auto p1 = find_closest_optimized(xL, yL);
auto p2 = find_closest_optimized(xR, yR);
auto minPair = (p1.first <= p2.first) ? p1 : p2;
auto yS = std::vector<point_t>();
std::copy_if(std::begin(yP), std::end(yP), std::back_inserter(yS), [&minPair, &xM](const point_t& p) {
return std::abs(xM - p.first) < minPair.first;
});
auto result = minPair;
for (auto i = std::begin(yS); i != (std::end(yS) - 1); ++i) {
for (auto k = i + 1; k != std::end(yS) &&
((k->second - i->second) < minPair.first); ++k) {
auto newDistance = std::abs(distance_between(*k, *i));
if (newDistance < result.first) {
result = { newDistance, { *k, *i } };
}
}
}
return result;
}
 
void print_point(const point_t& point) {
std::cout << "(" << point.first
<< ", " << point.second
<< ")";
}
 
int main(int argc, char * argv[]) {
std::default_random_engine re(std::chrono::system_clock::to_time_t(
std::chrono::system_clock::now()));
std::uniform_real_distribution<double> urd(-500.0, 500.0);
std::vector<point_t> points(100);
std::generate(std::begin(points), std::end(points), [&urd, &re]() {
return point_t { 1000 + urd(re), 1000 + urd(re) };
});
auto answer = find_closest_brute(points);
std::sort(std::begin(points), std::end(points), [](const point_t& a, const point_t& b) {
return a.first < b.first;
});
auto xP = points;
std::sort(std::begin(points), std::end(points), [](const point_t& a, const point_t& b) {
return a.second < b.second;
});
auto yP = points;
std::cout << "Min distance (brute): " << answer.first << " ";
print_point(answer.second.first);
std::cout << ", ";
print_point(answer.second.second);
answer = find_closest_optimized(xP, yP);
std::cout << "\nMin distance (optimized): " << answer.first << " ";
print_point(answer.second.first);
std::cout << ", ";
print_point(answer.second.second);
return 0;
}
Output:
Min distance (brute): 6.95886 (932.735, 1002.7), (939.216, 1000.17)
Min distance (optimized): 6.95886 (932.735, 1002.7), (939.216, 1000.17)

[edit] Clojure

 
(defn distance [[x1 y1] [x2 y2]]
(let [dx (- x2 x1), dy (- y2 y1)]
(Math/sqrt (+ (* dx dx) (* dy dy)))))
 
(defn brute-force [points]
(let [n (count points)]
(when (< 1 n)
(apply min-key first
(for [i (range 0 (dec n)), :let [p1 (nth points i)],
j (range (inc i) n), :let [p2 (nth points j)]]
[(distance p1 p2) p1 p2])))))
 
(defn combine [yS [dmin pmin1 pmin2]]
(apply min-key first
(conj (for [[p1 p2] (partition 2 1 yS)
 :let [[_ py1] p1 [_ py2] p2]
 :while (< (- py1 py2) dmin)]
[(distance p1 p2) p1 p2])
[dmin pmin1 pmin2])))
 
(defn closest-pair
([points]
(closest-pair
(sort-by first points)
(sort-by second points)))
([xP yP]
(if (< (count xP) 4)
(brute-force xP)
(let [[xL xR] (partition-all (Math/ceil (/ (count xP) 2)) xP)
[xm _] (last xL)
{yL true yR false} (group-by (fn [[px _]] (<= px xm)) yP)
dL&pairL (closest-pair xL yL)
dR&pairR (closest-pair xR yR)
[dmin pmin1 pmin2] (min-key first dL&pairL dR&pairR)
{yS true} (group-by (fn [[px _]] (< (Math/abs (- xm px)) dmin)) yP)]
(combine yS [dmin pmin1 pmin2])))))
 

[edit] Common Lisp

Points are conses whose cars are x coördinates and whose cdrs are y coördinates. This version includes the optimizations given in the McGill description of the algorithm.

(defun point-distance (p1 p2)
(destructuring-bind (x1 . y1) p1
(destructuring-bind (x2 . y2) p2
(let ((dx (- x2 x1)) (dy (- y2 y1)))
(sqrt (+ (* dx dx) (* dy dy)))))))
 
(defun closest-pair-bf (points)
(let ((pair (list (first points) (second points)))
(dist (point-distance (first points) (second points))))
(dolist (p1 points (values pair dist))
(dolist (p2 points)
(unless (eq p1 p2)
(let ((pdist (point-distance p1 p2)))
(when (< pdist dist)
(setf (first pair) p1
(second pair) p2
dist pdist))))))))
 
(defun closest-pair (points)
(labels
((cp (xp &aux (length (length xp)))
(if (<= length 3)
(multiple-value-bind (pair distance) (closest-pair-bf xp)
(values pair distance (sort xp '< :key 'cdr)))
(let* ((xr (nthcdr (1- (floor length 2)) xp))
(xm (/ (+ (caar xr) (caadr xr)) 2)))
(psetf xr (rest xr)
(rest xr) '())
(multiple-value-bind (lpair ldist yl) (cp xp)
(multiple-value-bind (rpair rdist yr) (cp xr)
(multiple-value-bind (dist pair)
(if (< ldist rdist)
(values ldist lpair)
(values rdist rpair))
(let* ((all-ys (merge 'vector yl yr '< :key 'cdr))
(ys (remove-if #'(lambda (p)
(> (abs (- (car p) xm)) dist))
all-ys))
(ns (length ys)))
(dotimes (i ns)
(do ((k (1+ i) (1+ k)))
((or (= k ns)
(> (- (cdr (aref ys k))
(cdr (aref ys i)))
dist)))
(let ((pd (point-distance (aref ys i)
(aref ys k))))
(when (< pd dist)
(setf dist pd
(first pair) (aref ys i)
(second pair) (aref ys k))))))
(values pair dist all-ys)))))))))
(multiple-value-bind (pair distance)
(cp (sort (copy-list points) '< :key 'car))
(values pair distance))))

[edit] C#

We provide a small helper class for distance comparisons:

class Segment
{
public Segment(PointF p1, PointF p2)
{
P1 = p1;
P2 = p2;
}
 
public readonly PointF P1;
public readonly PointF P2;
 
public float Length()
{
return (float)Math.Sqrt(LengthSquared());
}
 
public float LengthSquared()
{
return (P1.X - P2.X) * (P1.X - P2.X)
+ (P1.Y - P2.Y) * (P1.Y - P2.Y);
}
}

Brute force:

Segment Closest_BruteForce(List<PointF> points)
{
int n = points.Count;
var result = Enumerable.Range( 0, n-1)
.SelectMany( i => Enumerable.Range( i+1, n-(i+1) )
.Select( j => new Segment( points[i], points[j] )))
.OrderBy( seg => seg.LengthSquared())
.First();
 
return result;
}


And divide-and-conquer.

 
public static Segment MyClosestDivide(List<PointF> points)
{
return MyClosestRec(points.OrderBy(p => p.X).ToList());
}
 
private static Segment MyClosestRec(List<PointF> pointsByX)
{
int count = pointsByX.Count;
if (count <= 4)
return Closest_BruteForce(pointsByX);
 
// left and right lists sorted by X, as order retained from full list
var leftByX = pointsByX.Take(count/2).ToList();
var leftResult = MyClosestRec(leftByX);
 
var rightByX = pointsByX.Skip(count/2).ToList();
var rightResult = MyClosestRec(rightByX);
 
var result = rightResult.Length() < leftResult.Length() ? rightResult : leftResult;
 
// There may be a shorter distance that crosses the divider
// Thus, extract all the points within result.Length either side
var midX = leftByX.Last().X;
var bandWidth = result.Length();
var inBandByX = pointsByX.Where(p => Math.Abs(midX - p.X) <= bandWidth);
 
// Sort by Y, so we can efficiently check for closer pairs
var inBandByY = inBandByX.OrderBy(p => p.Y).ToArray();
 
int iLast = inBandByY.Length - 1;
for (int i = 0; i < iLast; i++ )
{
var pLower = inBandByY[i];
 
for (int j = i + 1; j <= iLast; j++)
{
var pUpper = inBandByY[j];
 
// Comparing each point to successivly increasing Y values
// Thus, can terminate as soon as deltaY is greater than best result
if ((pUpper.Y - pLower.Y) >= result.Length())
break;
 
if (Segment.Length(pLower, pUpper) < result.Length())
result = new Segment(pLower, pUpper);
}
}
 
return result;
}
 

However, the difference in speed is still remarkable.

var randomizer = new Random(10);
var points = Enumerable.Range( 0, 10000).Select( i => new PointF( (float)randomizer.NextDouble(), (float)randomizer.NextDouble())).ToList();
Stopwatch sw = Stopwatch.StartNew();
var r1 = Closest_BruteForce(points);
sw.Stop();
Debugger.Log(1, "", string.Format("Time used (Brute force) (float): {0} ms", sw.Elapsed.TotalMilliseconds));
Stopwatch sw2 = Stopwatch.StartNew();
var result2 = Closest_Recursive(points);
sw2.Stop();
Debugger.Log(1, "", string.Format("Time used (Divide & Conquer): {0} ms",sw2.Elapsed.TotalMilliseconds));
Assert.Equal(r1.Length(), result2.Length());
Output:
Time used (Brute force) (float): 145731.8935 ms
Time used (Divide & Conquer): 1139.2111 ms

Non Linq Brute Force:

 
Segment Closest_BruteForce(List<PointF> points)
{
Trace.Assert(points.Count >= 2);
 
int count = points.Count;
 
// Seed the result - doesn't matter what points are used
// This just avoids having to do null checks in the main loop below
var result = new Segment(points[0], points[1]);
var bestLength = result.Length();
 
for (int i = 0; i < count; i++)
for (int j = i + 1; j < count; j++)
if (Segment.Length(points[i], points[j]) < bestLength)
{
result = new Segment(points[i], points[j]);
bestLength = result.Length();
}
 
return result;
}

Targeted Search: Much simpler than divide and conquer, and actually runs faster for the random points. Key optimization is that if the distance along the X axis is greater than the best total length you already have, you can terminate the inner loop early. However, as only sorts in the X direction, it degenerates into an N^2 algorithm if all the points have the same X.

 
Segment Closest(List<PointF> points)
{
Trace.Assert(points.Count >= 2);
 
int count = points.Count;
points.Sort((lhs, rhs) => lhs.X.CompareTo(rhs.X));
 
var result = new Segment(points[0], points[1]);
var bestLength = result.Length();
 
for (int i = 0; i < count; i++)
{
var from = points[i];
 
for (int j = i + 1; j < count; j++)
{
var to = points[j];
 
var dx = to.X - from.X;
if (dx >= bestLength)
{
break;
}
 
if (Segment.Length(from, to) < bestLength)
{
result = new Segment(from, to);
bestLength = result.Length();
}
}
}
 
return result;
}
 

[edit] D

[edit] Compact Versions

import std.stdio, std.typecons, std.math, std.algorithm,
std.random, std.traits, std.range, std.complex;
 
auto bruteForceClosestPair(T)(in T[] points) pure nothrow @nogc {
// return pairwise(points.length.iota, points.length.iota)
// .reduce!(min!((i, j) => abs(points[i] - points[j])));
auto minD = Unqual!(typeof(T.re)).infinity;
T minI, minJ;
foreach (immutable i, const p1; points.dropBackOne)
foreach (const p2; points[i + 1 .. $]) {
immutable dist = abs(p1 - p2);
if (dist < minD) {
minD = dist;
minI = p1;
minJ = p2;
}
}
return tuple(minD, minI, minJ);
}
 
auto closestPair(T)(T[] points) pure nothrow {
static Tuple!(typeof(T.re), T, T) inner(in T[] xP, /*in*/ T[] yP)
pure nothrow {
if (xP.length <= 3)
return xP.bruteForceClosestPair;
const Pl = xP[0 .. $ / 2];
const Pr = xP[$ / 2 .. $];
immutable xDiv = Pl.back.re;
auto Yr = yP.partition!(p => p.re <= xDiv);
immutable dl_pairl = inner(Pl, yP[0 .. yP.length - Yr.length]);
immutable dr_pairr = inner(Pr, Yr);
immutable dm_pairm = dl_pairl[0]<dr_pairr[0] ? dl_pairl : dr_pairr;
immutable dm = dm_pairm[0];
const nextY = yP.filter!(p => abs(p.re - xDiv) < dm).array;
 
if (nextY.length > 1) {
auto minD = typeof(T.re).infinity;
size_t minI, minJ;
foreach (immutable i; 0 .. nextY.length - 1)
foreach (immutable j; i + 1 .. min(i + 8, nextY.length)) {
immutable double dist = abs(nextY[i] - nextY[j]);
if (dist < minD) {
minD = dist;
minI = i;
minJ = j;
}
}
return dm <= minD ? dm_pairm :
typeof(return)(minD, nextY[minI], nextY[minJ]);
} else
return dm_pairm;
}
 
points.sort!q{ a.re < b.re };
const xP = points.dup;
points.sort!q{ a.im < b.im };
return inner(xP, points);
}
 
void main() {
alias C = complex;
auto pts = [C(5,9), C(9,3), C(2), C(8,4), C(7,4), C(9,10), C(1,9),
C(8,2), C(0,10), C(9,6)];
pts.writeln;
writeln("bruteForceClosestPair: ", pts.bruteForceClosestPair);
writeln(" closestPair: ", pts.closestPair);
 
rndGen.seed = 1;
Complex!double[10_000] points;
foreach (ref p; points)
p = C(uniform(0.0, 1000.0) + uniform(0.0, 1000.0));
writeln("bruteForceClosestPair: ", points.bruteForceClosestPair);
writeln(" closestPair: ", points.closestPair);
}
Output:
[5+9i, 9+3i, 2+0i, 8+4i, 7+4i, 9+10i, 1+9i, 8+2i, 0+10i, 9+6i]
bruteForceClosestPair: Tuple!(double, Complex!double, Complex!double)(1, 8+4i, 7+4i)
          closestPair: Tuple!(double, Complex!double, Complex!double)(1, 7+4i, 8+4i)
bruteForceClosestPair: Tuple!(double, Complex!double, Complex!double)(1.76951e-05, 1040.2+0i, 1040.2+0i)
          closestPair: Tuple!(double, Complex!double, Complex!double)(1.76951e-05, 1040.2+0i, 1040.2+0i)

About 1.87 seconds run-time for data generation and brute force version, and about 0.03 seconds for data generation and divide & conquer (10_000 points in both cases) with ldc2 compiler.

[edit] Faster Brute-force Version

import std.stdio, std.random, std.math, std.typecons, std.complex,
std.traits;
 
Nullable!(Tuple!(size_t, size_t))
bfClosestPair2(T)(in Complex!T[] points) pure nothrow @nogc {
auto minD = Unqual!(typeof(points[0].re)).infinity;
if (points.length < 2)
return typeof(return)();
 
size_t minI, minJ;
foreach (immutable i; 0 .. points.length - 1)
foreach (immutable j; i + 1 .. points.length) {
auto dist = (points[i].re - points[j].re) ^^ 2;
if (dist < minD) {
dist += (points[i].im - points[j].im) ^^ 2;
if (dist < minD) {
minD = dist;
minI = i;
minJ = j;
}
}
}
 
return typeof(return)(tuple(minI, minJ));
}
 
void main() {
alias C = Complex!double;
auto rng = 31415.Xorshift;
C[10_000] pts;
foreach (ref p; pts)
p = C(uniform(0.0, 1000.0, rng), uniform(0.0, 1000.0, rng));
 
immutable ij = pts.bfClosestPair2;
if (ij.isNull)
return;
writefln("Closest pair: Distance: %f p1, p2: %f, %f",
abs(pts[ij[0]] - pts[ij[1]]), pts[ij[0]], pts[ij[1]]);
}
Output:
Closest pair: Distance: 0.019212  p1, p2: 9.74223+119.419i, 9.72306+119.418i

About 0.12 seconds run-time for brute-force version 2 (10_000 points) with with LDC2 compiler.

[edit] F#

Brute force:

 
let closest_pairs (xys: Point []) =
let n = xys.Length
seq { for i in 0..n-2 do
for j in i+1..n-1 do
yield xys.[i], xys.[j] }
|> Seq.minBy (fun (p0, p1) -> (p1 - p0).LengthSquared)
 

For example:

 
closest_pairs
[|Point(0.0, 0.0); Point(1.0, 0.0); Point (2.0, 2.0)|]
 

gives:

 
(0,0, 1,0)
 

Divide And Conquer:

 
 
open System;
open System.Drawing;
open System.Diagnostics;
 
let Length (seg : (PointF * PointF) option) =
match seg with
| None -> System.Single.MaxValue
| Some(line) ->
let f = fst line
let t = snd line
 
let dx = f.X - t.X
let dy = f.Y - t.Y
sqrt (dx*dx + dy*dy)
 
 
let Shortest a b =
if Length(a) < Length(b) then
a
else
b
 
 
let rec ClosestBoundY from maxY (ptsByY : PointF list) =
match ptsByY with
| [] -> None
| hd :: tl ->
if hd.Y > maxY then
None
else
let toHd = Some(from, hd)
let bestToRest = ClosestBoundY from maxY tl
Shortest toHd bestToRest
 
 
let rec ClosestWithinRange ptsByY maxDy =
match ptsByY with
| [] -> None
| hd :: tl ->
let fromHd = ClosestBoundY hd (hd.Y + maxDy) tl
let fromRest = ClosestWithinRange tl maxDy
Shortest fromHd fromRest
 
 
// Cuts pts half way through it's length
// Order is not maintained in result lists however
let Halve pts =
let rec ShiftToFirst first second n =
match (n, second) with
| 0, _ -> (first, second) // finished the split, so return current state
| _, [] -> (first, []) // not enough items, so first takes the whole original list
| n, hd::tl -> ShiftToFirst (hd :: first) tl (n-1) // shift 1st item from second to first, then recurse with n-1
 
let n = (List.length pts) / 2
ShiftToFirst [] pts n
 
 
let rec ClosestPair (pts : PointF list) =
if List.length pts < 2 then
None
else
let ptsByX = pts |> List.sortBy(fun(p) -> p.X)
 
let (left, right) = Halve ptsByX
let leftResult = ClosestPair left
let rightResult = ClosestPair right
 
let bestInHalf = Shortest leftResult rightResult
let bestLength = Length bestInHalf
 
let divideX = List.head(right).X
let inBand = pts |> List.filter(fun(p) -> Math.Abs(p.X - divideX) < bestLength)
 
let byY = inBand |> List.sortBy(fun(p) -> p.Y)
let bestCross = ClosestWithinRange byY bestLength
Shortest bestInHalf bestCross
 
 
let GeneratePoints n =
let rand = new Random()
[1..n] |> List.map(fun(i) -> new PointF(float32(rand.NextDouble()), float32(rand.NextDouble())))
 
let timer = Stopwatch.StartNew()
let pts = GeneratePoints (50 * 1000)
let closest = ClosestPair pts
let takenMs = timer.ElapsedMilliseconds
 
printfn "Closest Pair '%A'. Distance %f" closest (Length closest)
printfn "Took %d [ms]" takenMs
 

[edit] Fantom

(Based on the Ruby example.)

 
class Point
{
Float x
Float y
 
// create a random point
new make (Float x := Float.random * 10, Float y := Float.random * 10)
{
this.x = x
this.y = y
}
 
Float distance (Point p)
{
((x-p.x)*(x-p.x) + (y-p.y)*(y-p.y)).sqrt
}
 
override Str toStr () { "($x, $y)" }
}
 
class Main
{
// use brute force approach
static Point[] findClosestPair1 (Point[] points)
{
if (points.size < 2) return points // list too small
Point[] closestPair := [points[0], points[1]]
Float closestDistance := points[0].distance(points[1])
 
(1..<points.size).each |Int i|
{
((i+1)..<points.size).each |Int j|
{
Float trydistance := points[i].distance(points[j])
if (trydistance < closestDistance)
{
closestPair = [points[i], points[j]]
closestDistance = trydistance
}
}
}
 
return closestPair
}
 
// use recursive divide-and-conquer approach
static Point[] findClosestPair2 (Point[] points)
{
if (points.size <= 3) return findClosestPair1(points)
points.sort |Point a, Point b -> Int| { a.x <=> b.x }
bestLeft := findClosestPair2 (points[0..(points.size/2)])
bestRight := findClosestPair2 (points[(points.size/2)..-1])
 
Float minDistance
Point[] closePoints := [,]
if (bestLeft[0].distance(bestLeft[1]) < bestRight[0].distance(bestRight[1]))
{
minDistance = bestLeft[0].distance(bestLeft[1])
closePoints = bestLeft
}
else
{
minDistance = bestRight[0].distance(bestRight[1])
closePoints = bestRight
}
yPoints := points.findAll |Point p -> Bool|
{
(points.last.x - p.x).abs < minDistance
}.sort |Point a, Point b -> Int| { a.y <=> b.y }
 
closestPair := [,]
closestDist := Float.posInf
 
for (Int i := 0; i < yPoints.size - 1; ++i)
{
for (Int j := (i+1); j < yPoints.size; ++j)
{
if ((yPoints[j].y - yPoints[i].y) >= minDistance)
{
break
}
else
{
dist := yPoints[i].distance (yPoints[j])
if (dist < closestDist)
{
closestDist = dist
closestPair = [yPoints[i], yPoints[j]]
}
}
}
}
if (closestDist < minDistance)
return closestPair
else
return closePoints
}
 
public static Void main (Str[] args)
{
Int numPoints := 10 // default value, in case a number not given on command line
if ((args.size > 0) && (args[0].toInt(10, false) != null))
{
numPoints = args[0].toInt(10, false)
}
 
Point[] points := [,]
numPoints.times { points.add (Point()) }
 
Int t1 := Duration.now.toMillis
echo (findClosestPair1(points.dup))
Int t2 := Duration.now.toMillis
echo ("Time taken: ${(t2-t1)}ms")
echo (findClosestPair2(points.dup))
Int t3 := Duration.now.toMillis
echo ("Time taken: ${(t3-t2)}ms")
}
}
 
Output:
$ fan closestPoints 1000
[(1.4542885676006445, 8.238581003965352), (1.4528464044751888, 8.234724407229772)]
Time taken: 88ms
[(1.4528464044751888, 8.234724407229772), (1.4542885676006445, 8.238581003965352)]
Time taken: 80ms
$ fan closestPoints 10000
[(3.454790171891945, 5.307252398266497), (3.4540208686702245, 5.308350223433488)]
Time taken: 6248ms
[(3.454790171891945, 5.307252398266497), (3.4540208686702245, 5.308350223433488)]
Time taken: 228ms

[edit] Fortran

See Closest pair problem/Fortran

[edit] Go

Brute force

package main
 
import (
"fmt"
"math"
"math/rand"
)
 
type xy struct {
x, y float64
}
 
const n = 1000
const scale = 1.
 
func d(p1, p2 xy) float64 {
dx := p2.x - p1.x
dy := p2.y - p1.y
return math.Sqrt(dx*dx + dy*dy)
}
 
func main() {
points := make([]xy, n)
for i := range points {
points[i] = xy{rand.Float64(), rand.Float64() * scale}
}
p1, p2 := closestPair(points)
fmt.Println(p1, p2)
fmt.Println("distance:", d(p1, p2))
}
 
func closestPair(points []xy) (p1, p2 xy) {
if len(points) < 2 {
panic("at least two points expected")
}
min := 2 * scale
for i, q1 := range points[:len(points)-1] {
for _, q2 := range points[i+1:] {
if dq := d(q1, q2); dq < min {
p1, p2 = q1, q2
min = dq
}
}
}
return
}

O(n)

// implementation following algorithm described in
// http://www.cs.umd.edu/~samir/grant/cp.pdf
package main
 
import (
"fmt"
"math"
"math/rand"
)
 
// number of points to search for closest pair
const n = 1e6
 
// size of bounding box for points.
// x and y will be random with uniform distribution in the range [0,scale).
const scale = 1.
 
// point struct
type xy struct {
x, y float64 // coordinates
key int64 // an annotation used in the algorithm
}
 
// Euclidian distance
func d(p1, p2 xy) float64 {
dx := p2.x - p1.x
dy := p2.y - p1.y
return math.Sqrt(dx*dx + dy*dy)
}
 
func main() {
points := make([]xy, n)
for i := range points {
points[i] = xy{rand.Float64(), rand.Float64() * scale, 0}
}
p1, p2 := closestPair(points)
fmt.Println(p1, p2)
fmt.Println("distance:", d(p1, p2))
}
 
func closestPair(s []xy) (p1, p2 xy) {
if len(s) < 2 {
panic("2 points required")
}
var dxi float64
// step 0
for s1, i := s, 1; ; i++ {
// step 1: compute min distance to a random point
// (for the case of random data, it's enough to just try
// to pick a different point)
rp := i % len(s1)
xi := s1[rp]
dxi = 2 * scale
for p, xn := range s1 {
if p != rp {
if dq := d(xi, xn); dq < dxi {
dxi = dq
}
}
}
 
// step 2: filter
invB := 3 / dxi // b is size of a mesh cell
mx := int64(scale*invB) + 1 // mx is number of cells along a side
// construct map as a histogram:
// key is index into mesh. value is count of points in cell
hm := make(map[int64]int)
for ip, p := range s1 {
key := int64(p.x*invB)*mx + int64(p.y*invB)
s1[ip].key = key
hm[key]++
}
// construct s2 = s1 less the points without neighbors
var s2 []xy
nx := []int64{-mx - 1, -mx, -mx + 1, -1, 0, 1, mx - 1, mx, mx + 1}
for i, p := range s1 {
nn := 0
for _, ofs := range nx {
nn += hm[p.key+ofs]
if nn > 1 {
s2 = append(s2, s1[i])
break
}
}
}
 
// step 3: done?
if len(s2) == 0 {
break
}
s1 = s2
}
// step 4: compute answer from approximation
invB := 1 / dxi
mx := int64(scale*invB) + 1
hm := make(map[int64][]int)
for i, p := range s {
key := int64(p.x*invB)*mx + int64(p.y*invB)
s[i].key = key
hm[key] = append(hm[key], i)
}
nx := []int64{-mx - 1, -mx, -mx + 1, -1, 0, 1, mx - 1, mx, mx + 1}
var min = scale * 2
for ip, p := range s {
for _, ofs := range nx {
for _, iq := range hm[p.key+ofs] {
if ip != iq {
if d1 := d(p, s[iq]); d1 < min {
min = d1
p1, p2 = p, s[iq]
}
}
}
}
}
return p1, p2
}

[edit] Groovy

Point class:

class Point {
final Number x, y
Point(Number x = 0, Number y = 0) { this.x = x; this.y = y }
Number distance(Point that) { ((this.x - that.x)**2 + (this.y - that.y)**2)**0.5 }
String toString() { "{x:${x}, y:${y}}" }
}

Brute force solution. Incorporates X-only and Y-only pre-checks in two places to cut down on the square root calculations:

def bruteClosest(Collection pointCol) {
assert pointCol
List l = pointCol
int n = l.size()
assert n > 1
if (n == 2) return [distance:l[0].distance(l[1]), points:[l[0],l[1]]]
def answer = [distance: Double.POSITIVE_INFINITY]
(0..<(n-1)).each { i ->
((i+1)..<n).findAll { j ->
(l[i].x - l[j].x).abs() < answer.distance &&
(l[i].y - l[j].y).abs() < answer.distance
}.each { j ->
if ((l[i].x - l[j].x).abs() < answer.distance &&
(l[i].y - l[j].y).abs() < answer.distance) {
def dist = l[i].distance(l[j])
if (dist < answer.distance) {
answer = [distance:dist, points:[l[i],l[j]]]
}
}
}
}
answer
}

Elegant (divide-and-conquer reduction) solution. Incorporates X-only and Y-only pre-checks in two places (four if you count the inclusion of the brute force solution) to cut down on the square root calculations:

def elegantClosest(Collection pointCol) {
assert pointCol
List xList = (pointCol as List).sort { it.x }
List yList = xList.clone().sort { it.y }
reductionClosest(xList, xList)
}
 
def reductionClosest(List xPoints, List yPoints) {
// assert xPoints && yPoints
// assert (xPoints as Set) == (yPoints as Set)
int n = xPoints.size()
if (n < 10) return bruteClosest(xPoints)
 
int nMid = Math.ceil(n/2)
List xLeft = xPoints[0..<nMid]
List xRight = xPoints[nMid..<n]
Number xMid = xLeft[-1].x
List yLeft = yPoints.findAll { it.x <= xMid }
List yRight = yPoints.findAll { it.x > xMid }
if (xRight[0].x == xMid) {
yLeft = xLeft.collect{ it }.sort { it.y }
yRight = xRight.collect{ it }.sort { it.y }
}
 
Map aLeft = reductionClosest(xLeft, yLeft)
Map aRight = reductionClosest(xRight, yRight)
Map aMin = aRight.distance < aLeft.distance ? aRight : aLeft
List yMid = yPoints.findAll { (xMid - it.x).abs() < aMin.distance }
int nyMid = yMid.size()
if (nyMid < 2) return aMin
 
Map answer = aMin
(0..<(nyMid-1)).each { i ->
((i+1)..<nyMid).findAll { j ->
(yMid[j].x - yMid[i].x).abs() < aMin.distance &&
(yMid[j].y - yMid[i].y).abs() < aMin.distance &&
yMid[j].distance(yMid[i]) < aMin.distance
}.each { k ->
if ((yMid[k].x - yMid[i].x).abs() < answer.distance && (yMid[k].y - yMid[i].y).abs() < answer.distance) {
def ikDist = yMid[i].distance(yMid[k])
if ( ikDist < answer.distance) {
answer = [distance:ikDist, points:[yMid[i],yMid[k]]]
}
}
}
}
answer
}

Benchmark/Test:

def random = new Random()
 
(1..4).each {
def point10 = (0..<(10**it)).collect { new Point(random.nextInt(1000001) - 500000,random.nextInt(1000001) - 500000) }
 
def startE = System.currentTimeMillis()
def closestE = elegantClosest(point10)
def elapsedE = System.currentTimeMillis() - startE
println """
${10**it} POINTS
-----------------------------------------
Elegant reduction:
elapsed: ${elapsedE/1000} s
closest: ${closestE}
"""

 
 
def startB = System.currentTimeMillis()
def closestB = bruteClosest(point10)
def elapsedB = System.currentTimeMillis() - startB
println """Brute force:
elapsed: ${elapsedB/1000} s
closest: ${closestB}
 
Speedup ratio (B/E): ${elapsedB/elapsedE}
=========================================
"""

}

Results:

10 POINTS
-----------------------------------------
Elegant reduction:
elapsed: 0.019 s
closest: [distance:85758.5249173515, points:[{x:310073, y:-27339}, {x:382387, y:18761}]]

Brute force:
elapsed: 0.001 s
closest: [distance:85758.5249173515, points:[{x:310073, y:-27339}, {x:382387, y:18761}]]

Speedup ratio (B/E): 0.0526315789
=========================================


100 POINTS
-----------------------------------------
Elegant reduction:
elapsed: 0.019 s
closest: [distance:3166.229934796271, points:[{x:-343735, y:-244394}, {x:-341099, y:-246148}]]

Brute force:
elapsed: 0.027 s
closest: [distance:3166.229934796271, points:[{x:-343735, y:-244394}, {x:-341099, y:-246148}]]

Speedup ratio (B/E): 1.4210526316
=========================================


1000 POINTS
-----------------------------------------
Elegant reduction:
elapsed: 0.241 s
closest: [distance:374.22586762542215, points:[{x:411817, y:-83016}, {x:412038, y:-82714}]]

Brute force:
elapsed: 0.618 s
closest: [distance:374.22586762542215, points:[{x:411817, y:-83016}, {x:412038, y:-82714}]]

Speedup ratio (B/E): 2.5643153527
=========================================


10000 POINTS
-----------------------------------------
Elegant reduction:
elapsed: 1.957 s
closest: [distance:79.00632886041473, points:[{x:187928, y:-452338}, {x:187929, y:-452259}]]

Brute force:
elapsed: 51.567 s
closest: [distance:79.00632886041473, points:[{x:187928, y:-452338}, {x:187929, y:-452259}]]

Speedup ratio (B/E): 26.3500255493
=========================================

[edit] Haskell

BF solution:

import Data.List
import System.Random
import Control.Monad
import Control.Arrow
import Data.Ord
 
vecLeng [[a,b],[p,q]] = sqrt $ (a-p)^2+(b-q)^2
 
findClosestPair = foldl1' ((minimumBy (comparing vecLeng). ). (. return). (:)) .
concatMap (\(x:xs) -> map ((x:).return) xs) . init . tails
 
testCP = do
g <- newStdGen
let pts :: [[Double]]
pts = take 1000. unfoldr (Just. splitAt 2) $ randomRs(-1,1) g
print . (id &&& vecLeng ) . findClosestPair $ pts
Output:
*Main> testCP
([[0.8347201880148426,0.40774840545089647],[0.8348731214261784,0.4087113189531284]],9.749825850154334e-4)
(4.02 secs, 488869056 bytes)

[edit] Icon and Unicon

This is a brute force solution. It combines reading the points with computing the closest pair seen so far.

record point(x,y)
 
procedure main()
minDist := 0
minPair := &null
every (points := [],p1 := readPoint()) do {
if *points == 1 then minDist := dSquared(p1,points[1])
every minDist >=:= dSquared(p1,p2 := !points) do minPair := [p1,p2]
push(points, p1)
}
 
if \minPair then {
write("(",minPair[1].x,",",minPair[1].y,") -> ",
"(",minPair[2].x,",",minPair[2].y,")")
}
else write("One or fewer points!")
end
 
procedure readPoint() # Skips lines that don't have two numbers on them
suspend !&input ? point(numeric(tab(upto(', '))), numeric((move(1),tab(0))))
end
 
procedure dSquared(p1,p2) # Compute the square of the distance
return (p2.x-p1.x)^2 + (p2.y-p1.y)^2 # (sufficient for closeness)
end

[edit] J

Solution of the simpler (brute-force) problem:

vecl   =:  +/"1&.:*:                  NB. length of each of vectors
dist =: <@:vecl@:({: -"1 }:)\ NB. calculate all distances among vectors
minpair=: ({~ > {.@($ #: I.@,)@:= <./@;)dist NB. find one pair of the closest points
closestpairbf =: (; vecl@:-/)@minpair NB. the pair and their distance

Examples of use:

   ]pts=:10 2 ?@$ 0
0.654682 0.925557
0.409382 0.619391
0.891663 0.888594
0.716629 0.9962
0.477721 0.946355
0.925092 0.81822
0.624291 0.142924
0.211332 0.221507
0.293786 0.691701
0.839186 0.72826
 
closestpairbf pts
+-----------------+---------+
|0.891663 0.888594|0.0779104|
|0.925092 0.81822| |
+-----------------+---------+

The program also works for higher dimensional vectors:

   ]pts=:10 4 ?@$ 0
0.559164 0.482993 0.876 0.429769
0.217911 0.729463 0.97227 0.132175
0.479206 0.169165 0.495302 0.362738
0.316673 0.797519 0.745821 0.0598321
0.662585 0.726389 0.658895 0.653457
0.965094 0.664519 0.084712 0.20671
0.840877 0.591713 0.630206 0.99119
0.221416 0.114238 0.0991282 0.174741
0.946262 0.505672 0.776017 0.307362
0.262482 0.540054 0.707342 0.465234
 
closestpairbf pts
+------------------------------------+--------+
|0.217911 0.729463 0.97227 0.132175|0.708555|
|0.316673 0.797519 0.745821 0.0598321| |
+------------------------------------+--------+

[edit] Java

Both the brute-force and the divide-and-conquer methods are implemented.

Code:

import java.util.*;
 
public class ClosestPair
{
public static class Point
{
public final double x;
public final double y;
 
public Point(double x, double y)
{
this.x = x;
this.y = y;
}
 
public String toString()
{ return "(" + x + ", " + y + ")"; }
}
 
public static class Pair
{
public Point point1 = null;
public Point point2 = null;
public double distance = 0.0;
 
public Pair()
{ }
 
public Pair(Point point1, Point point2)
{
this.point1 = point1;
this.point2 = point2;
calcDistance();
}
 
public void update(Point point1, Point point2, double distance)
{
this.point1 = point1;
this.point2 = point2;
this.distance = distance;
}
 
public void calcDistance()
{ this.distance = distance(point1, point2); }
 
public String toString()
{ return point1 + "-" + point2 + " : " + distance; }
}
 
public static double distance(Point p1, Point p2)
{
double xdist = p2.x - p1.x;
double ydist = p2.y - p1.y;
return Math.hypot(xdist, ydist);
}
 
public static Pair bruteForce(List<? extends Point> points)
{
int numPoints = points.size();
if (numPoints < 2)
return null;
Pair pair = new Pair(points.get(0), points.get(1));
if (numPoints > 2)
{
for (int i = 0; i < numPoints - 1; i++)
{
Point point1 = points.get(i);
for (int j = i + 1; j < numPoints; j++)
{
Point point2 = points.get(j);
double distance = distance(point1, point2);
if (distance < pair.distance)
pair.update(point1, point2, distance);
}
}
}
return pair;
}
 
public static void sortByX(List<? extends Point> points)
{
Collections.sort(points, new Comparator<Point>() {
public int compare(Point point1, Point point2)
{
if (point1.x < point2.x)
return -1;
if (point1.x > point2.x)
return 1;
return 0;
}
}
);
}
 
public static void sortByY(List<? extends Point> points)
{
Collections.sort(points, new Comparator<Point>() {
public int compare(Point point1, Point point2)
{
if (point1.y < point2.y)
return -1;
if (point1.y > point2.y)
return 1;
return 0;
}
}
);
}
 
public static Pair divideAndConquer(List<? extends Point> points)
{
List<Point> pointsSortedByX = new ArrayList<Point>(points);
sortByX(pointsSortedByX);
List<Point> pointsSortedByY = new ArrayList<Point>(points);
sortByY(pointsSortedByY);
return divideAndConquer(pointsSortedByX, pointsSortedByY);
}
 
private static Pair divideAndConquer(List<? extends Point> pointsSortedByX, List<? extends Point> pointsSortedByY)
{
int numPoints = pointsSortedByX.size();
if (numPoints <= 3)
return bruteForce(pointsSortedByX);
 
int dividingIndex = numPoints >>> 1;
List<? extends Point> leftOfCenter = pointsSortedByX.subList(0, dividingIndex);
List<? extends Point> rightOfCenter = pointsSortedByX.subList(dividingIndex, numPoints);
 
List<Point> tempList = new ArrayList<Point>(leftOfCenter);
sortByY(tempList);
Pair closestPair = divideAndConquer(leftOfCenter, tempList);
 
tempList.clear();
tempList.addAll(rightOfCenter);
sortByY(tempList);
Pair closestPairRight = divideAndConquer(rightOfCenter, tempList);
 
if (closestPairRight.distance < closestPair.distance)
closestPair = closestPairRight;
 
tempList.clear();
double shortestDistance =closestPair.distance;
double centerX = rightOfCenter.get(0).x;
for (Point point : pointsSortedByY)
if (Math.abs(centerX - point.x) < shortestDistance)
tempList.add(point);
 
for (int i = 0; i < tempList.size() - 1; i++)
{
Point point1 = tempList.get(i);
for (int j = i + 1; j < tempList.size(); j++)
{
Point point2 = tempList.get(j);
if ((point2.y - point1.y) >= shortestDistance)
break;
double distance = distance(point1, point2);
if (distance < closestPair.distance)
{
closestPair.update(point1, point2, distance);
shortestDistance = distance;
}
}
}
return closestPair;
}
 
public static void main(String[] args)
{
int numPoints = (args.length == 0) ? 1000 : Integer.parseInt(args[0]);
List<Point> points = new ArrayList<Point>();
Random r = new Random();
for (int i = 0; i < numPoints; i++)
points.add(new Point(r.nextDouble(), r.nextDouble()));
System.out.println("Generated " + numPoints + " random points");
long startTime = System.currentTimeMillis();
Pair bruteForceClosestPair = bruteForce(points);
long elapsedTime = System.currentTimeMillis() - startTime;
System.out.println("Brute force (" + elapsedTime + " ms): " + bruteForceClosestPair);
startTime = System.currentTimeMillis();
Pair dqClosestPair = divideAndConquer(points);
elapsedTime = System.currentTimeMillis() - startTime;
System.out.println("Divide and conquer (" + elapsedTime + " ms): " + dqClosestPair);
if (bruteForceClosestPair.distance != dqClosestPair.distance)
System.out.println("MISMATCH");
}
}
Output:
java ClosestPair 10000
Generated 10000 random points
Brute force (1594 ms): (0.9246533850872104, 0.098709007587097)-(0.924591196030625, 0.09862206991823985) : 1.0689077146927108E-4
Divide and conquer (250 ms): (0.924591196030625, 0.09862206991823985)-(0.9246533850872104, 0.098709007587097) : 1.0689077146927108E-4

[edit] JavaScript

Using bruteforce algorithm, the bruteforceClosestPair method below expects an array of objects with x- and y-members set to numbers, and returns an object containing the members distance and points.

function distance(p1, p2) {
var dx = Math.abs(p1.x - p2.x);
var dy = Math.abs(p1.y - p2.y);
return Math.sqrt(dx*dx + dy*dy);
}
 
function bruteforceClosestPair(arr) {
if (arr.length < 2) {
return Infinity;
} else {
var minDist = distance(arr[0], arr[1]);
var minPoints = arr.slice(0, 2);
 
for (var i=0; i<arr.length-1; i++) {
for (var j=i+1; j<arr.length; j++) {
if (distance(arr[i], arr[j]) < minDist) {
minDist = distance(arr[i], arr[j]);
minPoints = [ arr[i], arr[j] ];
}
}
}
return {
distance: minDist,
points: minPoints
};
}
}

[edit] jq

Works with: jq version 1.4

The solution presented here is essentially a direct translation into jq of the pseudo-code presented in the task description, but "closest_pair" is added so that any list of [x,y] points can be presented, and extra lines are added to ensure that xL and yL have the same lengths.

Infrastructure:

# Iteratively evaluate "next" while the condition is satisfied;
# emit the last value that satisfies the condition, or else the original input.
def dowhile(condition; next):
def w: if condition then ., (next|w) else empty end;
def last(f): reduce f as $i (.; $i);
last(w);
 
# Euclidean 2d distance
def dist(x;y):
[x[0] - y[0], x[1] - y[1]] | map(.*.) | add | sqrt;
 
# P is an array of points, [x,y].
# Emit the solution in the form [dist, [P1, P2]]
def bruteForceClosestPair(P):
(P|length) as $length
| if $length < 2 then null
else
reduce range(0; $length-1) as $i
( null;
reduce range($i+1; $length) as $j
(.;
dist(P[$i]; P[$j]) as $d
| if . == null or $d < .[0] then [$d, [ P[$i], P[$j] ] ] else . end ) )
end;
 
def closest_pair:
 
def abs: if . < 0 then -. else . end;
def ceil: floor as $floor
| if . == $floor then $floor else $floor + 1 end;
 
# xP is an array [P(1), .. P(N)] sorted by x coordinate, and
# yP is an array [P(1), .. P(N)] sorted by y coordinate (ascending order).
# if N <= 3 then return closest points of xP using the brute-force algorithm.
def closestPair(xP; yP):
if xP|length <= 3 then bruteForceClosestPair(xP)
else
((xP|length)/2|ceil) as $N
| xP[0:$N] as $xL
| xP[$N:] as $xR
| xP[$N-1][0] as $xm # middle
| (yP | map(select(.[0] <= $xm ))) as $yL0 # might be too long
| (yP | map(select(.[0] > $xm ))) as $yR0 # might be too short
| (if $yL0|length == $N then $yL0 else $yL0[0:$N] end) as $yL
| (if $yL0|length == $N then $yR0 else $yL0[$N:] + $yR0 end) as $yR
| closestPair($xL; $yL) as $pairL # [dL, pairL]
| closestPair($xR; $yR) as $pairR # [dR, pairR]
| (if $pairL[0] < $pairR[0] then $pairL else $pairR end) as $pair # [ dmin, pairMin]
| (yP | map(select( (($xm - .[0])|abs) < $pair[0]))) as $yS
| ($yS | length) as $nS
| $pair[0] as $dmin
| reduce range(0; $nS - 1) as $i
( [0, $pair]; # state: [k, [d, [P1,P2]]]
.[0] = $i + 1
| dowhile( .[0] as $k | $k < $nS and ($yS[$k][1] - $yS[$i][1]) < $dmin;
.[0] as $k
| dist($yS[$k]; $yS[$i]) as $d
| if $d < .[1][0]
then [$k+1, [ $d, [$yS[$k], $yS[$i]]]]
else .[0] += 1
end) )
| .[1]
end;
closestPair( sort_by(.[0]); sort_by(.[1])) ;

Example from the Mathematica section:

def data:
[[0.748501, 4.09624],
[3.00302, 5.26164],
[3.61878, 9.52232],
[7.46911, 4.71611],
[5.7819, 2.69367],
[2.34709, 8.74782],
[2.87169, 5.97774],
[6.33101, 0.463131],
[7.46489, 4.6268],
[1.45428, 0.087596] ];
 
data | closest_pair
Output:
$jq -M -c -n -f closest_pair.jq
[0.0894096443343775,[[7.46489,4.6268],[7.46911,4.71611]]]

[edit] Liberty BASIC

NB array terms can not be READ directly.

 
N =10
 
dim x( N), y( N)
 
firstPt =0
secondPt =0
 
for i =1 to N
read f: x( i) =f
read f: y( i) =f
next i
 
minDistance =1E6
 
for i =1 to N -1
for j =i +1 to N
dxSq =( x( i) -x( j))^2
dySq =( y( i) -y( j))^2
D =abs( ( dxSq +dySq)^0.5)
if D <minDistance then
minDistance =D
firstPt =i
secondPt =j
end if
next j
next i
 
print "Distance ="; minDistance; " between ( "; x( firstPt); ", "; y( firstPt); ") and ( "; x( secondPt); ", "; y( secondPt); ")"
 
end
 
data 0.654682, 0.925557
data 0.409382, 0.619391
data 0.891663, 0.888594
data 0.716629, 0.996200
data 0.477721, 0.946355
data 0.925092, 0.818220
data 0.624291, 0.142924
data 0.211332, 0.221507
data 0.293786, 0.691701
data 0.839186, 0.72826
 
 
Distance =0.77910191e-1 between ( 0.891663, 0.888594) and ( 0.925092, 0.81822)

[edit] Mathematica

nearestPair[data_] := 
Block[{pos, dist = N[Outer[EuclideanDistance, data, data, 1]]},
pos = Position[dist, Min[DeleteCases[Flatten[dist], 0.]]];
data[[pos[[1]]]]]
Output:
nearestPair[{{0.748501, 4.09624}, {3.00302, 5.26164}, {3.61878, 
  9.52232}, {7.46911, 4.71611}, {5.7819, 2.69367}, {2.34709, 
  8.74782}, {2.87169, 5.97774}, {6.33101, 0.463131}, {7.46489, 
  4.6268}, {1.45428, 0.087596}}]

{{7.46911, 4.71611}, {7.46489, 4.6268}}

[edit] MATLAB

This solution is an almost direct translation of the above pseudo-code into MATLAB.

function [closest,closestpair] = closestPair(xP,yP)
 
N = numel(xP);
 
if(N <= 3)
 
%Brute force closestpair
if(N < 2)
closest = +Inf;
closestpair = {};
else
closest = norm(xP{1}-xP{2});
closestpair = {xP{1},xP{2}};
 
for i = ( 1:N-1 )
for j = ( (i+1):N )
if ( norm(xP{i} - xP{j}) < closest )
closest = norm(xP{i}-xP{j});
closestpair = {xP{i},xP{j}};
end %if
end %for
end %for
end %if (N < 2)
else
 
halfN = ceil(N/2);
 
xL = { xP{1:halfN} };
xR = { xP{halfN+1:N} };
xm = xP{halfN}(1);
 
%cellfun( @(p)le(p(1),xm),yP ) is the same as { p ∈ yP : px ≤ xm }
yLIndicies = cellfun( @(p)le(p(1),xm),yP );
 
yL = { yP{yLIndicies} };
yR = { yP{~yLIndicies} };
 
[dL,pairL] = closestPair(xL,yL);
[dR,pairR] = closestPair(xR,yR);
 
if dL < dR
dmin = dL;
pairMin = pairL;
else
dmin = dR;
pairMin = pairR;
end
 
%cellfun( @(p)lt(norm(xm-p(1)),dmin),yP ) is the same as
%{ p ∈ yP : |xm - px| < dmin }
yS = {yP{ cellfun( @(p)lt(norm(xm-p(1)),dmin),yP ) }};
nS = numel(yS);
 
closest = dmin;
closestpair = pairMin;
 
for i = (1:nS-1)
k = i+1;
 
while( (k<=nS) && (yS{k}(2)-yS{i}(2) < dmin) )
 
if norm(yS{k}-yS{i}) < closest
closest = norm(yS{k}-yS{i});
closestpair = {yS{k},yS{i}};
end
 
k = k+1;
end %while
end %for
end %if (N <= 3)
end %closestPair
Output:
[distance,pair]=closestPair({[0 -.3],[1 1],[1.5 2],[2 2],[3 3]},{[0 -.3],[1 1],[1.5 2],[2 2],[3 3]})
 
distance =
 
0.500000000000000
 
 
pair =
 
[1x2 double] [1x2 double] %The pair is [1.5 2] and [2 2] which is correct

[edit] Objective-C

See Closest-pair problem/Objective-C

[edit] OCaml

 
 
type point = { x : float; y : float }
 
 
let cmpPointX (a : point) (b : point) = compare a.x b.x
let cmpPointY (a : point) (b : point) = compare a.y b.y
 
 
let distSqrd (seg : (point * point) option) =
match seg with
| None -> max_float
| Some(line) ->
let a = fst line in
let b = snd line in
 
let dx = a.x -. b.x in
let dy = a.y -. b.y in
 
dx*.dx +. dy*.dy
 
 
let dist seg =
sqrt (distSqrd seg)
 
 
let shortest l1 l2 =
if distSqrd l1 < distSqrd l2 then
l1
else
l2
 
 
let halve l =
let n = List.length l in
BatList.split_at (n/2) l
 
 
let rec closestBoundY from maxY (ptsByY : point list) =
match ptsByY with
| [] -> None
| hd :: tl ->
if hd.y > maxY then
None
else
let toHd = Some(from, hd) in
let bestToRest = closestBoundY from maxY tl in
shortest toHd bestToRest
 
 
let rec closestInRange ptsByY maxDy =
match ptsByY with
| [] -> None
| hd :: tl ->
let fromHd = closestBoundY hd (hd.y +. maxDy) tl in
let fromRest = closestInRange tl maxDy in
shortest fromHd fromRest
 
 
let rec closestPairByX (ptsByX : point list) =
if List.length ptsByX < 2 then
None
else
let (left, right) = halve ptsByX in
let leftResult = closestPairByX left in
let rightResult = closestPairByX right in
 
let bestInHalf = shortest leftResult rightResult in
let bestLength = dist bestInHalf in
 
let divideX = (List.hd right).x in
let inBand = List.filter(fun(p) -> abs_float(p.x -. divideX) < bestLength) ptsByX in
 
let byY = List.sort cmpPointY inBand in
let bestCross = closestInRange byY bestLength in
shortest bestInHalf bestCross
 
 
let closestPair pts =
let ptsByX = List.sort cmpPointX pts in
closestPairByX ptsByX
 
 
let parsePoint str =
let sep = Str.regexp_string "," in
let tokens = Str.split sep str in
let xStr = List.nth tokens 0 in
let yStr = List.nth tokens 1 in
 
let xVal = (float_of_string xStr) in
let yVal = (float_of_string yStr) in
 
{ x = xVal; y = yVal }
 
 
let loadPoints filename =
let ic = open_in filename in
let result = ref [] in
try
while true do
let s = input_line ic in
if s <> "" then
let p = parsePoint s in
result := p :: !result;
done;
!result
with End_of_file ->
close_in ic;
!result
;;
 
let loaded = (loadPoints "Points.txt") in
let start = Sys.time() in
let c = closestPair loaded in
let taken = Sys.time() -. start in
Printf.printf "Took %f [s]\n" taken;
 
match c with
| None -> Printf.printf "No closest pair\n"
| Some(seg) ->
let a = fst seg in
let b = snd seg in
 
Printf.printf "(%f, %f) (%f, %f) Dist %f\n" a.x a.y b.x b.y (dist c)
 
 

[edit] Oz

Translation of pseudocode:

declare
fun {Distance X1#Y1 X2#Y2}
{Sqrt {Pow X2-X1 2.0} + {Pow Y2-Y1 2.0}}
end
 
%% brute force
fun {BFClosestPair Points=P1|P2|_}
Ps = {List.toTuple unit Points} %% for efficient random access
N = {Width Ps}
MinDist = {NewCell {Distance P1 P2}}
MinPoints = {NewCell P1#P2}
in
for I in 1..N-1 do
for J in I+1..N do
IJDist = {Distance Ps.I Ps.J}
in
if IJDist < @MinDist then
MinDist := IJDist
MinPoints := Ps.I#Ps.J
end
end
end
@MinPoints
end
 
%% divide and conquer
fun {ClosestPair Points}
case {ClosestPair2
{Sort Points {LessThanBy X}}
{Sort Points {LessThanBy Y}}}
of Distance#Pair then
Pair
end
end
 
%% XP: points sorted by X, YP: sorted by Y
%% returns a pair Distance#Pair
fun {ClosestPair2 XP YP}
N = {Length XP} = {Length YP}
in
if N =< 3 then
P = {BFClosestPair XP}
in
{Distance P.1 P.2}#P
else
XL XR
{List.takeDrop XP (N div 2) ?XL ?XR}
XM = {Nth XP (N div 2)}.X
YL YR
{List.partition YP fun {$ P} P.X =< XM end ?YL ?YR}
DL#PairL = {ClosestPair2 XL YL}
DR#PairR = {ClosestPair2 XR YR}
DMin#PairMin = if DL < DR then DL#PairL else DR#PairR end
YSList = {Filter YP fun {$ P} {Abs XM-P.X} < DMin end}
YS = {List.toTuple unit YSList} %% for efficient random access
NS = {Width YS}
Closest = {NewCell DMin}
ClosestPair = {NewCell PairMin}
in
for I in 1..NS-1 do
for K in I+1..NS while:YS.K.Y - YS.I.Y < DMin do
DistKI = {Distance YS.K YS.I}
in
if DistKI < @Closest then
Closest := DistKI
ClosestPair := YS.K#YS.I
end
end
end
@Closest#@ClosestPair
end
end
 
%% To access components when points are represented as pairs
X = 1
Y = 2
 
%% returns a less-than predicate that accesses feature F
fun {LessThanBy F}
fun {$ A B}
A.F < B.F
end
end
 
fun {Random Min Max}
Min +
{Int.toFloat {OS.rand}} * (Max-Min)
/ {Int.toFloat {OS.randLimits _}}
end
 
fun {RandomPoint}
{Random 0.0 100.0}#{Random 0.0 100.0}
end
 
Points = {MakeList 5}
in
{ForAll Points RandomPoint}
{Show Points}
{Show {ClosestPair Points}}

[edit] PARI/GP

Naive quadratic solution.

closestPair(v)={
my(r=norml2(v[1]-v[2]),at=[1,2]);
for(a=1,#v-1,
for(b=a+1,#v,
if(norml2(v[a]-v[b])<r,
at=[a,b];
r=norml2(v[a]-v[b])
)
)
);
[v[at[1]],v[at[2]]]
};

[edit] Perl

#! /usr/bin/perl
use strict;
use POSIX qw(ceil);
 
sub dist
{
my ( $a, $b) = @_;
return sqrt( ($a->[0] - $b->[0])**2 +
($a->[1] - $b->[1])**2 );
}
 
sub closest_pair_simple
{
my $ra = shift;
my @arr = @$ra;
my $inf = 1e600;
return $inf if scalar(@arr) < 2;
my ( $a, $b, $d ) = ($arr[0], $arr[1], dist($arr[0], $arr[1]));
while( @arr ) {
my $p = pop @arr;
foreach my $l (@arr) {
my $t = dist($p, $l);
($a, $b, $d) = ($p, $l, $t) if $t < $d;
}
}
return ($a, $b, $d);
}
 
sub closest_pair
{
my $r = shift;
my @ax = sort { $a->[0] <=> $b->[0] } @$r;
my @ay = sort { $a->[1] <=> $b->[1] } @$r;
return closest_pair_real(\@ax, \@ay);
}
 
sub closest_pair_real
{
my ($rx, $ry) = @_;
my @xP = @$rx;
my @yP = @$ry;
my $N = @xP;
return closest_pair_simple($rx) if scalar(@xP) <= 3;
 
my $inf = 1e600;
my $midx = ceil($N/2)-1;
 
my @PL = @xP[0 .. $midx];
my @PR = @xP[$midx+1 .. $N-1];
 
my $xm = ${$xP[$midx]}[0];
 
my @yR = ();
my @yL = ();
foreach my $p (@yP) {
if ( ${$p}[0] <= $xm ) {
push @yR, $p;
} else {
push @yL, $p;
}
}
 
my ($al, $bl, $dL) = closest_pair_real(\@PL, \@yR);
my ($ar, $br, $dR) = closest_pair_real(\@PR, \@yL);
 
my ($m1, $m2, $dmin) = ($al, $bl, $dL);
($m1, $m2, $dmin) = ($ar, $br, $dR) if $dR < $dL;
 
my @yS = ();
foreach my $p (@yP) {
push @yS, $p if abs($xm - ${$p}[0]) < $dmin;
}
 
if ( @yS ) {
my ( $w1, $w2, $closest ) = ($m1, $m2, $dmin);
foreach my $i (0 .. ($#yS - 1)) {
 
my $k = $i + 1;
while ( ($k <= $#yS) && ( (${$yS[$k]}[1] - ${$yS[$i]}[1]) < $dmin) ) {
my $d = dist($yS[$k], $yS[$i]);
($w1, $w2, $closest) = ($yS[$k], $yS[$i], $d) if $d < $closest;
$k++;
}
 
}
return ($w1, $w2, $closest);
 
} else {
return ($m1, $m2, $dmin);
}
}
 
 
 
my @points = ();
my $N = 5000;
 
foreach my $i (1..$N) {
push @points, [rand(20)-10.0, rand(20)-10.0];
}
 
 
my ($a, $b, $d) = closest_pair_simple(\@points);
print "$d\n";
 
my ($a1, $b1, $d1) = closest_pair(\@points);
#print "$d1\n";

Time for the brute-force algorithm gave 40.63user 0.12system 0:41.06elapsed, while the divide&conqueer algorithm gave 0.37user 0.00system 0:00.38elapsed with 5000 points.

[edit] Perl 6

Translation of: Perl 5

We avoid taking square roots in the slow method because the squares are just as comparable. (This doesn't always work in the fast method because of distance assumptions in the algorithm.)

sub MAIN ($N = 5000) {
my @points = (^$N).map: { [rand * 20 - 10, rand * 20 - 10] }
 
my ($af, $bf, $df) = closest_pair(@points);
say "fast $df at [$af], [$bf]";
 
my ($as, $bs, $ds) = closest_pair_simple(@points);
say "slow $ds at [$as], [$bs]";
}
 
sub dist-squared($a,$b) {
($a[0] - $b[0]) ** 2 +
($a[1] - $b[1]) ** 2;
}
 
sub closest_pair_simple(@arr is copy) {
return Inf if @arr < 2;
my ($a, $b, $d) = @arr[0,1], dist-squared(|@arr[0,1]);
while @arr {
my $p = pop @arr;
for @arr -> $l {
my $t = dist-squared($p, $l);
($a, $b, $d) = $p, $l, $t if $t < $d;
}
}
return $a, $b, sqrt $d;
}
 
sub closest_pair(@r) {
my @ax = @r.sort: { .[0] }
my @ay = @r.sort: { .[1] }
return closest_pair_real(@ax, @ay);
}
 
sub closest_pair_real(@rx, @ry) {
return closest_pair_simple(@rx) if @rx <= 3;
 
my @xP = @rx;
my @yP = @ry;
my $N = @xP;
 
my $midx = ceiling($N/2)-1;
 
my @PL = @xP[0 .. $midx];
my @PR = @xP[$midx+1 ..^ $N];
 
my $xm = @xP[$midx][0];
 
my @yR;
my @yL;
push ($_[0] <= $xm ?? @yR !! @yL), $_ for @yP;
 
my ($al, $bl, $dL) = closest_pair_real(@PL, @yR);
my ($ar, $br, $dR) = closest_pair_real(@PR, @yL);
 
my ($m1, $m2, $dmin) = $dR < $dL
?? ($ar, $br, $dR)
!! ($al, $bl, $dL);
 
my @yS = @yP.grep: { abs($xm - .[0]) < $dmin }
 
if @yS {
my ($w1, $w2, $closest) = $m1, $m2, $dmin;
for 0 ..^ @yS.end -> $i {
for $i+1 ..^ @yS -> $k {
last unless @yS[$k][1] - @yS[$i][1] < $dmin;
my $d = sqrt dist-squared(@yS[$k], @yS[$i]);
($w1, $w2, $closest) = @yS[$k], @yS[$i], $d if $d < $closest;
}
 
}
return $w1, $w2, $closest;
 
} else {
return $m1, $m2, $dmin;
}
}

[edit] PicoLisp

(de closestPairBF (Lst)
(let Min T
(use (Pt1 Pt2)
(for P Lst
(for Q Lst
(or
(== P Q)
(>=
(setq N
(let (A (- (car P) (car Q)) B (- (cdr P) (cdr Q)))
(+ (* A A) (* B B)) ) )
Min )
(setq Min N Pt1 P Pt2 Q) ) ) )
(list Pt1 Pt2 (sqrt Min)) ) ) )

Test:

: (scl 6)
-> 6

: (closestPairBF
   (quote
      (0.654682 . 0.925557)
      (0.409382 . 0.619391)
      (0.891663 . 0.888594)
      (0.716629 . 0.996200)
      (0.477721 . 0.946355)
      (0.925092 . 0.818220)
      (0.624291 . 0.142924)
      (0.211332 . 0.221507)
      (0.293786 . 0.691701)
      (0.839186 . 0.728260) ) )
-> ((891663 . 888594) (925092 . 818220) 77910)

[edit] PL/I

 
/* Closest Pair Problem */
closest: procedure options (main);
declare n fixed binary;
 
get list (n);
begin;
declare 1 P(n),
2 x float,
2 y float;
declare (i, ii, j, jj) fixed binary;
declare (distance, min_distance initial (0) ) float;
 
get list (P);
min_distance = sqrt( (P.x(1) - P.x(2))**2 + (P.y(1) - P.y(2))**2 );
ii = 1; jj = 2;
do i = 1 to n;
do j = 1 to n;
distance = sqrt( (P.x(i) - P.x(j))**2 + (P.y(i) - P.y(j))**2 );
if distance > 0 then
if distance < min_distance then
do;
min_distance = distance;
ii = i; jj = j;
end;
end;
end;
put skip edit ('The minimum distance ', min_distance,
' is between the points [', P.x(ii),
',', P.y(ii), '] and [', P.x(jj), ',', P.y(jj), ']' )
(a, f(6,2));
end;
end closest;
 

[edit] PureBasic

Brute force version

Procedure.d bruteForceClosestPair(Array P.coordinate(1))
Protected N=ArraySize(P()), i, j
Protected mindistance.f=Infinity(), t.d
Shared a, b
If N<2
a=0: b=0
Else
For i=0 To N-1
For j=i+1 To N
t=Pow(Pow(P(i)\x-P(j)\x,2)+Pow(P(i)\y-P(j)\y,2),0.5)
If mindistance>t
mindistance=t
a=i: b=j
EndIf
Next
Next
EndIf
ProcedureReturn mindistance
EndProcedure
 

Implementation can be as

Structure coordinate
x.d
y.d
EndStructure
 
Dim DataSet.coordinate(9)
Define i, x.d, y.d, a, b
 
;- Load data from datasection
Restore DataPoints
For i=0 To 9
Read.d x: Read.d y
DataSet(i)\x=x
DataSet(i)\y=y
Next i
 
If OpenConsole()
PrintN("Mindistance= "+StrD(bruteForceClosestPair(DataSet()),6))
PrintN("Point 1= "+StrD(DataSet(a)\x,6)+": "+StrD(DataSet(a)\y,6))
PrintN("Point 2= "+StrD(DataSet(b)\x,6)+": "+StrD(DataSet(b)\y,6))
Print(#CRLF$+"Press ENTER to quit"): Input()
EndIf
 
DataSection
DataPoints:
Data.d 0.654682, 0.925557, 0.409382, 0.619391, 0.891663, 0.888594
Data.d 0.716629, 0.996200, 0.477721, 0.946355, 0.925092, 0.818220
Data.d 0.624291, 0.142924, 0.211332, 0.221507, 0.293786, 0.691701, 0.839186, 0.72826
EndDataSection
Output:
Mindistance= 0.077910
Point 1= 0.891663: 0.888594
Point 2= 0.925092: 0.818220

Press ENTER to quit

[edit] Python

"""
Compute nearest pair of points using two algorithms
 
First algorithm is 'brute force' comparison of every possible pair.
Second, 'divide and conquer', is based on:
www.cs.iupui.edu/~xkzou/teaching/CS580/Divide-and-conquer-closestPair.ppt
"""

 
from random import randint, randrange
from operator import itemgetter, attrgetter
 
infinity = float('inf')
 
# Note the use of complex numbers to represent 2D points making distance == abs(P1-P2)
 
def bruteForceClosestPair(point):
numPoints = len(point)
if numPoints < 2:
return infinity, (None, None)
return min( ((abs(point[i] - point[j]), (point[i], point[j]))
for i in range(numPoints-1)
for j in range(i+1,numPoints)),
key=itemgetter(0))
 
def closestPair(point):
xP = sorted(point, key= attrgetter('real'))
yP = sorted(point, key= attrgetter('imag'))
return _closestPair(xP, yP)
 
def _closestPair(xP, yP):
numPoints = len(xP)
if numPoints <= 3:
return bruteForceClosestPair(xP)
Pl = xP[:numPoints/2]
Pr = xP[numPoints/2:]
Yl, Yr = [], []
xDivider = Pl[-1].real
for p in yP:
if p.real <= xDivider:
Yl.append(p)
else:
Yr.append(p)
dl, pairl = _closestPair(Pl, Yl)
dr, pairr = _closestPair(Pr, Yr)
dm, pairm = (dl, pairl) if dl < dr else (dr, pairr)
# Points within dm of xDivider sorted by Y coord
closeY = [p for p in yP if abs(p.real - xDivider) < dm]
numCloseY = len(closeY)
if numCloseY > 1:
# There is a proof that you only need compare a max of 7 next points
closestY = min( ((abs(closeY[i] - closeY[j]), (closeY[i], closeY[j]))
for i in range(numCloseY-1)
for j in range(i+1,min(i+8, numCloseY))),
key=itemgetter(0))
return (dm, pairm) if dm <= closestY[0] else closestY
else:
return dm, pairm
 
def times():
''' Time the different functions
'''

import timeit
 
functions = [bruteForceClosestPair, closestPair]
for f in functions:
print 'Time for', f.__name__, timeit.Timer(
'%s(pointList)' % f.__name__,
'from closestpair import %s, pointList' % f.__name__).timeit(number=1)
 
 
 
pointList = [randint(0,1000)+1j*randint(0,1000) for i in range(2000)]
 
if __name__ == '__main__':
pointList = [(5+9j), (9+3j), (2+0j), (8+4j), (7+4j), (9+10j), (1+9j), (8+2j), 10j, (9+6j)]
print pointList
print ' bruteForceClosestPair:', bruteForceClosestPair(pointList)
print ' closestPair:', closestPair(pointList)
for i in range(10):
pointList = [randrange(11)+1j*randrange(11) for i in range(10)]
print '\n', pointList
print ' bruteForceClosestPair:', bruteForceClosestPair(pointList)
print ' closestPair:', closestPair(pointList)
print '\n'
times()
times()
times()
Output:
followed by timing comparisons

(Note how the two algorithms agree on the minimum distance, but may return a different pair of points if more than one pair of points share that minimum separation):

[(5+9j), (9+3j), (2+0j), (8+4j), (7+4j), (9+10j), (1+9j), (8+2j), 10j, (9+6j)]
  bruteForceClosestPair: (1.0, ((8+4j), (7+4j)))
            closestPair: (1.0, ((8+4j), (7+4j)))

[(10+6j), (7+0j), (9+4j), (4+8j), (7+5j), (6+4j), (1+9j), (6+4j), (1+3j), (5+0j)]
 bruteForceClosestPair: (0.0, ((6+4j), (6+4j)))
           closestPair: (0.0, ((6+4j), (6+4j)))

[(4+10j), (8+5j), (10+3j), (9+7j), (2+5j), (6+7j), (6+2j), (9+6j), (3+8j), (5+1j)]
 bruteForceClosestPair: (1.0, ((9+7j), (9+6j)))
           closestPair: (1.0, ((9+7j), (9+6j)))

[(10+0j), (3+10j), (10+7j), (1+8j), (5+10j), (8+8j), (4+7j), (6+2j), (6+10j), (9+3j)]
 bruteForceClosestPair: (1.0, ((5+10j), (6+10j)))
           closestPair: (1.0, ((5+10j), (6+10j)))

[(3+7j), (5+3j), 0j, (2+9j), (2+5j), (9+6j), (5+9j), (4+3j), (3+8j), (8+7j)]
 bruteForceClosestPair: (1.0, ((3+7j), (3+8j)))
           closestPair: (1.0, ((4+3j), (5+3j)))

[(4+3j), (10+9j), (2+7j), (7+8j), 0j, (3+10j), (10+2j), (7+10j), (7+3j), (1+4j)]
 bruteForceClosestPair: (2.0, ((7+8j), (7+10j)))
           closestPair: (2.0, ((7+8j), (7+10j)))

[(9+2j), (9+8j), (6+4j), (7+0j), (10+2j), (10+0j), (2+7j), (10+7j), (9+2j), (1+5j)]
 bruteForceClosestPair: (0.0, ((9+2j), (9+2j)))
           closestPair: (0.0, ((9+2j), (9+2j)))

[(3+3j), (8+2j), (4+0j), (1+1j), (9+10j), (5+0j), (2+3j), 5j, (5+0j), (7+0j)]
 bruteForceClosestPair: (0.0, ((5+0j), (5+0j)))
           closestPair: (0.0, ((5+0j), (5+0j)))

[(1+5j), (8+3j), (8+10j), (6+8j), (10+9j), (2+0j), (2+7j), (8+7j), (8+4j), (1+2j)]
 bruteForceClosestPair: (1.0, ((8+3j), (8+4j)))
           closestPair: (1.0, ((8+3j), (8+4j)))

[(8+4j), (8+6j), (8+0j), 0j, (10+7j), (10+6j), 6j, (1+3j), (1+8j), (6+9j)]
 bruteForceClosestPair: (1.0, ((10+7j), (10+6j)))
           closestPair: (1.0, ((10+7j), (10+6j)))

[(6+8j), (10+1j), 3j, (7+9j), (4+10j), (4+7j), (5+7j), (6+10j), (4+7j), (2+4j)]
 bruteForceClosestPair: (0.0, ((4+7j), (4+7j)))
           closestPair: (0.0, ((4+7j), (4+7j)))


Time for bruteForceClosestPair 4.57953371169
Time for closestPair 0.122539596513
Time for bruteForceClosestPair 5.13221177552
Time for closestPair 0.124602707886
Time for bruteForceClosestPair 4.83609397284
Time for closestPair 0.119326618327
>>> 

[edit] R

Works with: R version 2.8.1+

Brute force solution as per wikipedia pseudo-code

closest_pair_brute <-function(x,y,plotxy=F) { 
xy = cbind(x,y)
cp = bruteforce(xy)
cat("\n\nShortest path found = \n From:\t\t(",cp[1],',',cp[2],")\n To:\t\t(",cp[3],',',cp[4],")\n Distance:\t",cp[5],"\n\n",sep="")
if(plotxy) {
plot(x,y,pch=19,col='black',main="Closest Pair", asp=1)
points(cp[1],cp[2],pch=19,col='red')
points(cp[3],cp[4],pch=19,col='red')
}
distance <- function(p1,p2) {
x1 = (p1[1])
y1 = (p1[2])
x2 = (p2[1])
y2 = (p2[2])
sqrt((x2-x1)^2 + (y2-y1)^2)
}
bf_iter <- function(m,p,idx=NA,d=NA,n=1) {
dd = distance(p,m[n,])
if((is.na(d) || dd<=d) && p!=m[n,]){d = dd; idx=n;}
if(n == length(m[,1])) { c(m[idx,],d) }
else bf_iter(m,p,idx,d,n+1)
}
bruteforce <- function(pmatrix,n=1,pd=c(NA,NA,NA,NA,NA)) {
p = pmatrix[n,]
ppd = c(p,bf_iter(pmatrix,p))
if(ppd[5]<pd[5] || is.na(pd[5])) pd = ppd
if(n==length(pmatrix[,1])) pd
else bruteforce(pmatrix,n+1,pd)
}
}

Quicker brute force solution for R that makes use of the apply function native to R for dealing with matrices. It expects x and y to take the form of separate vectors.

closestPair<-function(x,y)
{
distancev <- function(pointsv)
{
x1 <- pointsv[1]
y1 <- pointsv[2]
x2 <- pointsv[3]
y2 <- pointsv[4]
sqrt((x1 - x2)^2 + (y1 - y2)^2)
}
pairstocompare <- t(combn(length(x),2))
pointsv <- cbind(x[pairstocompare[,1]],y[pairstocompare[,1]],x[pairstocompare[,2]],y[pairstocompare[,2]])
pairstocompare <- cbind(pairstocompare,apply(pointsv,1,distancev))
minrow <- pairstocompare[pairstocompare[,3] == min(pairstocompare[,3])]
if (!is.null(nrow(minrow))) {print("More than one point at this distance!"); minrow <- minrow[1,]}
cat("The closest pair is:\n\tPoint 1: ",x[minrow[1]],", ",y[minrow[1]],
"\n\tPoint 2: ",x[minrow[2]],", ",y[minrow[2]],
"\n\tDistance: ",minrow[3],"\n",sep="")
c(distance=minrow[3],x1.x=x[minrow[1]],y1.y=y[minrow[1]],x2.x=x[minrow[2]],y2.y=y[minrow[2]])
}


This is the quickest version, that makes use of the 'dist' function of R. It takes a two-column object of x,y-values as input, or creates such an object from seperate x and y-vectors.

closest.pairs <- function(x, y=NULL, ...){
# takes two-column object(x,y-values), or creates such an object from x and y values
if(!is.null(y)) x <- cbind(x, y)
 
distances <- dist(x)
min.dist <- min(distances)
point.pair <- combn(1:nrow(x), 2)[, which.min(distances)]
 
cat("The closest pair is:\n\t",
sprintf("Point 1: %.3f, %.3f \n\tPoint 2: %.3f, %.3f \n\tDistance: %.3f.\n",
x[point.pair[1],1], x[point.pair[1],2],
x[point.pair[2],1], x[point.pair[2],2],
min.dist),
sep="" )
c( x1=x[point.pair[1],1],y1=x[point.pair[1],2],
x2=x[point.pair[2],1],y2=x[point.pair[2],2],
distance=min.dist)
}
Example
x = (sample(-1000.00:1000.00,100))
y = (sample(-1000.00:1000.00,length(x)))
cp = closest.pairs(x,y)
#cp = closestPair(x,y)
plot(x,y,pch=19,col='black',main="Closest Pair", asp=1)
points(cp["x1.x"],cp["y1.y"],pch=19,col='red')
points(cp["x2.x"],cp["y2.y"],pch=19,col='red')
#closest_pair_brute(x,y,T)
 
Performance
system.time(closest_pair_brute(x,y), gcFirst = TRUE)
Shortest path found =
From: (32,-987)
To: (25,-993)
Distance: 9.219544
 
user system elapsed
0.35 0.02 0.37
 
system.time(closest.pairs(x,y), gcFirst = TRUE)
The closest pair is:
Point 1: 32.000, -987.000
Point 2: 25.000, -993.000
Distance: 9.220.
 
user system elapsed
0.08 0.00 0.10
 
system.time(closestPair(x,y), gcFirst = TRUE)
The closest pair is:
Point 1: 32, -987
Point 2: 25, -993
Distance: 9.219544
 
user system elapsed
0.17 0.00 0.19
 
 

Using dist function for brute force, but divide and conquer (as per pseudocode) for speed:

closest.pairs.bruteforce <- function(x, y=NULL)
{
if (!is.null(y))
{
x <- cbind(x,y)
}
d <- dist(x)
cp <- x[combn(1:nrow(x), 2)[, which.min(d)],]
list(p1=cp[1,], p2=cp[2,], d=min(d))
}
 
closest.pairs.dandc <- function(x, y=NULL)
{
if (!is.null(y))
{
x <- cbind(x,y)
}
if (sd(x[,"x"]) < sd(x[,"y"]))
{
x <- cbind(x=x[,"y"],y=x[,"x"])
swap <- TRUE
}
else
{
swap <- FALSE
}
xp <- x[order(x[,"x"]),]
.cpdandc.rec <- function(xp,yp)
{
n <- dim(xp)[1]
if (n <= 4)
{
closest.pairs.bruteforce(xp)
}
else
{
xl <- xp[1:floor(n/2),]
xr <- xp[(floor(n/2)+1):n,]
cpl <- .cpdandc.rec(xl)
cpr <- .cpdandc.rec(xr)
if (cpl$d<cpr$d) cp <- cpl else cp <- cpr
cp
}
}
cp <- .cpdandc.rec(xp)
 
yp <- x[order(x[,"y"]),]
xm <- xp[floor(dim(xp)[1]/2),"x"]
ys <- yp[which(abs(xm - yp[,"x"]) <= cp$d),]
nys <- dim(ys)[1]
if (!is.null(nys) && nys > 1)
{
for (i in 1:(nys-1))
{
k <- i + 1
while (k <= nys && ys[i,"y"] - ys[k,"y"] < cp$d)
{
d <- sqrt((ys[k,"x"]-ys[i,"x"])^2 + (ys[k,"y"]-ys[i,"y"])^2)
if (d < cp$d) cp <- list(p1=ys[i,],p2=ys[k,],d=d)
k <- k + 1
}
}
}
if (swap)
{
list(p1=cbind(x=cp$p1["y"],y=cp$p1["x"]),p2=cbind(x=cp$p2["y"],y=cp$p2["x"]),d=cp$d)
}
else
{
cp
}
}
 
# Test functions
cat("How many points?\n")
n <- scan(what=integer(),n=1)
x <- rnorm(n)
y <- rnorm(n)
tstart <- proc.time()[3]
cat("Closest pairs divide and conquer:\n")
print(cp <- closest.pairs.dandc(x,y))
cat(sprintf("That took %.2f seconds.\n",proc.time()[3] - tstart))
plot(x,y)
points(c(cp$p1["x"],cp$p2["x"]),c(cp$p1["y"],cp$p2["y"]),col="red")
tstart <- proc.time()[3]
cat("\nClosest pairs brute force:\n")
print(closest.pairs.bruteforce(x,y))
cat(sprintf("That took %.2f seconds.\n",proc.time()[3] - tstart))
 
Output:
How many points?
1: 500
Read 1 item
Closest pairs divide and conquer:
$p1
         x          y 
1.68807938 0.05876328 

$p2
         x          y 
1.68904694 0.05878173 

$d
[1] 0.0009677302

That took 0.43 seconds.

Closest pairs brute force:
$p1
         x          y 
1.68807938 0.05876328 

$p2
         x          y 
1.68904694 0.05878173 

$d
[1] 0.0009677302

That took 6.38 seconds.

[edit] Racket

The brute force solution using complex numbers to represent pairs.

 
#lang racket
(define (dist z0 z1) (magnitude (- z1 z0)))
(define (dist* zs) (apply dist zs))
 
(define (closest-pair zs)
(if (< (length zs) 2)
-inf.0
(first
(sort (for/list ([z0 zs])
(list z0 (argmin (λ(z) (if (= z z0) +inf.0 (dist z z0))) zs)))
< #:key dist*))))
 
(define result (closest-pair '(0+1i 1+2i 3+4i)))
(displayln (~a "Closest points: " result))
(displayln (~a "Distance: " (dist* result)))
 
Output:
 
Closest points: (0+1i 1+2i)
Distance: 1.4142135623730951
 

[edit] REXX

(This version of the REXX program is modeled after the psuedo-code.)

/*REXX program to solve the  closest pair problem.                      */
parse arg N low high seed .; if n=='' then n=100
if seed\=='' then call random ,,seed /*use seed, makes rand repeatable*/
if low=='' | low==',' then low=0
if high=='' | high==',' then high=20000
w=length(high); w=w + (w//2==0)
do j=1 for N /*gen random points*/
@.j.xx=random(low,high)
@.j.yy=random(low,high)
end /*j*/
nearA=1
nearB=2
minDD=(@.nearA.xx-@.nearB.xx)**2 + (@.nearA.yy-@.nearB.yy)**2
 
do j=1 for N-1
do k=j+1 to N
dd=(@.j.xx-@.k.xx)**2 + (@.j.yy-@.k.yy)**2
if dd\=0 then if dd<minDD then do; minDD=dd
nearA=j
nearB=k
end
end /*k*/
end /*j*/
 
say 'For' N "points:"; say
say ' 'center('x',w,"═")' ' center('y',w,"═")
say 'The points ['right(@.nearA.xx,w)"," right(@.nearA.yy,w)"]" ' and'
say ' ['right(@.nearB.xx,w)"," right(@.nearB.yy,w)"]"; say
say 'the minimum distance between them is: ' sqrt(abs(minDD))
exit /*stick a fork in it, we're done.*/
/*───────────────────────────────────sqrt subroutine────────────────────*/
sqrt: procedure; parse arg x; if x=0 then return 0;d=digits();numeric digits 11
g=.sqrtG(); do j=0 while p>9; m.j=p; p=p%2+1; end
do k=j+5 to 0 by -1; if m.k>11 then numeric digits m.k; g=.5*(g+x/g); end
numeric digits d; return g/1
.sqrtG: numeric form; m.=11; p=d+d%4+2
parse value format(x,2,1,,0) 'E0' with g 'E' _ .; return g*.5'E'_%2
Output:
when using the input of: 200
For 100 points:

                 ══x══  ══y══
The points      [11341, 19534]   and
                [11136, 19498]

the minimum distance between them is:  208.136974

[edit] Ruby

Point = Struct.new(:x, :y)
 
def distance(p1, p2)
Math.hypot(p1.x - p2.x, p1.y - p2.y)
end
 
def closest_bruteforce(points)
mindist, minpts = Float::MAX, []
points.length.times do |i|
(i+1).upto(points.length - 1) do |j|
dist = distance(points[i], points[j])
if dist < mindist
mindist = dist
minpts = [points[i], points[j]]
end
end
end
[mindist, minpts]
end
 
def closest_recursive(points)
if points.length <= 3
return closest_bruteforce(points)
end
xP = points.sort_by {|p| p.x}
mid = (points.length / 2.0).ceil
pL = xP[0,mid]
pR = xP[mid..-1]
dL, pairL = closest_recursive(pL)
dR, pairR = closest_recursive(pR)
if dL < dR
dmin, dpair = dL, pairL
else
dmin, dpair = dR, pairR
end
yP = xP.find_all {|p| (pL[-1].x - p.x).abs < dmin}.sort_by {|p| p.y}
closest = Float::MAX
closestPair = []
0.upto(yP.length - 2) do |i|
(i+1).upto(yP.length - 1) do |k|
break if (yP[k].y - yP[i].y) >= dmin
dist = distance(yP[i], yP[k])
if dist < closest
closest = dist
closestPair = [yP[i], yP[k]]
end
end
end
if closest < dmin
[closest, closestPair]
else
[dmin, dpair]
end
end
 
 
points = Array.new(100) {Point.new(rand, rand)}
p ans1 = closest_bruteforce(points)
p ans2 = closest_recursive(points)
fail "bogus!" if ans1[0] != ans2[0]
 
require 'benchmark'
 
points = Array.new(10000) {Point.new(rand, rand)}
Benchmark.bm(12) do |x|
x.report("bruteforce") {ans1 = closest_bruteforce(points)}
x.report("recursive") {ans2 = closest_recursive(points)}
end
[0.00522229060545241, [#<struct Point x=0.43887011964135, y=0.00656904813877568>, #<struct Point x=0.433711197400243, y=0.00575797448120408>]]
[0.00522229060545241, [#<struct Point x=0.433711197400243, y=0.00575797448120408>, #<struct Point x=0.43887011964135, y=0.00656904813877568>]]
                  user     system      total        real
bruteforce  133.437000   0.000000 133.437000 (134.633000)
recursive     0.516000   0.000000   0.516000 (  0.559000)

[edit] Run BASIC

Courtesy http://dkokenge.com/rbp

n =10                              ' 10 data points input
dim x(n)
dim y(n)
 
pt1 = 0 ' 1st point
pt2 = 0 ' 2nd point
 
for i =1 to n ' read in data
read x(i)
read y(i)
next i
 
minDist = 1000000
 
for i =1 to n -1
for j =i +1 to n
distXsq =(x(i) -x(j))^2
disYsq =(y(i) -y(j))^2
d =abs((dxSq +disYsq)^0.5)
if d <minDist then
minDist =d
pt1 =i
pt2 =j
end if
next j
next i
 
print "Distance ="; minDist; " between ("; x(pt1); ", "; y(pt1); ") and ("; x(pt2); ", "; y(pt2); ")"
 
end
 
data 0.654682, 0.925557
data 0.409382, 0.619391
data 0.891663, 0.888594
data 0.716629, 0.996200
data 0.477721, 0.946355
data 0.925092, 0.818220
data 0.624291, 0.142924
data 0.211332, 0.221507
data 0.293786, 0.691701
data 0.839186, 0.72826

[edit] Scala

import scala.util._
 
object ClosestPair{
class Point(x:Double, y:Double) extends Pair(x,y){
def distance(p:Point)=math.hypot(_1-p._1, _2-p._2)
}
 
def closestPairBF(a:Array[Point])={
var minDist=a(0) distance a(1)
var minPoints=(a(0), a(1))
 
for(p1<-a; p2<-a if p1.ne(p2); dist=p1 distance p2 if(dist<minDist)){
minDist=dist;
minPoints=(p1, p2)
}
(minPoints, minDist)
}
 
def main(args: Array[String]): Unit = {
val a=Array.fill(1000)(new Point(Random.nextDouble, Random.nextDouble))
val (points, dist)=closestPairBF(a)
println("min: "+points._1+" - "+points._2+" ->"+dist)
}
}

[edit] Seed7

This is the brute force algorithm:

const type: point is new struct
var float: x is 0.0;
var float: y is 0.0;
end struct;
 
const func float: distance (in point: p1, in point: p2) is
return sqrt((p1.x-p2.x)**2+(p1.y-p2.y)**2);
 
const func array point: closest_pair (in array point: points) is func
result
var array point: result is 0 times point.value;
local
var float: dist is 0.0;
var float: minDistance is Infinity;
var integer: i is 0;
var integer: j is 0;
var integer: savei is 0;
var integer: savej is 0;
begin
for i range 1 to pred(length(points)) do
for j range succ(i) to length(points) do
dist := distance(points[i], points[j]);
if dist < minDistance then
minDistance := dist;
savei := i;
savej := j;
end if;
end for;
end for;
if minDistance <> Infinity then
result := [] (points[savei], points[savej]);
end if;
end func;

[edit] Smalltalk

See Closest-pair problem/Smalltalk

[edit] Tcl

Each point is represented as a list of two floating-point numbers, the first being the x coordinate, and the second being the y.

package require Tcl 8.5
 
# retrieve the x-coordinate
proc x p {lindex $p 0}
# retrieve the y-coordinate
proc y p {lindex $p 1}
 
proc distance {p1 p2} {
expr {hypot(([x $p1]-[x $p2]), ([y $p1]-[y $p2]))}
}
 
proc closest_bruteforce {points} {
set n [llength $points]
set mindist Inf
set minpts {}
for {set i 0} {$i < $n - 1} {incr i} {
for {set j [expr {$i + 1}]} {$j < $n} {incr j} {
set p1 [lindex $points $i]
set p2 [lindex $points $j]
set dist [distance $p1 $p2]
if {$dist < $mindist} {
set mindist $dist
set minpts [list $p1 $p2]
}
}
}
return [list $mindist $minpts]
}
 
proc closest_recursive {points} {
set n [llength $points]
if {$n <= 3} {
return [closest_bruteforce $points]
}
set xP [lsort -real -increasing -index 0 $points]
set mid [expr {int(ceil($n/2.0))}]
set PL [lrange $xP 0 [expr {$mid-1}]]
set PR [lrange $xP $mid end]
set procname [lindex [info level 0] 0]
lassign [$procname $PL] dL pairL
lassign [$procname $PR] dR pairR
if {$dL < $dR} {
set dmin $dL
set dpair $pairL
} else {
set dmin $dR
set dpair $pairR
}
 
set xM [x [lindex $PL end]]
foreach p $xP {
if {abs($xM - [x $p]) < $dmin} {
lappend S $p
}
}
set yP [lsort -real -increasing -index 1 $S]
set closest Inf
set nP [llength $yP]
for {set i 0} {$i <= $nP-2} {incr i} {
set yPi [lindex $yP $i]
for {set k [expr {$i+1}]; set yPk [lindex $yP $k]} {
$k < $nP-1 && ([y $yPk]-[y $yPi]) < $dmin
} {incr k; set yPk [lindex $yP $k]} {
set dist [distance $yPk $yPi]
if {$dist < $closest} {
set closest $dist
set closestPair [list $yPi $yPk]
}
}
}
expr {$closest < $dmin ? [list $closest $closestPair] : [list $dmin $dpair]}
}
 
# testing
set N 10000
for {set i 1} {$i <= $N} {incr i} {
lappend points [list [expr {rand()*100}] [expr {rand()*100}]]
}
 
# instrument the number of calls to [distance] to examine the
# efficiency of the recursive solution
trace add execution distance enter comparisons
proc comparisons args {incr ::comparisons}
 
puts [format "%-10s  %9s  %9s  %s" method compares time closest]
foreach method {bruteforce recursive} {
set ::comparisons 0
set time [time {set ::dist($method) [closest_$method $points]} 1]
puts [format "%-10s  %9d  %9d  %s" $method $::comparisons [lindex $time 0] [lindex $::dist($method) 0]]
}
Output:
method      compares      time closest
bruteforce  49995000 512967207 0.0015652738546658382
recursive      14613    488094 0.0015652738546658382

Note that the lindex and llength commands are both O(1).

[edit] Ursala

The brute force algorithm is easy. Reading from left to right, clop is defined as a function that forms the Cartesian product of its argument, and then extracts the member whose left side is a minimum with respect to the floating point comparison relation after deleting equal pairs and attaching to the left of each remaining pair the sum of the squares of the differences between corresponding coordinates.

#import flo
 
clop = @iiK0 fleq$-&l+ *EZF ^\~& plus+ sqr~~+ minus~~bbI

The divide and conquer algorithm following the specification given above is a little more hairy but not much longer. The eudist library function is used to compute the distance between points.

#import std
#import flo
 
clop =
 
^(fleq-<&l,fleq-<&r); @blrNCCS ~&lrbhthPX2X+ ~&a^& fleq$-&l+ leql/8?al\^(eudist,~&)*altK33htDSL -+
^C/~&rr ^(eudist,~&)*tK33htDSL+ @rlrlPXPlX ~| fleq^\~&lr abs+ minus@llPrhPX,
^/~&ar @farlK30K31XPGbrlrjX3J ^/~&arlhh @W lesser fleq@bl+-

test program:

test_data =
 
<
(1.547290e+00,3.313053e+00),
(5.250805e-01,-7.300260e+00),
(7.062114e-02,1.220251e-02),
(-4.473024e+00,-5.393712e+00),
(-2.563714e+00,-3.595341e+00),
(-2.132372e+00,2.358850e+00),
(2.366238e+00,-9.678425e+00),
(-1.745694e+00,3.276434e+00),
(8.066843e+00,-9.101268e+00),
(-8.256901e+00,-8.717900e+00),
(7.397744e+00,-5.366434e+00),
(2.060291e-01,2.840891e+00),
(-6.935319e+00,-5.192438e+00),
(9.690418e+00,-9.175753e+00),
(3.448993e+00,2.119052e+00),
(-7.769218e+00,4.647406e-01)>
 
#cast %eeWWA
 
example = clop test_data
Output:

The output shows the minimum distance and the two points separated by that distance. (If the brute force algorithm were used, it would have displayed the square of the distance.)

9.957310e-01: (
   (-2.132372e+00,2.358850e+00),
   (-1.745694e+00,3.276434e+00))

[edit] XPL0

The brute force method is simpler than the recursive solution and is perfectly adequate, even for a thousand points.

include c:\cxpl\codes;          \intrinsic 'code' declarations
 
proc ClosestPair(P, N); \Show closest pair of points in array P
real P; int N;
real Dist2, MinDist2;
int I, J, SI, SJ;
[MinDist2:= 1e300;
for I:= 0 to N-2 do
[for J:= I+1 to N-1 do
[Dist2:= sq(P(I,0)-P(J,0)) + sq(P(I,1)-P(J,1));
if Dist2 < MinDist2 then \squared distances are sufficient for compares
[MinDist2:= Dist2;
SI:= I; SJ:= J;
];
];
];
IntOut(0, SI); Text(0, " -- "); IntOut(0, SJ); CrLf(0);
RlOut(0, P(SI,0)); Text(0, ","); RlOut(0, P(SI,1));
Text(0, " -- ");
RlOut(0, P(SJ,0)); Text(0, ","); RlOut(0, P(SJ,1));
CrLf(0);
];
 
real Data;
[Format(1, 6);
Data:= [[0.654682, 0.925557], \0 test data from BASIC examples
[0.409382, 0.619391], \1
[0.891663, 0.888594], \2
[0.716629, 0.996200], \3
[0.477721, 0.946355], \4
[0.925092, 0.818220], \5
[0.624291, 0.142924], \6
[0.211332, 0.221507], \7
[0.293786, 0.691701], \8
[0.839186, 0.728260]]; \9
ClosestPair(Data, 10);
]
Output:
2 -- 5
0.891663,0.888594 -- 0.925092,0.818220

[edit] zkl

An ugly solution in both time and space.

class Point{
fcn init(_x,_y){ var[const] x=_x, y=_y; }
fcn distance(p){(p.x-x).hypot(p.y-y)}
fcn toString{String("Point(",x,",",y,")")}
}
 
// find closest two points using brute ugly force:
// find all combinations of two points, measure distance, pick smallest
fcn closestPoints(points){
pairs:=Utils.Helpers.pickNFrom(2,points);
triples:=pairs.apply(fcn([(p1,p2)]){T(p1,p2,p1.distance(p2))});
triples.reduce(fcn([(_,_,d1)]p1,[(_,_,d2)]p2){
if(d1 < d2) p1 else p2});
}
points:=T( 5.0, 9.0,  9.0, 3.0,
2.0, 0.0, 8.0, 4.0,
7.0, 4.0, 9.0, 10.0,
1.0, 9.0, 8.0, 2.0,
0.0, 10.0, 9.0, 6.0 ).pump(List,T.fp(Void.Read,1),Point);
 
closestPoints(points).println(); //-->L(Point(8,4),Point(7,4),1)
 
points:=T( 0.654682, 0.925557, 0.409382, 0.619391,
0.891663, 0.888594, 0.716629, 0.9962,
0.477721, 0.946355, 0.925092, 0.81822,
0.624291, 0.142924, 0.211332, 0.221507,
0.293786, 0.691701, 0.839186, 0.72826)
.pump(List,T.fp(Void.Read,1),Point);
closestPoints(points).println();
//-->L(Point(0.891663,0.888594),Point(0.925092,0.81822),0.0779102)
Personal tools
Namespaces

Variants
Actions
Community
Explore
Misc
Toolbox