Skip to main content
replaced http://mathematica.stackexchange.com/ with https://mathematica.stackexchange.com/
Source Link

I suggest to use linked lists and recursionlinked 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.

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.

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.

Source Link
Leonid Shifrin
  • 114.9k
  • 16
  • 333
  • 424

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