Closest-pair problem: Difference between revisions
Line 264: | Line 264: | ||
And divide-and-conquer. Notice that the code has been written for brevity and to demonstrate the algorithm, not true optimization. There are further language-specific optimizations that could be applied. |
And divide-and-conquer. Notice that the code has been written for brevity and to demonstrate the algorithm, not true optimization. There are further language-specific optimizations that could be applied. |
||
<lang csharp> |
<lang csharp> |
||
public static Segment MyClosestDivide(List<PointF> points) |
|||
{ |
{ |
||
return MyClosestRec(points.OrderBy(p => p.X).ToList()); |
|||
} |
|||
public static Segment MyClosestRec(List<PointF> pointsByX) |
|||
int split = points.Count() / 2; |
|||
{ |
|||
var ordered = points.OrderBy(point => point.X); |
|||
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 leftResult = MyClosestRec(leftByX); |
|||
⚫ | |||
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 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 leftMin = Closest_Recursive(pointsOnLeft); |
|||
{ |
|||
float leftDist = leftMin.Length(); |
|||
var pUpper = inBandByY[j]; |
|||
var rightMin = Closest_Recursive(pointsOnRight); |
|||
float rightDist = rightMin.Length(); |
|||
// Comparing each point to successivly increasing Y values |
|||
float minDist = Math.Min(leftDist, rightDist); |
|||
// Thus, can terminate as soon as deltaY is greater than best result |
|||
⚫ | |||
if ((pUpper.Y - pLower.Y) >= result.Length()) |
|||
var closeY = pointsOnRight.TakeWhile(point => point.X - xDivider < minDist).OrderBy(point => point.Y); |
|||
break; |
|||
if (Segment.Length(pLower, pUpper) < result.Length()) |
|||
var crossingPairs = pointsOnLeft.SkipWhile(point => xDivider - point.X > minDist) |
|||
result = new Segment(pLower, pUpper); |
|||
.SelectMany(p1 => closeY.SkipWhile(i => i.Y < p1.Y - minDist) |
|||
} |
|||
.TakeWhile(i => i.Y < p1.Y + minDist) |
|||
} |
|||
.Select(p2 => new Segment( p1, p2 ))); |
|||
return result; |
|||
return crossingPairs.Union( new [] { leftMin, rightMin }) |
|||
} |
|||
.OrderBy(segment => segment.Length()).First(); |
|||
</lang> |
|||
However, the difference in speed is still remarkable. |
However, the difference in speed is still remarkable. |
||
{{Lines_too_long}} |
|||
<lang csharp>var randomizer = new Random(10); |
<lang csharp>var randomizer = new Random(10); |
||
var points = Enumerable.Range( 0, 10000).Select( i => new PointF( (float)randomizer.NextDouble(), (float)randomizer.NextDouble())).ToList(); |
var points = Enumerable.Range( 0, 10000).Select( i => new PointF( (float)randomizer.NextDouble(), (float)randomizer.NextDouble())).ToList(); |
Revision as of 23:37, 29 May 2011
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 return ∞ else 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
- Closest pair of points problem
- Closest Pair (McGill)
- Closest Pair (UCBS)
- Closest pair (WUStL)
- Closest pair (IUPUI)
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) <lang Ada>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;</lang>
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
C
Clojure
<lang 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])))))
</lang>
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.
<lang lisp>(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))))</lang>
C#
We provide a small helper class for distance comparisons: <lang csharp>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); }
}</lang>
Brute force: <lang csharp>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;
}</lang>
And divide-and-conquer. Notice that the code has been written for brevity and to demonstrate the algorithm, not true optimization. There are further language-specific optimizations that could be applied.
<lang csharp>
public static Segment MyClosestDivide(List<PointF> points)
{
return MyClosestRec(points.OrderBy(p => p.X).ToList());
}
public 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; } </lang>
However, the difference in speed is still remarkable. <lang csharp>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());</lang>
Output: <lang>Time used (Brute force) (float): 145731.8935 ms Time used (Divide & Conquer): 1139.2111 ms</lang>
Non Linq Brute Force: <lang csharp>
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; }</lang>
Targeted Search: Similar concept to recursive divide-and-conquer above, but I think much simpler. Also, much faster in my tests. 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.
<lang csharp>
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; }
</lang>
D
D version 2 implementation, designed for compactness: <lang d>import std.stdio,std.typecons,std.math,std.algorithm,std.array,std.random;
auto bruteForceClosestPair(T)(T[] points) {
auto minD = typeof(T.re).infinity; typeof(points[0]) minI, minJ; foreach (i; 0 .. points.length-1) foreach (j; i+1 .. points.length) { double dist = abs(points[i] - points[j]); if (dist < minD) { minD = dist; minI = points[i]; minJ = points[j]; } } return tuple(minD, minI, minJ);
}
auto closestPair(T)(T[] points) {
static Tuple!(typeof(T.re),T,T) inner(T[] xP, T[] yP) { if (xP.length <= 3) return bruteForceClosestPair(xP); auto Pl = xP[0 .. xP.length/2]; auto Pr = xP[xP.length/2 .. $]; auto xDiv = Pl[$ - 1].re; auto Yr = partition!((T p){ return p.re <= xDiv; })(yP); auto Yl = yP[0 .. yP.length - Yr.length]; auto dl_pairl = inner(Pl, Yl); auto dr_pairr = inner(Pr, Yr); auto dm_pairm = (dl_pairl[0] < dr_pairr[0]) ? dl_pairl : dr_pairr; auto dm = dm_pairm[0]; auto nextY = array(filter!((T p){return abs(p.re - xDiv) < dm;})(yP)); if (nextY.length > 1) { auto minD = typeof(T.re).infinity; size_t minI, minJ; foreach (i; 0 .. nextY.length-1) foreach (j; i+1 .. min(i+8, nextY.length)) { double dist = abs(nextY[i] - nextY[j]); if (dist < minD) { minD = dist; minI = i; minJ = j; } } return dm <= minD ? dm_pairm : tuple(minD,nextY[minI],nextY[minJ]); } else return dm_pairm; }
sort!("a.re < b.re")(points); auto xP = points.dup; sort!("a.im < b.im")(points); return inner(xP, points);
}
void main() {
auto pts = [5+9i, 9+3i, 2, 8+4i, 7+4i, 9+10i, 1+9i, 8+2i, 10i, 9+6i]; writeln(pts); writeln("bruteForceClosestPair: ", bruteForceClosestPair(pts)); writeln(" closestPair: ", closestPair(pts));
auto rnd = Random(1); // set seed auto points = new cdouble[10_000]; foreach (ref p; points) p = uniform(0.0, 1000.0, rnd) + uniform(0.0, 1000.0, rnd) * 1i; writeln("bruteForceClosestPair: ", bruteForceClosestPair(points)); writeln(" closestPair: ", closestPair(points));
}</lang> Output D2 version:
[5+9i, 9+3i, 2+0i, 8+4i, 7+4i, 9+10i, 1+9i, 8+2i, 0+10i, 9+6i] bruteForceClosestPair: Tuple!(double,cdouble,cdouble)(1, 8+4i, 7+4i) closestPair: Tuple!(double,cdouble,cdouble)(1, 8+4i, 7+4i) bruteForceClosestPair: Tuple!(double,cdouble,cdouble)(0.0947756, 643.463+788.31i, 643.407+788.387i) closestPair: Tuple!(double,cdouble,cdouble)(0.0947756, 643.407+788.387i, 643.463+788.31i)
Time gave 0:04.03 elapsed for brute force version, and 0:00.06 elapsed for the divide & conquer one (10_000 points).
A more efficient implementation of the brute-force algorithm, D version 1 with Tango library, LDC compiler: <lang d>import tango.stdc.stdio, tango.stdc.stdlib, tango.math.Math;
int bfClosestPair2(cdouble[] points, out size_t i1, out size_t i2) {
auto minD = typeof(points[0].re).infinity; if (points.length < 2) { i1 = i2 = size_t.max; return -1; } size_t minI, minJ; for (int i = 0; i < points.length-1; i++) { auto points_i_re = points[i].re; auto points_i_im = points[i].im; for (int j = i+1; j < points.length; j++) { auto dre = points_i_re - points[j].re; auto dist = dre * dre; if (dist < minD) { auto dim = points_i_im - points[j].im; dist += dim * dim; if (dist < minD) { minD = dist; minI = i; minJ = j; } } } } i1 = minI; i2 = minJ; return 0;
}
void main() {
srand(31415); auto pts = new cdouble[10_000]; foreach (ref p; pts) p = 1000.0 * (cast(double)rand() / (RAND_MAX + 1.0)) + 1000.0i * (cast(double)rand() / (RAND_MAX + 1.0)); size_t i, j; int err = bfClosestPair2(pts, i, j); if (err < 0) return; double d = sqrt((pts[i].re - pts[j].re) * (pts[i].re - pts[j].re) + (pts[i].im - pts[j].im) * (pts[i].im - pts[j].im)); printf("Closest pair: dist: %lf p1, p2: (%lf, %lf), (%lf, %lf)\n", d, pts[i].re, pts[i].im, pts[j].re, pts[j].im);
}</lang> Time gave 0:00.15 elapsed for brute-force version 2 (10_000 points).
F#
Brute force: <lang fsharp> 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)
</lang> For example: <lang fsharp> closest_pairs
[|Point(0.0, 0.0); Point(1.0, 0.0); Point (2.0, 2.0)|]
</lang> gives: <lang fsharp> (0,0, 1,0) </lang>
Fantom
(Based on the Ruby example.)
<lang fantom> 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") }
} </lang>
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
Fortran
See Closest pair problem/Fortran
Haskell
BF solution: <lang Haskell>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</lang>
Output example: <lang Haskell>*Main> testCP ([[0.8347201880148426,0.40774840545089647],[0.8348731214261784,0.4087113189531284]],9.749825850154334e-4) (4.02 secs, 488869056 bytes)</lang>
J
Solution of the simpler (brute-force) problem: <lang j>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</lang> Examples of use: <lang j> ]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| | +-----------------+---------+</lang> The program also works for higher dimensional vectors: <lang j> ]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| | +------------------------------------+--------+</lang>
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.
<lang javascript>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 }; }
}</lang>
MATLAB
This solution is an almost direct translation of the above pseudo-code into MATLAB. <lang 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</lang>
Sample Output: <lang MATLAB>[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</lang>
Objective-C
See Closest-pair problem/Objective-C
Oz
Translation of pseudocode: <lang oz>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}}</lang>
Perl
<lang 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";</lang>
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.
Perl 6
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.) <lang perl6>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; }
}</lang>
PicoLisp
<lang 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)) ) ) )</lang>
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)
PL/I
<lang> /* 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; </lang>
PureBasic
Brute force version <lang PureBasic>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 </lang>
Implementation can be as <lang PureBasic>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</lang> Output
Mindistance= 0.077910 Point 1= 0.891663: 0.888594 Point 2= 0.925092: 0.818220 Press ENTER to quit
Python
<lang 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()</lang>
Sample 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 >>>
R
Brute force solution as per wikipedia pseudo-code <lang R>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) }
}</lang>
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. <lang R>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]]) }</lang>
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.
<lang R>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) }</lang>
Example<lang R>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
</lang>
Ruby
<lang 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</lang>
[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)
Scala
<lang 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) }
}</lang>
Seed7
This is the brute force algorithm:
<lang seed7>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;</lang>
Smalltalk
See Closest-pair problem/Smalltalk
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. <lang Tcl>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]]
}</lang> Example 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).
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. <lang Ursala>#import flo
clop = @iiK0 fleq$-&l+ *EZF ^\~& plus+ sqr~~+ minus~~bbI</lang>
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.
<lang Ursala>#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+-</lang>
test program: <lang Ursala>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</lang> 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))