Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

When I first started learning about rule-based programming with Mathematica, I tried to translate this algorithm for computing the convex hull of a set of 2-D points in $O(n \log(n))$ time, to use rule-based replacement. For convenience, I'll paste the pseudocode of the algorithm below:

Input: a list P of points in the plane.

  • Sort the points of P by x-coordinate (in case of a tie, sort by y-coordinate).

  • Initialize U and L as empty lists. The lists will hold the vertices of upper and lower hulls respectively.

  • for i = 1, 2, ..., n:

    while L contains at least two points and the sequence of last two points of L and the point P[i] does not make a counter-clockwise turn:
    remove the last point from L
    append P[i] to L

  • for i = n, n-1, ..., 1:

    while U contains at least two points and the sequence of last two points of U and the point P[i] does not make a counter-clockwise turn:
    remove the last point from U
    append P[i] to U

  • Remove the last point of each list (it's the same as the first point of the other list).

  • Concatenate L and U to obtain the convex hull of P. Points in the result will be listed in counter-clockwise order.

(Actually, I used a version of this algorithm from the book "Computational Geometry" by de Berg & others, but it's pretty much the same thing except for the "counter-clockwise turn condition" that uses determinant instead of cross-product, which isn't really relevant to this question.)

This was my implementation, which has the same end result as the stated algorithm, but as you can see doesn't operate quite the same way:

convexHullPM[pts_] := Module[{li, isRightTurn}, 
  isRightTurn[p1_, p2_, p3_] := 
     Sign[Det@MapThread[Prepend, {{p1, p2, p3}, {1, 1, 1}}]] == -1;
  li = Sort[pts];
  li = Join[li, Reverse@Most@li];
  li //. {pre___, a_List, b_List, c_List, post___} :>
         {pre, a, c, post} /; Not[isRightTurn[a, b, c]]]

It works, giving the same sequence of points as Mathematica's own ConvexHull function from the Computational Geometry Package.

However it's slow (referring to its time complexity, as the input size increases); timings with increasingly larger sets of input points seem to indicate $O(n^2)$ complexity, which seems about right - because the way I constructed the rule in Mathematica, it starts matching from the start every time a replacement is made. Whereas the listed algorithm has the notion of a current point, and points are removed from the set (if need be) from the most recent points.

So what would be a better rule-based implementation?

share|improve this question
    
Probably you need ReplacePart and FixedPoint. –  swish May 25 '13 at 10:33
    
@swish ReplacePart won't help, since it copies the entire list even when only one or a few points are modified. Besides, it does not delete points, so Delete will likely be more appropriate (but will have the same problem). –  Leonid Shifrin May 25 '13 at 13:45
1  
A few pointers: 1. you can compactly implement your test, like so (though I have already incorporated the negation): notRightTurn = NonNegative[Det[PadLeft[{##}, {3, 3}, 1]]] &;; 2. Here's a shorter way to generate li: n = Length[pts]; li = ArrayPad[Sort[pts], {{0, n - 1}}, "Reflected"] –  J. M. May 25 '13 at 15:14
    
@J.M. Thanks. Mathematica seems to have more functions than I have brain cells! –  Aky May 25 '13 at 15:22

1 Answer 1

up vote 10 down vote accepted

Preamble

I suggest to use linked lists and recursion. Note that this won't be the fastest possible solution, but it will be idiomatic, rule-based, will have the right asymptotic complexity, and will have decent performance. I want to stress that I am here after these features, rather than raw performance, as those seems to be at the core of your question.

Linked lists

Here is the linked list API we will need:

ClearAll[ll, toLinkedList, fromLinkedList];
SetAttributes[ll, HoldAllComplete];
toLinkedList[l_List] := Fold[ll[#2, #1] &, ll[], Reverse@l];
fromLinkedList[lst_ll] := List @@ Flatten[lst, Infinity, ll];

This has been explained in more details in the linked post. We will need an additional function to reverse a linked list. One possible (idiomatic) implementation can be:

ClearAll[llrev];
llrev[lst_ll] := llrev[lst, ll[]];
llrev[ll[], accum_] := accum;
llrev[ll[head_, tail_], accum_] := llrev[tail, ll[head, accum]];

(of course, one could also use Composition[toLinkedList, Reverse, fromLinkedList], which would also be somewhat faster, but I like the recursive version).

Some simple examples how this is used:

tst = toLinkedList[Range[5]]

(* ll[1, ll[2, ll[3, ll[4, ll[5, ll[]]]]]] *)

fromLinkedList[tst]

(* {1, 2, 3, 4, 5} *)

llrev[tst]

(* ll[5, ll[4, ll[3, ll[2, ll[1, ll[]]]]]] *)

Implementation

Now we are ready to rewrite your algorithm. One subtlety here is that we'd either need to consider 4 adjacent points in this approach, or use FixedPoint (many runs). I would pick the FixedPoint route. Here is the code:

ClearAll[isRightTurn];
isRightTurn[p1_, p2_, p3_] := 
   Sign[Det@MapThread[Prepend, {{p1, p2, p3}, {1, 1, 1}}]] == -1;

and the main function:

ClearAll[convexHullLL];
convexHullLL[pts_List] :=
   convexHullLL[toLinkedList[Join[#, Reverse@Most@#] &[Sort[pts]]]];

convexHullLL[pts_ll] :=
  fromLinkedList@FixedPoint[
    llrev[convexHullLL[#, ll[]]] &,
    pts
  ];

convexHullLL[ll[], accum_] := accum; 
convexHullLL[ll[a_, ll[b_, tail : ll[c_, _]]], accum_]/;Not[isRightTurn[a, b, c]] :=
   convexHullLL[ll[a, tail], accum];
convexHullLL[ll[head_, tail_], accum_] := convexHullLL[tail, ll[head, accum]];

The main recursive body here is typical for the solutions based on linked lists, and the ideas behind this technique were explained in more details in the linked post.

Examples and benchmarks

We can test that it gives the same results as your code:

rpts = RandomInteger[{-10,10},{100,2}];
convexHullPM[rpts]==convexHullLL[rpts]

(* True *)

But it has a different complexity:

rptsLarge = RandomReal[{-5,5},{2*10^3,2}];
(cnv1 = convexHullPM[rptsLarge]);//AbsoluteTiming
(cnv2 = convexHullLL[rptsLarge]);//AbsoluteTiming
cnv1==cnv2

 (*
    {2.935547,Null}
    {0.111328,Null}
    True
 *)
share|improve this answer
    
Thanks for your reply, Leonid. I'll go through it carefully, along with your linked list discussion, when my mind is fresh. You're right that my question is more about idiomaticity (if that's even a word) rather than raw performance. –  Aky May 25 '13 at 14:12
    
@Aky Sure, was glad to help. Good if it works for you. –  Leonid Shifrin May 25 '13 at 20:22

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.