new files

Task/Closest-pair-problem/0DESCRIPTION

@@ -0,0 +1,63 @@
+{{Wikipedia|Closest pair of points problem}}
+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 [[wp:Closest pair of points problem|Closest pair of points problem]] in the ''planar'' case.
+
+The straightforward solution is a O(n<sup>2</sup>) 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 &lt; 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&amp;conquer approach, as explained also at [[wp:Closest pair of points problem#Planar_case|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⌉)<sub>x</sub>
+   yL ← { p ∈ yP : p<sub>x</sub> ≤ xm }
+   yR ← { p ∈ yP : p<sub>x</sub> &gt; xm }
+   (dL, pairL) ← ''closestPair'' of (xL, yL)
+   (dR, pairR) ← ''closestPair'' of (xR, yR)
+   (dmin, pairMin) ← (dR, pairR)
+   '''if''' dL &lt; dR '''then'''
+     (dmin, pairMin) ← (dL, pairL)
+   '''endif'''
+   yS ← { p ∈ yP : |xm - p<sub>x</sub>| &lt; 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)<sub>y</sub> - yS(i)<sub>y</sub> &lt; dmin
+       '''if''' |yS(k) - yS(i)| &lt; 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'''
+* [[wp:Closest pair of points problem|Closest pair of points problem]]
+* [http://www.cs.mcgill.ca/~cs251/ClosestPair/ClosestPairDQ.html Closest Pair (McGill)]
+* [http://www.cs.ucsb.edu/~suri/cs235/ClosestPair.pdf Closest Pair (UCSB)]
+* [http://classes.cec.wustl.edu/~cse241/handouts/closestpair.pdf Closest pair (WUStL)]
+* [http://www.cs.iupui.edu/~xkzou/teaching/CS580/Divide-and-conquer-closestPair.ppt Closest pair (IUPUI)]

Task/Closest-pair-problem/1META.yaml

@@ -0,0 +1,2 @@
+---
+note: Classic CS problems and programs

Task/Closest-pair-problem/Ada/closest-pair-problem.ada

@@ -0,0 +1,67 @@
+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;

Task/Closest-pair-problem/Clojure/closest-pair-problem.clj

@@ -0,0 +1,36 @@
+(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])))))

Task/Closest-pair-problem/Go/closest-pair-problem-1.go

@@ -0,0 +1,46 @@
+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
+}

Task/Closest-pair-problem/Go/closest-pair-problem-2.go

@@ -0,0 +1,117 @@
+// 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
+}

Task/Closest-pair-problem/Haskell/closest-pair-problem-1.hs

@@ -0,0 +1,16 @@
+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

Task/Closest-pair-problem/Haskell/closest-pair-problem-2.hs

@@ -0,0 +1,3 @@
+*Main> testCP
+([[0.8347201880148426,0.40774840545089647],[0.8348731214261784,0.4087113189531284]],9.749825850154334e-4)
+(4.02 secs, 488869056 bytes)

Task/Closest-pair-problem/Java/closest-pair-problem.java

@@ -0,0 +1,186 @@
+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");
+  }
+}

Task/Closest-pair-problem/JavaScript/closest-pair-problem.js

@@ -0,0 +1,27 @@
+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
+    };
+  }
+}

Task/Closest-pair-problem/Perl/closest-pair-problem.pl

@@ -0,0 +1,107 @@
+#! /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";

Task/Closest-pair-problem/PicoLisp/closest-pair-problem.l

@@ -0,0 +1,14 @@
+(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)) ) ) )

Task/Closest-pair-problem/Python/closest-pair-problem.py

@@ -0,0 +1,87 @@
+"""
+  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()

Task/Closest-pair-problem/R/closest-pair-problem-1.r

@@ -0,0 +1,30 @@
+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)
+    }
+}

Task/Closest-pair-problem/R/closest-pair-problem-2.r

@@ -0,0 +1,20 @@
+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]])
+  }

Task/Closest-pair-problem/R/closest-pair-problem-3.r

@@ -0,0 +1,18 @@
+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)
+     }

Task/Closest-pair-problem/R/closest-pair-problem-4.r

@@ -0,0 +1,36 @@
+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

Task/Closest-pair-problem/R/closest-pair-problem-5.r

@@ -0,0 +1,88 @@
+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))

Task/Closest-pair-problem/REXX/closest-pair-problem.rexx

@@ -0,0 +1,37 @@
+/*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

Task/Closest-pair-problem/Ruby/closest-pair-problem.rb

@@ -0,0 +1,68 @@
+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

Task/Closest-pair-problem/Scala/closest-pair-problem.scala

@@ -0,0 +1,24 @@
+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)
+   }
+}

Task/Closest-pair-problem/Tcl/closest-pair-problem.tcl

@@ -0,0 +1,90 @@
+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]]
+}