# rule-based implementation of an algorithm

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?

-
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
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"] –  Guess who it is. 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

### 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.

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;


(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[]]]]]] *)

(* {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[pts_ll] :=
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];


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
*)

-
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