Closest Pair, Part 1

February 3, 2015

We will represent points using a pair, with the x coordinate in the car and the y coordinate in the cdr. We begin with a function to generate a random list of points, for testing:

(define max-coord 1000000)

(define (rand-points n)
  (define (r) (randint max-coord))
  (do ((n n (- n 1))
       (ps (list) (cons (cons (r) (r)) ps)))
      ((zero? n) ps)))

Now the algorithm computes the distance between each pair of points, keeping track of the minimum:

(define (dist p q)
  (define (square x) (* x x))
  (sqrt (+ (square (- (car p) (car q)))
           (square (- (cdr p) (cdr q))))))

(define (closest-pair ps)
  (let ((min-dist (* max-coord max-coord))
        (min-pair (list)))
    (do ((ps ps (cdr ps))) ((null? (cdr ps)))
      (do ((qs (cdr ps) (cdr qs))) ((null? qs))
        (let ((d (dist (car ps) (car qs))))
          (when (< d min-dist)
            (set! min-dist d)
            (set! min-pair (cons (car ps) (list (car qs))))))))
    (values min-dist min-pair)))

We used random numbers from the Standard Prelude. You can run the program at http://programmingpraxis.codepad.org/a8jmQeXn, where you will also see some examples. We will see a solution for the closest pair problem that takes time O(n log n) in the next exercise.

Pages: 1 2

