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