Mathematica 9 is now available
Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2004
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2004

[Date Index] [Thread Index] [Author Index]

Search the Archive

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


  • Prev by Date: Re: LegendreP (Symbolic) is different in Mathematica5 than previous versions (M4, M3 ..)
  • Next by Date: Extrapolation in mathematica
  • Previous by thread: Re: Need code to calculate the Lower Envelope for a set of (non collinear) points.
  • Next by thread: Re: Need code to calculate the Lower Envelope for a set of (non collinear) points.