11 Responses to “Closest Pair, Part 1”

  1. matthew said

    Here’s one I made earlier:

    https://github.com/matthewarcus/misc/blob/master/closest.cpp

    Divide and conquer, splitting region alternately horizontally and vertically (not sure if that’s really a benefit over always splitting the same way, but it’s more fun).

  2. use strict;
    use warnings;
    
    sub closest {
      my @m;
      foreach my $i (0..(@_-1)) {
        foreach my $j (0..($i-1)) {
          my $d = ($_[$i][0]-$_[$j][0])**2 + ($_[$i][1]-$_[$j][1])**2;
          @m = ($d,$i,$j) unless @m && $m[0] < $d;
        }
      }
      return @m;
    }
    
    my @p = map { [rand,rand] } 1..100; # Generate a random list of co-ordinates..
    my ($d2,$p1,$p2) = closest( @p );    # Get details
    printf "(%f,%f) - (%f,%f) - %f\n", @{$p[$p1]},@{$p[$p2]}, sqrt $d2; # Render them!
    
    
  3. jeff said

    my attempt in python:

    import random
    import math
    
    
    def closest(amt, pts):
    
        dlist = []
        for i in range(amt - 1):
            pt = pts.pop(0)
            x = pt[0]
            y = pt[1]
    
            for q in range(len(pts)):
                xlen = x - pts[q][0]
                ylen = y - pts[q][1]
                distance = math.sqrt((xlen ** 2) + (ylen ** 2))
                dlist.append([distance, i + 1, q + 2 + i])
    
        for i, d in enumerate(dlist):
            print "Point {} to {} = {}".format(d[1], d[2], d[0])
    
        dlist.sort()
    
        msg = "\nClosest points are {} and {} \nDistance = {}"
        print msg.format(dlist[0][1], dlist[0][2], dlist[0][0])
    
    
    def pointmaker(amt):
    
        pts = []
    
        for i in range(amt):
            x = round(random.uniform(0, 100), 2)
            y = round(random.uniform(0, 100), 2)
    
            point = (x, y)
            pts.append(point)
    
        for i, pt in enumerate(pts):
            print "Point {} - {}".format(i + 1, pt)
        print""
        closest(amt, pts)
    
    pointmaker(8)
    

    >>>
    Point 1 – (3.44, 55.26)
    Point 2 – (1.13, 95.38)
    Point 3 – (96.25, 6.12)
    Point 4 – (51.36, 25.47)
    Point 5 – (98.23, 12.24)
    Point 6 – (63.23, 13.53)
    Point 7 – (9.95, 20.72)
    Point 8 – (21.43, 99.51)

    Point 1 to 2 = 40.1864467203
    Point 1 to 3 = 105.016359202
    Point 1 to 4 = 56.4249102791
    Point 1 to 5 = 104.095458595
    Point 1 to 6 = 72.9125297874
    Point 1 to 7 = 35.1481393533
    Point 1 to 8 = 47.7671707347
    Point 2 to 3 = 130.442178761
    Point 2 to 4 = 86.0840345244
    Point 2 to 5 = 127.830628568
    Point 2 to 6 = 102.741581164
    Point 2 to 7 = 75.1791726477
    Point 2 to 8 = 20.7158610731
    Point 3 to 4 = 48.8828661189
    Point 3 to 5 = 6.43232461867
    Point 3 to 6 = 33.8412248596
    Point 3 to 7 = 87.5262817672
    Point 3 to 8 = 119.665051289
    Point 4 to 5 = 48.7014352971
    Point 4 to 6 = 16.8362852197
    Point 4 to 7 = 41.6815378795
    Point 4 to 8 = 79.8606692935
    Point 5 to 6 = 35.0237647891
    Point 5 to 7 = 88.6863506973
    Point 5 to 8 = 116.250990964
    Point 6 to 7 = 53.7629472778
    Point 6 to 8 = 95.6023033195
    Point 7 to 8 = 79.6219473512

    Closest points are 3 and 5
    Distance = 6.43232461867
    >>>

  4. Mike said

    Brute force version in Python:

    def closest_pair(points):
        """returns the distance between a closest pair of points and
           a pair of such points.
           """
        
        min_dd = float("+inf")
        min_pair = None
    
        for p, q in combinations(points, 2):
            dd = (q[0] - p[0])**2 + (q[1] - p[1])**2
            if dd < min_dd:
                min_dd = dd
                min_pairs= (p, q)
    
        return sqrt(min_dd), min_pair
    

    A modified version that returns a list of all pairs that have the minimum distance.

    def all_closest_pairs(points):
        """returns the distance between the closest pair of points, and
           a list of such point pairs.
           """
        
        min_dd = float("+inf")
        min_pairs = []
    
        for p, q in combinations(points, 2):
            dd = (q[0] - p[0])**2 + (q[1] - p[1])**2
            if dd < min_dd:
                min_dd = dd
                min_pairs = [(p, q)]
    
            elif dd == min_dd:
                min_pairs.append((p, q))
    
        return sqrt(min_dd), min_pairs
    
  5. import math
    def closestPoints(n):
    	#sort the points based on x coordinate.
    	#compare distances between consecutive points.
    	#sort points based on y coordinate.
    	#compare distances between consecutive points.
    	#take the smallest of the two.
    	factor = 1.618
    	step = int(len(n)/factor)
    	while step > 0:
    		for i in range(0,len(n)-step):
    			j = i
    			while j+step < len(n) and n[j+step][0] < n[j][0]:
    				n[j+step][0],n[j][0],n[j+step][1],n[j][1] = (n[j][0],
    				n[j+step][0],n[j][1],n[j+step][1])
    				j+= step
    		step = int(step/factor)
    	minx = 10**100
    	for i in range(0,len(n)-1):
    		if math.sqrt((n[i][0]-n[i+1][0])**2 +(n[i][1]-n[i+1][1])**2) < minx:
    			minx =  math.sqrt((n[i][0]-n[i+1][0])**2 +(n[i][1]-n[i+1][1])**2)
    
    	while step > 0:
    		for i in range(0,len(n)-step):
    			j = i
    			while j+step < len(n) and n[j+step][1] < n[j][1]:
    				n[j+step][0],n[j][0],n[j+step][1],n[j][1] = (n[j][0],
    				n[j+step][0],n[j][1],n[j+step][1])
    				j+= step
    		step = int(step/factor)
    	miny = 10**100
    	for i in range(0,len(n)-1):
    		if math.sqrt((n[i][0]-n[i+1][0])**2 +(n[i][1]-n[i+1][1])**2) < minx:
    			miny =  math.sqrt((n[i][0]-n[i+1][0])**2 +(n[i][1]-n[i+1][1])**2)
    	if minx < miny:
    		return minx
    	else:
    		return miny
    
  6. small error, i have < minx in the second run, should be < miny.

  7. Paul said

    Here the brute force method in Python.

    from math import sqrt
    from random import uniform as uni
    
    def closest_pair(points):
         d, p, q = min(((p[0] - q[0]) ** 2 + (p[1] - q[1]) ** 2, p, q)
                   for i, p in enumerate(points)
                       for q in points[i+1:])
         return (sqrt(d), p, q)
    
    def random_points(n, xmin=0, xmax=1, ymin=0, ymax=1):
        return [(uni(xmin, xmax), uni(ymin, ymax)) for _ in xrange(n)]
    
    print closest_pair(random_points(1000))
    
  8. Ken said

    Brute force in C# done very, very, very Sloppy

    Cat Class I don’t know why I called it that

    class Cat
    {
    List PointsList = new List();
    double Shortest_Distance = 10000;
    string ClosestPair;

    public void CreatePoints()
    {
    int i = 0;

    try
    {
    Console.WriteLine(“Add point (X,Y)”);
    string Input;
    do
    {
    Input = Console.ReadLine();

    if (Input != “End”)
    {
    string[] NewPointData = Input.Split(‘,’);
    PointsList.Add(new Point(Convert.ToInt32(NewPointData[0]), Convert.ToInt32(NewPointData[1])));
    i++;
    }
    } while (Input != “End”);

    Search();
    }
    catch (Exception ex)
    {
    Console.WriteLine(ex);
    }
    }

    public void Search()
    {
    for (int i = 0; i < PointsList.Count; i++)
    {
    for (int j = 0; j < PointsList.Count; j++)
    {
    if (i != j)
    {
    float DeltaX = PointsList[j].Get_X – PointsList[i].Get_X;
    float DeltaY = PointsList[j].Get_Y – PointsList[i].Get_Y;

    if (DeltaY < 0)
    DeltaY *= -1;
    if (DeltaX < 0)
    DeltaX *= -1;

    double TempCsqr = Math.Sqrt((Math.Pow(DeltaX, 2) + Math.Pow(DeltaY, 2)));

    Console.WriteLine("Current: " + TempCsqr);

    if (TempCsqr < Shortest_Distance)
    {
    Shortest_Distance = TempCsqr;
    ClosestPair = "(" + PointsList[j].Get_String + ")(" + PointsList[i].Get_String + ")";
    }
    }
    Console.WriteLine ("Shortest: " + Shortest_Distance + "\n Closest Pair is" + ClosestPair);
    }
    }
    }
    }

    Points Class:

    namespace Points
    {
    class Point
    {
    int X, Y;

    public Point(int _X, int _Y)
    {
    X = _X;
    Y = _Y;
    }

    public int Get_X
    {
    get { return X; }
    }

    public int Get_Y
    {
    get { return Y; }
    }

    public string Get_String
    {
    get { return X.ToString() + "," + Y.ToString(); }
    }
    }
    }

    INPUT:

    Add point (X,Y)
    1,2
    2,1
    1,5
    1,9
    6,5
    End

    OUTPUT:

    Shortest: 10000
    Closest Pair is

    Current: 1.4142135623731
    Shortest: 1.4142135623731
    Closest Pair is(2,1)(1,2)

    Current: 3
    Shortest: 1.4142135623731
    Closest Pair is(2,1)(1,2)

    Current: 7
    Shortest: 1.4142135623731
    Closest Pair is(2,1)(1,2)

    Current: 5.8309518948453
    Shortest: 1.4142135623731
    Closest Pair is(2,1)(1,2)

    Current: 1.4142135623731
    Shortest: 1.4142135623731
    Closest Pair is(2,1)(1,2)
    Shortest: 1.4142135623731
    Closest Pair is(2,1)(1,2)

    Current: 4.12310562561766
    Shortest: 1.4142135623731
    Closest Pair is(2,1)(1,2)

    Current: 8.06225774829855
    Shortest: 1.4142135623731
    Closest Pair is(2,1)(1,2)

    Current: 5.65685424949238
    Shortest: 1.4142135623731
    Closest Pair is(2,1)(1,2)

    Current: 3
    Shortest: 1.4142135623731
    Closest Pair is(2,1)(1,2)

    Current: 4.12310562561766
    Shortest: 1.4142135623731
    Closest Pair is(2,1)(1,2)
    Shortest: 1.4142135623731
    Closest Pair is(2,1)(1,2)

    Current: 4
    Shortest: 1.4142135623731
    Closest Pair is(2,1)(1,2)

    Current: 5
    Shortest: 1.4142135623731
    Closest Pair is(2,1)(1,2)

    Current: 7
    Shortest: 1.4142135623731
    Closest Pair is(2,1)(1,2)

    Current: 8.06225774829855
    Shortest: 1.4142135623731
    Closest Pair is(2,1)(1,2)

    Current: 4
    Shortest: 1.4142135623731
    Closest Pair is(2,1)(1,2)
    Shortest: 1.4142135623731
    Closest Pair is(2,1)(1,2)

    Current: 6.40312423743285
    Shortest: 1.4142135623731
    Closest Pair is(2,1)(1,2)

    Current: 5.8309518948453
    Shortest: 1.4142135623731
    Closest Pair is(2,1)(1,2)

    Current: 5.65685424949238
    Shortest: 1.4142135623731
    Closest Pair is(2,1)(1,2)

    Current: 5
    Shortest: 1.4142135623731
    Closest Pair is(2,1)(1,2)

    Current: 6.40312423743285
    Shortest: 1.4142135623731
    Closest Pair is(2,1)(1,2)
    Shortest: 1.4142135623731
    Closest Pair is(2,1)(1,2)

  9. Claire said

    C++ Version:
    #include
    #include
    #include
    #include
    using namespace std;

    #define ArraySize 4

    struct Point
    {
    Point(double xx, double yy):x(xx),y(yy){}

    static double distance(const Point &p1, const Point &p2)
    {
    double dx = p2.x – p1.x;
    double dy = p2.y – p1.y;
    return sqrt(dx*dx + dy*dy);
    }
    double x;
    double y;
    };

    void CalcDistance(vector pointsInput, vector& points, int size)
    {
    int num = 0;
    int min = 0,p1 = 0, p2 = 0;
    for( int i = 0; i< size – 1; i++)
    {
    for ( int j = i+1; j< size; j++)
    {
    points.push_back(Point::distance(pointsInput[i],pointsInput[j]));
    cout<<"point "<<i+1 <<" and Point "<<j+1 <<" = " <<points.at(num) <<endl;
    if(num == 0 )
    {
    min = points.at(0);
    }
    else if( points.at(num) < min )
    {
    min = points.at(num);
    p1 = i + 1;
    p2 = j + 1;
    }
    num++;
    }
    }

    cout<<"The closest pair are "<<"point "<<p1 <<" and point "<<p2<<endl;
    cout<< "The distance is "<< min << endl;
    }

    int main(int argc, char *argv[])
    {
    vector myPoints;
    int count = 0;
    double x;
    double y;

    while( count > x >> y;
    Point point(x,y);
    myPoints.push_back(point);
    count++;
    }

    vector output;
    CalcDistance(myPoints,output,ArraySize);
    return 0;
    }

  10. Claire said

    The comment is really weird, some pieces of my codes are missing. Here is the new code link http://codepad.org/mNXrJb9z

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: