Embedded Files
Posted by Uncle Bob on 08/11/2009
OK, all you Clojure gurus, I need your help. I need to speed this algorithm up.
Just for fun I thought I’d write the Jarvis March algorithm in Clojure. This algorithm is for finding the convex hull of a collection of points. The basic idea is really simple. Imagine that all the points in the collection are represented by nails in a board. Find the left-most nail and tie a string to it. Then wind the string around the nails. The string will only touch nails that are part of the convex hull.
The details are not really that much more difficult. You start at the left-most point and calculate the angle from vertical required to touch every other point. The point with the minimum angle is the next point. You keep going around looking finding points with the minimum angle that is greater than the previous angle. When you get back to the starting point you are done.
Calculating angles can be time consuming, so I use a psuedo-angle algorithm. It doesn’t calculate the actual angle, rather it is a function that increases with the true angle, and goes from [0, 4).
The code is pretty simple.
(ns convexHull.convex-hull (:use clojure.contrib.math)) (defn quadrant-one-pseudo-angle [dx dy] (/ dx (+ dy dx))) (defn pseudo-angle [[dx dy]] (cond (and (= dx 0) (= dy 0)) 0 (and (>= dx 0) (> dy 0)) (quadrant-one-pseudo-angle dx dy) (and (> dx 0) (<= dy 0)) (+ 1 (quadrant-one-pseudo-angle (abs dy) dx)) (and (<= dx 0) (< dy 0)) (+ 2 (quadrant-one-pseudo-angle (abs dx) (abs dy))) (and (< dx 0) (>= dy 0)) (+ 3 (quadrant-one-pseudo-angle dy (abs dx))) :else nil)) (defn point-min [[x1 y1 :as p1] [x2 y2 :as p2]] (cond (< x1 x2) p1 (= x1 x2) (if (< y1 y2) p1 p2) :else p2)) (defn find-min-point [points] (reduce point-min points)) (defn delta-point [[x1 y1] [x2 y2]] [(- x1 x2) (- y1 y2)]) (defn angle-and-point [point base] [(pseudo-angle (delta-point point base)) point]) (defn min-angle-and-point [ap1 ap2] (if (< (first ap1) (first ap2)) ap1 ap2)) (defn find-point-with-least-angle-from [base angle points] (reduce min-angle-and-point (remove #(< (first %) angle) (map #(angle-and-point % base) (remove (fn [p] (= base p)) points))))) (defn hull [points] (println "Start") (let [starting-point (find-min-point points)] (println starting-point) (loop [hull-list [starting-point] angle 0 last-point starting-point] (let [[angle next-point] (find-point-with-least-angle-from last-point angle points)] (if (= next-point (first hull-list)) hull-list (recur (conj hull-list next-point) angle next-point))))))
I execute it with this:
(ns convexHull.time-hull (:use convexHull.convex-hull)) (def r (java.util.Random.)) (defn rands [] (repeatedly #(.nextGaussian r))) (defn points [] (take 400000 (partition 2 (rands)))) (let [hull-points (time (hull (points)))] (printf "Points: %d\n" (count hull-points)) (doseq [x hull-points] (println x)))
This takes way too long to run. The equivalent java program will do a million points in half a second. This one is taking 25 seconds to do a half-million points. It won’t even do a million points. It slows way way down and then runs out of memory. (There must be some kind of disk caching going on or something.)
Anyway, I’d be interested in seeing how a real Clojure programmer would speed this program up.
Comments
Chouser about 13 hours later:
I’m no performance guru, but on a sequence that large (0.5 or 1 million items), you might do better with a vector instead, at least on recent HEAD versions of Clojure with chunked seq support.
(defn points [] (vec (take 400000 (partition 2 (rands)))))
This means points will no longer return a lazy seq, but now ‘reduce’, ‘remove’ and ‘map’ will go in chunks. I’m seeing about a 26% improvement with that change.
Leonel about 14 hours later:
I’m no Clojure expert, but here are a few quick tips:
First, function points returns a lazy seq, so keep in mind that counting it means exausting it, which will account for a few seconds.
Second, function rands relies on reflection to invoke method nextGaussian on the randomizer. Use (set! *warn-on-reflection true) and the Clojure compiler will alert you when reflection is taking place. Seehttp://clojure.org/java_interop#typehints for how to use Type Hints to avoid reflection.
Alternatively, use r as a local variable of function rands. I don’t know how much faster it is than looking at a global var, though, but at least it eliminates the reflection call.
(defn rands[] (let [r (java.util.Random.)] (repeatedly #(.nextGaussian r))))
Finally, use into-array to consume the seq before using it.
(let [hull-points (time (hull (into-array (points))))] (printf “Points: %d\n” (count hull-points)) (doseq [x hull-points] (println x)))
With all this combined, I reduced the time from 36s to 29s (on my box). It’s not yet the half-second we get in pure Java, though.
This is for now. I’ll give it a deeper look later.
LPetit about 18 hours later:
Hi Uncle Bob,
Would you mind sharing with us the java code you’re talking about, so that we exactly know what we are comparing ?
Thanks,
— Laurent
Jonathan Smith 1 day later:
Hi Bob,
I am by no means a Clojure Guru, but I have come up with a few simple improvements that make it about 2-3x faster.
The most important thing I did was to remove the ‘delta-point’ function. Essentially the problem is that the delta point function causes you to create 2x as many vectors as you would have otherwise.
I suspect that if there is a way to remove/reduce the data-structure creating properties of ‘angle and point’ as well, you would see even better performance from the code.
I also incorporated some of the changes from the people posting above and changed the ‘point-min’ function to avoid Destructuring, which gets us another second or so.
(I did some miscellaneous inlining using macros as well, but I believe Java would have done that inlining by itself anyway…)
code is here: http://github.com/JonathanSmith/misc/tree/master
Under ‘march.clj’
Anyway, some things to think about. :-)
-Jon
Uncle Bob 1 day later:
I wrote the Java code a couple of years ago. Here it is:
import java.util.*; public class JarvisMarch { Points pts; private Points hullPoints = null; private List<Double> hy; private List<Double> hx; private int startingPoint; private double currentAngle; private static final double MAX_ANGLE = 4; public JarvisMarch(Points pts) { this.pts = pts; } /** * The Jarvis March, sometimes known as the Gift Wrap Algorithm. * The next point is the point with the next largest angle. * <p/> * Imagine wrapping a string around a set of nails in a board. Tie the string to the leftmost nail * and hold the string vertical. Now move the string clockwise until you hit the next, then the next, then * the next. When the string is vertical again, you will have found the hull. */ public int calculateHull() { initializeHull(); startingPoint = getStartingPoint(); currentAngle = 0; addToHull(startingPoint); for (int p = getNextPoint(startingPoint); p != startingPoint; p = getNextPoint(p)) addToHull(p); buildHullPoints(); return hullPoints.x.length; } public int getStartingPoint() { return pts.startingPoint(); } private int getNextPoint(int p) { double minAngle = MAX_ANGLE; int minP = startingPoint; for (int i = 0; i < pts.x.length; i++) { if (i != p) { double thisAngle = relativeAngle(i, p); if (thisAngle >= currentAngle && thisAngle <= minAngle) { minP = i; minAngle = thisAngle; } } } currentAngle = minAngle; return minP; } private double relativeAngle(int i, int p) {return pseudoAngle(pts.x[i] - pts.x[p], pts.y[i] - pts.y[p]);} private void initializeHull() { hx = new LinkedList<Double>(); hy = new LinkedList<Double>(); } private void buildHullPoints() { double[] ax = new double[hx.size()]; double[] ay = new double[hy.size()]; int n = 0; for (Iterator<Double> ix = hx.iterator(); ix.hasNext();) ax[n++] = ix.next(); n = 0; for (Iterator<Double> iy = hy.iterator(); iy.hasNext();) ay[n++] = iy.next(); hullPoints = new Points(ax, ay); } private void addToHull(int p) { hx.add(pts.x[p]); hy.add(pts.y[p]); } /** * The PseudoAngle is a number that increases as the angle from vertical increases. * The current implementation has the maximum pseudo angle < 4. The pseudo angle for each quadrant is 1. * The algorithm is very simple. It just finds where the angle intesects a square and measures the * perimeter of the square at that point. The math is in my Sept '06 notebook. UncleBob. */ public static double pseudoAngle(double dx, double dy) { if (dx >= 0 && dy >= 0) return quadrantOnePseudoAngle(dx, dy); if (dx >= 0 && dy < 0) return 1 + quadrantOnePseudoAngle(Math.abs(dy), dx); if (dx < 0 && dy < 0) return 2 + quadrantOnePseudoAngle(Math.abs(dx), Math.abs(dy)); if (dx < 0 && dy >= 0) return 3 + quadrantOnePseudoAngle(dy, Math.abs(dx)); throw new Error("Impossible"); } public static double quadrantOnePseudoAngle(double dx, double dy) { return dx / (dy + dx); } public Points getHullPoints() { return hullPoints; } public static class Points { public double x[]; public double y[]; public Points(double[] x, double[] y) { this.x = x; this.y = y; } // The starting point is the point with the lowest X // With ties going to the lowest Y. This guarantees // that the next point over is clockwise. int startingPoint() { double minY = y[0]; double minX = x[0]; int iMin = 0; for (int i = 1; i < x.length; i++) { if (x[i] < minX) { minX = x[i]; iMin = i; } else if (minX == x[i] && y[i] < minY) { minY = y[i]; iMin = i; } } return iMin; } } }
And here’s the driver code.
import java.util.Random; public class TimeJarvisMarch { public static void main(String[] args) { for (int i=0; i<10; i++) test(); } private static void test() { final int POINTS = 1000000; double x[] = new double[POINTS]; double y[] = new double[POINTS]; Random r = new Random(); for (int i=0; i<POINTS; i++) { x[i] = r.nextGaussian(); y[i] = r.nextGaussian(); } JarvisMarch.Points pts = new JarvisMarch.Points(x,y); JarvisMarch jm = new JarvisMarch(pts); double start = System.currentTimeMillis(); int n = jm.calculateHull(); double end = System.currentTimeMillis(); System.out.printf("%d points found %d vertices %f seconds\n", POINTS, n, (end-start)/1000.); } }
tbrownaw 1 day later:
Hmm, using into-array on the points and timing that separately shows a huge slowdown and out-of-memory while generating the point list.
(let [points-arr (time (into-array (points))) hull-points (time (hull points-arr))] (printf "Points: %d\n" (count hull-points)) (doseq [x hull-points] (println x)))
With 321k points it OOMs (with no visible disk activity), with 320k it takes about 17 seconds, with 200k it takes about 7 seconds, with 100k it takes maybe 4 seconds.
There also seems to be about a 40% speedup of the actual calculation by replacing find-point-with-least-angle-from with a version that doesn’t nest so many functions (kinda contrary to the purpose of functional languages, I suppose)
(defn find-point-with-least-angle-from [base angle points] (reduce (fn [accum p] (if (= p base) accum (let [curr (angle-and-point p base)] (cond (< (first curr) angle) accum (< (first accum) (first curr)) accum :else curr)) )) [5 ()] points))
Both versions of this appear to only slow down linearly with point count, but even a single call on 50k points takes more than a half second.
This is all with everything in one file run with “clojure -i ”, the docs I see don’t appear to indicate that being significantly different than pre-compiled (which was taking far to long to figure out how to do).
David Nolen 1 day later:
I have have a version that runs in about 2-2.5s on my machine. This works by:
1. Inlining almost everything with macros.
2. Using javax.vecmath.Vector2d
3. Using primitive arrays of javax.vecmath.Vector2d
4. Aggressive type hinting.
The one thing that I note that could be made faster is that find-point-with-least-angle-from could be made parallel. This gets called around 20 times for 400,000 random points. This would be trivial to implement with agents and you could basically half the time spent here for each core.
David Nolen 1 day later:
I’d also like to point out there is a Clojure version of this program that has exactly the same performance characteristics of the Java program.
(.test (new TimeJarvisMarch))
;)
David Nolen 1 day later:
github.com/swannodette/convex-hull/tree/master
Here’s an optimized version that run in about 2s-2.5s for 400,000 pts, and ~5.6s for 1,000,000.
leonardo 2 days later:
In D language too, plus comparisons: leonardo-m.livejournal.com/87388.html
David Nolen 3 days later:
I just updated my optimized version of this program. It can do 400,000 points in ~550ms on my machine. It probably can’t be made much faster than this. It was fun figuring out how to use macros to inline much of the logic.
I have a hunch that javac does some interesting magic to optimize this particular kind of code for the JVM, the kind of thing that the Clojure compiler could eventually do, but will take time.
github.com/swannodette/convex-hull/tree/master
Page updated
Google Sites
Report abuse