Jarvis March in Clojure

Jarvis March in Clojure 57

Posted by Uncle Bob Wed, 12 Aug 2009 02:10:21 GMT

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.

Posted in Uncle Bob's Blatherings
Meta 57 comments, permalink, rss, atom

Comments

Leave a response

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. See http://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.

http://github.com/swannodette/convex-hull/tree/master
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
me 5 months later:

Would it be possible to use a quadtree to eliminate some points from comparison?
Kooba Handbags 8 months later:

Living without an aim is like sailing without a compass. with a new http://www.handbags4buy.com/ idea is a crank until the idea succeeds.
moncler clearance 8 months later:

Very quietly I take my leave.To seek a dream in http://www.edhardy-buy.com/ starlight.
han 8 months later:

http://www.ipadvideoconverters.biz iPad Video Converter iPad Video Converter is a great software for iPad lovers to convert videos to iPad with super fast speed and high quality. Its easy-to-use interface makes iPad to videos conversion routine very simple. And also it can keep the original quality of video files.
liding rode 11 months later:

OK, I got it. Thanks for your sharing. That is very interesting Smile I love reading and I am always searching for informative information like this.
cheap vps about 1 year later:

cheap VPS
supplynflshop about 1 year later:

good
best blu ray converter about 1 year later:

Blu Ray Converter can easily convert blu ray dvd/common dvd/DVD IFO file/DVD folder/DVD ISO file to most common Video, HD Video, Flash Video, iPod, iPhone, iPad, Apple TV, PSP/PS3, Creative Zen, BlackBerry, Zune, Xbox, Mobile Devices, Archos, Common Audio, Supported Portable Devices.

Editing funcionts like trim and crop DVD/ Blu-Ray clips, customize profile list and different watermarks?
Designer Bags about 1 year later:

Nice topic!Thanks for sharing! It really helpful to me about those information.
Dog bark collar about 1 year later:

@David Nolen: Thanks for the link. It has some useful information.
Pandora about 1 year later:

I suppose the students who take our TDD course could claim to be Object Mentor Certified TDDers; and they’d be right.
pandora about 1 year later:

The beauty of getting on the internet is you’re able to see anything your provides, whereby a shop it could be challenging to kind nevertheless his / her stock to locate what’s right. Moms enjoy gifts which are unique. Try and found your current mother utilizing a reward that is personal, as an example a diamond ring utilizing their brand personalized within it. Platinum bands make for wonderful presents given that they immortalize youridentify within gold. The precious metal diamond necklace and a ring creates a fantastic items to your mother, and also presents which can be usually unnoticed .
soundtrack download about 1 year later:

Thanks for uniqueness, much appreciate this!
stamping parts about 1 year later:

NingBo YinZhou WeiSheng is the profession production in various stamping parts,deep drawn parts,welding parts,machining parts, and metal semi-processed goods and finished product.
Criminal Check about 1 year later:

Great stuff, worth reading. Thanks for this awesome information.
hermes replica watches about 1 year later:

if they don’t go ????? go if they feel they need to. hermes h bracelet chestnut to. They may not look at the hermes birkin replica the casket. Let them lead the way.”Loder cheap hermes birkin bag way.”Loder was so worried about her own hermes h bracelet replica own daughter’s classmates’ reaction to the funeral 40cm birkin funeral that she arranged a closed coffin 40cm birkin coffin with photos of her children laid hermes h bracelet yellow laid on top.The Loders didn’t want to replica bags hermes to make the
ghd australia about 1 year later:

Growth hormone deficiency in Australia, teach you how to protect your hair curl and style hair salon in any way to rectify.
Tenant Screening about 1 year later:

Ct Credit Bureau offers full tenant screening services to property managers, landlords, and others in the real estate and rental industry.
Criminal Records about 1 year later:

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.
mechanical Seal about 1 year later:

This article is very usefull for me! I can see that you are putting a lots of efforts into your blog. I will keep watching in your blog, thanks.
mobile screening equipment about 1 year later:

KEESTRACK?manufacturing expert of mobile crusher and mobile screening equipment. Company engaged in the research & development, manufacture, sales of track screening machine,mobile screening,cone crusher,jaw crusher,impact crusher,mobile screening equipment,mobile crushing equipment, till now we are proud to say we are at front of the industry.
ssara about 1 year later:

The concept is great. Hope it is well understood. childcare jobs
Designer Sunglasses about 1 year later:

Inspired ,MEN designer Sunglasses Women Replica Sunglass at cheap discount price
okey oyunu oyna about 1 year later:

very good.

internette görüntülü olarak okey oyunu oyna, gerçek kisilerle tanis, turnuva heyecanini yasa.
okey oyunu indir about 1 year later:

I like this comment, fers full tenant screening services to property managers.

