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