Re: Need code to calculate the Lower Envelope for a set of (non collinear) points.

*To*: mathgroup at smc.vnet.net*Subject*: [mg51493] Re: [mg51409] Need code to calculate the Lower Envelope for a set of (non collinear) points.*From*: DrBob <drbob at bigfoot.com>*Date*: Wed, 20 Oct 2004 01:21:14 -0400 (EDT)*References*: <79344CBB6297454F86DD3106ADCC76CC82EF83@mail3.nwfwmd.state.fl.us> <opsf2qbuawiz9bcq@monster.cox-internet.com> <opsf3h8v1tiz9bcq@monster.cox-internet.com>*Reply-to*: drbob at bigfoot.com*Sender*: owner-wri-mathgroup at wolfram.com

And here's the fastest (without ConvexHull) solution yet... if you don't mind using Real arithmetic to test the points: Clear[iterate, inverseSlope] inverseSlope[a_, b_] := Compile[{{x, _Real, 1}}, (x[[1]] - a)/(x[[2]] - b)] iterate@{} := {} iterate@{one_} := (Sow@one; {}) iterate[{{a_, b_}, rest__}] := ( Sow@{a, b}; Drop[{Null, rest}, First@Ordering[inverseSlope[a, b] /@ {rest}, -1]] ) Timing@First@Last@Reap@FixedPoint[iterate, a] {0. Second, {{4, 2}, {16, 5}, {44, 13}, { 92, 31}, {688, 269}, {992, 421}, {1000, 491}}} Bobby On Mon, 18 Oct 2004 21:23:09 -0500, DrBob <drbob at bigfoot.com> wrote: > The ConvexHull solutions depend on ConvexHull returning its result in counterclockwise order, and that seems to be undocumented. Can we count on it? > > So here's a solution without ConvexHull: > > rule = {x___, {a_, b_}, {c_, d_}, {e_, f_}, y___} /; > (d - b)/(c - a) > (f - b)/(e - a) :> > {x, {a, b}, {e, f}, y}; > Timing[a //. rule] > MultipleListPlot[{Last[%], a}, PlotJoined -> True, > PlotStyle -> {Hue[0.6]}] > > {0.047 Second, {{4, 2}, {8, 3}, {16, 5}, {44, 13}, {92, 31}, > {688, 269}, {992, 421}, {1000, 491}}} > > ">=" rather than ">" gives a correct result with fewer points: > > rule = {x___, {a_, b_}, {c_, d_}, {e_, f_}, y___} /; > (d - b)/(c - a) >= (f - b)/(e - a) :> > {x, {a, b}, {e, f}, y}; > Timing[a //. rule] > MultipleListPlot[{Last[%], a}, PlotJoined -> True, > PlotStyle -> {Hue[0.6]}] > > {0.047 Second, {{4, 2}, {16, 5}, {44, 13}, {92, 31}, > {688, 269}, {992, 421}, {1000, 491}}} > > Here's a method using FixedPoint: > > iterate@{lowerEnvelope_, {}} := {lowerEnvelope, {}} > iterate@{lowerEnvelope_, {a_}} := {Append[lowerEnvelope, a], {}} > iterate@{lowerEnvelope_, data_} := Module[{one = First@data}, > {Append[lowerEnvelope, one], Drop[data, First@Ordering[ > Divide @@ (# - one) & /@ Rest@data, -1]]} > ] > Timing@First@FixedPoint[iterate, {{}, a}, 9] > > {0.031 Second, {{4, 2}, {16, 5}, {44, 13}, > {92, 31}, {688, 269}, {992, 421}, {1000, 491}}} > > Finally, here are three almost identical solutions using Sow and Reap: > > Clear@iterate > iterate@{} := {} > iterate@{a_} := (Sow@a; {}) > iterate[data_] := Module[{one = First@data}, > Sow@one; > Drop[data, First@Ordering[Divide @@ (# - one) & /@ Rest@data, -1]] > ] > Timing@First@Last@Reap@FixedPoint[iterate, a] > > {0.016 Second, {{4, 2}, {16, > 5}, {44, 13}, {92, 31}, {688, 269}, {992, 421}, {1000, 491}}} > > Clear@iterate > iterate@{} := {} > iterate@{one_} := (Sow@one; {}) > iterate[data : {one_, rest__}] := (Sow@one; > Drop[data, First@Ordering[Divide @@ (# - one) & /@ {rest}, -1]] > ) > Timing@First@Last@Reap@FixedPoint[iterate, a] > > {0.016 Second, {{4, 2}, {16, > 5}, {44, 13}, {92, 31}, {688, 269}, {992, 421}, {1000, 491}}} > > Clear@iterate > iterate@{}:={} > iterate@{one_}:=(Sow@one;{}) > iterate[{one_,rest__}]:=(Sow@one; > Drop[{rest},-1+First@Ordering[Divide@@(#-one)&/@{rest},-1]] > ) > Timing@First@Last@Reap@FixedPoint[iterate,a] > > {0.016 Second,{{4,2},{16,5},{44,13},{92,31},{688,269},{992,421},{1000,491}}} > > I like the last solution best of all these, but ConvexHull still gives a faster solution. > > Bobby > > On Mon, 18 Oct 2004 11:20:08 -0500, DrBob <drbob at bigfoot.com> wrote: > >> Just for fun, here's a slightly different solution (for sorted data): >> >> LowerEnvelope[a_] := Module[{hull = ConvexHull@a}, >> a[[hull /. {{x___, 1, y___, Length@a, z___} :> {1, y, >> Length@a}, {x___, Length@a, y___, 1, z___} :> {1, z, x, Length@a}}]] >> ] >> >> It may be faster, as it uses /. once, rather than twice. >> >> Bobby >> >> On Mon, 18 Oct 2004 09:46:36 -0400, Gilmar Rodriguez <Gilmar.Rodriguez at nwfwmd.state.fl.us> wrote: >> >>> Dear Dr. Bob: >>> Your Lower Envelope Module works like a charm! >>> Intuitively I knew that using the Convex Hull would facilitate finding >>> the Lower Envelope, but my attempts to modify the code of the Convex >>> Hull program always failed. Your solution is both compact and clever! I >>> love it! Thank you!!!! >>> Very grateful to you: >>> Gilmar >>> >>> -----Original Message----- >>> From: DrBob [mailto:drbob at bigfoot.com] To: mathgroup at smc.vnet.net >>> Sent: Saturday, October 16, 2004 7:22 PM >>> To: Gilmar Rodriguez; mathgroup at smc.vnet.net >>> Subject: [mg51493] Re: [mg51409] Need code to calculate the Lower Envelope for a >>> set of (non collinear) points. >>> >>> I think this does it if the data is already sorted: >>> >>> LowerEnvelope[a_]:=Module[{hull=ConvexHull@a}, >>> a[[hull/.{x___,1,y___}:>{1,y,x} >>> /.{x__,Length@a,y___}:>{x,Length@a}]] >>> ] >>> LowerEnvelope@a >>> MultipleListPlot[{%,a},PlotJoined\[Rule]True,PlotStyle\[Rule]{Hue[.6]}] >>> >>> {{4,2},{8,3},{16,5},{44,13},{92,31},{688,269},{992,421},{1000,491}} >>> >>> If the data might NOT be pre-sorted, the function would be: >>> >>> LowerEnvelope[a_] := Module[{data = Sort@a, hull}, >>> hull = ConvexHull@data; >>> data[[hull /. {x___, 1, y___} :> {1, y, x} >>> /. {x__, Length@a, y___} :> {x, Length@a}]] >>> ] >>> >>> Testing with a random reordering of a still gives the same answer: >>> >>> LowerEnvelope[Sort[a, Random[] < 0.5 &]] >>> MultipleListPlot[{%, a}, PlotJoined -> True, PlotStyle -> {Hue[.6]}] >>> >>> {{4, 2}, {8, 3}, {16, 5}, {44, 13}, {92, 31}, {688, 269}, {992, 421}, >>> {1000, 491}} >>> >>> Bobby >>> >>> On Sat, 16 Oct 2004 04:20:45 -0400 (EDT), Gilmar Rodr?guez Pierluissi >>> <gilmar.rodriguez at nwfwmd.state.fl.us> wrote: >>> >>>> Dear Mathematica Solver Group: >>>> >>>> I'm looking for a program to calculate the Lower Envelope >>>> ("LE" for short)for a set of (non-collinear) points on the >>>> plane. Please download the following file containing an >>>> arbitrary set of such points(which I'm calling "A") by >>>> double-clicking the following shortcut: >>>> >>>> http://gilmarlily.netfirms.com/download/A.txt >>>> >>>> You can also down load: >>>> >>>> http://gilmarlily.netfirms.com/download/Lower Envelope.nb >>>> >>>> to evaluate the following steps using Mathematica (version 5): >>>> >>>> If your define your working directory as C:\Temporary, then kindly >>>> evaluate the following Mathematica commands: >>>> >>>> In[1]: A = ReadList["C:\\Temporary\\A.txt", {Number, Number}]; >>>> >>>> Next, plot the set A: >>>> >>>> In[2]: plt1 = ListPlot[A, PlotJoined -> True, PlotStyle -> {Hue[.7]}] >>>> >>>> I can manipulate the program ConvexHull to find the Lower Envelope >>>> for the set A as follows: >>>> >>>> In[3]: << DiscreteMath`ComputationalGeometry` >>>> >>>> In[4]: convexhull = ConvexHull[A] >>>> >>>> The following input gives a picture of the Convex Hull: >>>> >>>> In[5]: plt2 = ListPlot[Table[A[[convexhull[[i]]]], {i, >>>> 1, Length[convexhull]}], PlotJoined -> True, PlotStyle -> >>> {Hue[.6]}] >>>> >>>> Modifying the starting value of index i in In[5] above >>>> (starting at i=96 instead of i=1) gives a picture of the >>>> Lower Envelope of A: >>>> >>>> In[6]:plt3 = ListPlot[Table[A[[convexhull[[i]]]], {i, >>>> 96, Length[convexhull]}], PlotJoined -> True, PlotStyle -> >>> {Hue[.6]}] >>>> >>>> In[7]: Show[plt1, plt3] >>>> >>>> The Lower Envelope of A ("LEA" for short) is given by: >>>> >>>> In[8]: LEA = Table[A[[convexhull[[i]]]], {i, 96, Length[convexhull]}] >>>> >>>> So my question is: How can the code of the ConvexHull program be >>> modified, >>>> to get a program that calculates the LE of a set? >>>> >>>> The following (clumsy)alternative attempt: >>>> >>>> In[9]: LE[B_] := Module[{M}, {L = Length[B]; M = {}; >>>> AppendTo[M, B[[L]]]; {Xg, Yg} = B[[L]]; Do[If[B[[ >>>> L - i + 1]][[2]] < Yg, {{Xg, Yg} = B[[L - >>>> i + 1]], AppendTo[M, B[[L - i + 1]]]}], {i, 1, L - 1}]}; >>> Sort[M]] >>>> >>>> In[10]: LE[A] >>>> >>>> In[11]: plt4 = ListPlot[LE[A], PlotJoined -> True, PlotStyle -> >>> {Hue[.1]}] >>>> >>>> In[12]: Show[plt1, plt4] >>>> >>>> gives me something that is "not even close, and no cigar". >>>> >>>> What I need is an algorithm that gives me the Lower Envelope E of A >>>> as shown in In[8] above. Thank you for your help! >>>> >>>> >>>> >>>> >>> >>> >>> >> >> >> > > > -- DrBob at bigfoot.com www.eclecticdreams.net