rule-based implementation of an algorithm

Question

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?

Answer 1

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