Farkl? ki?ilerle online okey oyunu oynayabilece?iniz bu heyecan verici oyunu hemen bilgisayarlar?n?za indirin. Okey oyunu indir, bilgisayar?na kur ve oyanamaya hemen ba?la. okeyoyunuindir.com Thank you…
dentist in san antonio about 1 year later:

Your algorithm is right. But you have to modify it. You can use branch method. It is more effective. You can divide your work in two or three branch.
Jewellery about 1 year later:

Online UK costume and fashion jewellery shop with,
thailand about 1 year later:

Thank you for this information. I’ve been looking for something like this for quite a while. Keep up the good work, cheers! klimat thailand lången tips teknik blogg
engelska bildelar about 1 year later:

Thank you for this information. I’ve been looking for something like this for quite a while. Keep up the good work, cheers! reservdelar land rover reservdelar mg reservdelar mini
beats by dre store over 2 years later:

something like this for quite a while. Keep up the good work, cheers!high quality headphones new design headphones
bagsupplyer over 2 years later:

I really enjoy your articles. Very well written. New fashion women Gucci Purses with high quality for wholesale on line from China free shipping
african Mango dr oz over 2 years later:

It is not hard to understand why the question of value is still raised when you consider that today?s outsourcing models have some inherent limitations that reduce the overall gains companies can achieve .
Tips For Bowling over 2 years later:

As fathers commonly go, it is seldom a misfortune to be fatherless; and considering the general run of sons, as seldom a misfortune to be childless.
ysbearing/yisong@1stbearing.com over 2 years later:

Slewing bearing called slewing ring bearings, is a comprehensive load to bear a large bearing, can bear large axial, radial load and overturning moment. http://www.1stbearing.com
christian louboutin over 2 years later:
The professional design make you foot more comfortable. Even more tantalizing,this pattern make your legs look as long as you can,it will make you looked more attractive.Moveover,it has reasonable price.If you are a popular woman,do not miss it.

Technical details of Christian Louboutin Velours Scrunch Suede Boots Coffee:
```
Color: Coffee
Material: Suede
4(100mm) heel
Signature red sole x
```
Fashion, delicate, luxurious Christian louboutins shoes on sale, one of its series is Christian Louboutin Tall Boots, is urbanism collocation. This Christian louboutins shoes design makes people new and refreshing. Red soles shoes is personality, your charm will be wonderful performance.
Brisbane Limo Hire over 2 years later:

something like this for quite a while. Keep up the good work, cheers. Brisbane Limo Hire
Brisbane Limo Hire over 2 years later:

something like this for quite a while. Keep up the good work, cheers.Brisbane Limo Hire
tv guide listings over 2 years later:

Thanks for a good share. This is some quality articles you are posting right here. white beach
web design peterborough over 2 years later:

This site is very informative..Thanks for sharing this site with us..!
lilly pulitzer clearance over 2 years later:

Great share! Hope you will give us more of this kind! bubble skirt bubble skirt pattern how to make a bubble skirt bubble skirt tutorial
väder thailand over 3 years later:

This site is the best. I hope that more info of this kind will get to us soon.
bladeless fans over 3 years later:

Jarvis March in Clojure 52 good post54
louboutin sales over 3 years later:

Jarvis March in Clojure 53 hoo,good article!!I like the post!165
Injection mold over 3 years later:

Intertech Machinery Inc. provides the most precise Plastic Injection Mold and Rubber Molds from Taiwan. With applying excellent unscrewing device in molds,

Intertech is also very professional for making flip top Cap Molds in the world. Mold making is the core business of Intertech (Taiwan). With world level technology, Intertech enjoys a very good reputation for making Injection Mold and Plastic Molds for their worldwide customers.
hermes pink over 3 years later:

When you purchase a laptop, invariably you will opt for a simple leather bag to tote it around?
Mold Making over 3 years later:

With more than 20 years of experience, Intertech provides an extensive integrated operational ability from design to production of molds 100% made in Taiwan. Additional to our own mold making factory, we also cooperate with our team vendors to form a very strong working force in Taiwan.

For the overseas market, we work very closely with local representatives in order to take care of the technical communication and after-sales service to our customers. We also participate in the EUROMOLD & FAKUMA exhibitions and meet our customers every year in Europe. By concentrating on mold “niche markets”, we play a very useful mold maker role from the Far East whenever customers want to develop their new projects. We provide services from A to Z to our customers on a very economic cost and effect basis.

Mentor	twitter id
Uncle Bob	unclebobmartin
Brett Schuchert	schuchert
Michael Feathers	mfeathers
Bob Koss	bob_koss

Jarvis March in Clojure 57

Comments

Blog Search

Follow us on twitter

Categories

Blogroll

Syndicate