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

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

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: [mg51487] 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