Closest Pair, Part 2
February 13, 2015
; closest pair part 2
(define (take n xs)
(let loop ((n n) (xs xs) (ys '()))
(if (or (zero? n) (null? xs))
(reverse ys)
(loop (- n 1) (cdr xs)
(cons (car xs) ys)))))
(define (drop n xs)
(let loop ((n n) (xs xs))
(if (or (zero? n) (null? xs)) xs
(loop (- n 1) (cdr xs)))))
(define-syntax fold-of
(syntax-rules (range in is)
((_ "z" f b e) (set! b (f b e)))
((_ "z" f b e (v range fst pst stp) c ...)
(let* ((x fst) (p pst) (s stp)
(le? (if (positive? s) =)))
(do ((v x (+ v s))) ((le? p v) b)
(fold-of "z" f b e c ...))))
((_ "z" f b e (v range fst pst) c ...)
(let* ((x fst) (p pst) (s (if (< x p) 1 -1)))
(fold-of "z" f b e (v range x p s) c ...)))
((_ "z" f b e (v range pst) c ...)
(fold-of "z" f b e (v range 0 pst) c ...))
((_ "z" f b e (x in xs) c ...)
(do ((t xs (cdr t))) ((null? t) b)
(let ((x (car t)))
(fold-of "z" f b e c ...))))
((_ "z" f b e (x is y) c ...)
(let ((x y)) (fold-of "z" f b e c ...)))
((_ "z" f b e p? c ...)
(if p? (fold-of "z" f b e c ...)))
((_ f i e c ...)
(let ((b i)) (fold-of "z" f b e c ...)))))
(define-syntax list-of (syntax-rules ()
((_ arg ...) (reverse (fold-of
(lambda (d a) (cons a d)) '() arg ...)))))
(define sort #f)
(define merge #f)
(let ()
(define dosort
(lambda (pred? ls n)
(if (= n 1)
(list (car ls))
(let ((i (quotient n 2)))
(domerge pred?
(dosort pred? ls i)
(dosort pred? (list-tail ls i) (- n i)))))))
(define domerge
(lambda (pred? l1 l2)
(cond
((null? l1) l2)
((null? l2) l1)
((pred? (car l2) (car l1))
(cons (car l2) (domerge pred? l1 (cdr l2))))
(else (cons (car l1) (domerge pred? (cdr l1) l2))))))
(set! sort
(lambda (pred? l)
(if (null? l) l (dosort pred? l (length l)))))
(set! merge
(lambda (pred? l1 l2)
(domerge pred? l1 l2))))
(define rand #f)
(define randint #f)
(let ((two31 #x80000000) (a (make-vector 56 -1)) (fptr #f))
(define (mod-diff x y) (modulo (- x y) two31)) ; generic version
; (define (mod-diff x y) (logand (- x y) #x7FFFFFFF)) ; fast version
(define (flip-cycle)
(do ((ii 1 (+ ii 1)) (jj 32 (+ jj 1))) ((< 55 jj))
(vector-set! a ii (mod-diff (vector-ref a ii) (vector-ref a jj))))
(do ((ii 25 (+ ii 1)) (jj 1 (+ jj 1))) ((< 55 ii))
(vector-set! a ii (mod-diff (vector-ref a ii) (vector-ref a jj))))
(set! fptr 54) (vector-ref a 55))
(define (init-rand seed)
(let* ((seed (mod-diff seed 0)) (prev seed) (next 1))
(vector-set! a 55 prev)
(do ((i 21 (modulo (+ i 21) 55))) ((zero? i))
(vector-set! a i next) (set! next (mod-diff prev next))
(set! seed (+ (quotient seed 2) (if (odd? seed) #x40000000 0)))
(set! next (mod-diff next seed)) (set! prev (vector-ref a i)))
(flip-cycle) (flip-cycle) (flip-cycle) (flip-cycle) (flip-cycle)))
(define (next-rand)
(if (negative? (vector-ref a fptr)) (flip-cycle)
(let ((next (vector-ref a fptr))) (set! fptr (- fptr 1)) next)))
(define (unif-rand m)
(let ((t (- two31 (modulo two31 m))))
(let loop ((r (next-rand)))
(if (list a)))
((eq? (car seed) 'set) (set! fptr (caadr seed))
(set! a (list->vector (cdadr seed))))
(else (/ (init-rand (modulo (numerator
(inexact->exact (car seed))) two31)) two31)))))
(set! randint (lambda args
(cond ((null? (cdr args))
(if (< (car args) two31) (unif-rand (car args))
(floor (* (next-rand) (car args)))))
((< (car args) (cadr args))
(let ((span (- (cadr args) (car args))))
(+ (car args)
(if (< span two31) (unif-rand span)
(floor (* (next-rand) span))))))
(else (let ((span (- (car args) (cadr args))))
(- (car args)
(if (< span two31) (unif-rand span)
(floor (* (next-rand) span))))))))))
(define max-coord 1000000)
(define infinity (* max-coord max-coord))
(define (rand-points n)
(define (r) (randint max-coord))
(do ((n n (- n 1))
(ps (list) (cons (make-rectangular (r) (r)) ps)))
((zero? n) ps)))
(define (brute-force-closest-pair ps)
(let ((min-dist infinity) (min-pair (list)))
(do ((ps ps (cdr ps))) ((null? (cdr ps)))
(do ((qs (cdr ps) (cdr qs))) ((null? qs))
(let ((d (magnitude (- (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)))
(define (closest-pair-aux xs ys)
(if (< (length xs) 10) (brute-force-closest-pair xs)
(let* ((len (length xs)) (mid (quotient len 2))
(x-left (take mid xs)) (x-right (drop mid xs))
(y-left (list-of p (p in ys)
(< (real-part p) (real-part (car x-right))))))
(y-right (list-of p (p in ys)
(>= (real-part p) (real-part (car x-right))))))
(let-values (((left-dist left-pair) (closest-pair-aux x-left y-left))
((right-dist right-pair) (closest-pair-aux x-right y-right)))
(let ((min-dist (if (< left-dist right-dist) left-dist right-dist))
(min-pair (if (< left-dist right-dist) left-pair right-pair)))
(let ((close-points (list-of p (p in ys)
(< (magnitude (- p (car x-right))) min-dist))))
(if (< (length close-points) 2) (values min-dist min-pair)
(let-values (((close-dist close-pair)
(brute-force-closest-pair close-points)))
(if (< min-dist close-dist)
(values min-dist min-pair)
(values close-dist close-pair))))))))))
(define (closest-pair ps)
(closest-pair-aux
(sort (lambda (a b) (< (real-part a) (real-part b))) ps)
(sort (lambda (a b) (< (imag-part a) (imag-part b))) ps)))
(define ps (rand-points 1000))
(time (call-with-values
(lambda () (brute-force-closest-pair ps))
(lambda (min-dist min-pair)
(display min-dist) (display " ")
(display min-pair) (newline))))
(time (call-with-values
(lambda () (closest-pair ps))
(lambda (min-dist min-pair)
(display min-dist) (display " ")
(display min-pair) (newline))))
Like this:
Like Loading...
Related
Pages: 1 2 3
Posted by programmingpraxis
Programming Praxis 
I was jumping the gun a bit with my solution to Part I.
Do you use the fact that the ys list is sorted?
from collections import namedtuple from itertools import combinations, product class Point(namedtuple("Point_", "x y")): def mag(self): return (self.x ** 2 + self.y ** 2) ** 0.5 def __sub__(self, other): return Point(self.x - other.x, self.y - other.y) def closest_pair2(points): def cp(ps): if len(ps) <= 3: d, pr = min(((p - q).mag(), (p, q)) for p,q in combinations(ps, 2)) else: m = len(ps) // 2 mid = ps[m] d, pr = min(cp(ps[:m]), cp(ps[m:])) for p in ps[m-1::-1]: if p.x < mid.x - d: break for q in ps[m:]: if q.x > p.x + d: break if p.y - d <= q.y <= p.y + d: td = (p - q).mag() if td < d: d = td pr = p, q return d, pr return cp(sorted(points))Takes about 4.5 seconds to find closest pair in a set of 1000000 points–much faster than brute-force method.
Example:
>>> from random import randrange >>> from timeit import timeit >>> def random_point(a, b=None): return Point(randrange(a,b), randrange(a,b)) >>> points = set() >>> while len(points) < 1000000: points.add(random_point(5000000)) >>> timeit(stmt="print(closest_pair2(points))", setup="from __main__ import points, closest_pair2", number=1) (6.4031242374328485, (Point(x=185209, y=328746), Point(x=185214, y=328750))) 4.447064363699084 >>>