PRESENTATION

last edited on June 04, 2010 12:34 PM by admin

Typeset

%hide
#### Max and Min methods

def find_ymax(ls1, ls2):
#Maekake a check that tests if they are empty
Q=[]
W=[]
for i in range(len(ls1)):
g=ls1[i]
Q.append(g[1])
for i in range(len(ls2)):
g=ls2[i]
W.append(g[1])
ymax1 = max(Q)
ymax2 = max(W)
return max([ymax1, ymax2])

def find_ymin(ls1, ls2):
#Maekake a check that tests if they are empty
Q=[]
W=[]
for i in range(len(ls1)):
g=ls1[i]
Q.append(g[1])
for i in range(len(ls2)):
g=ls2[i]
W.append(g[1])
ymin1 = min(Q)
ymin2 = min(W)
return min(ymin1, ymin2)

def find_xmax(ls1, ls2):
#Maekake a check that tests if they are empty
Q=[]
W=[]
for i in range(len(ls1)):
g=ls1[i]
Q.append(g[0])
for i in range(len(ls2)):
g=ls2[i]
W.append(g[0])
xmax1 = max(Q)
xmax2 = max(W)
return max(xmax1, xmax2)

def find_xmin(ls1, ls2):
#Make a check that tests if they are empty
Q=[]
W=[]
for i in range(len(ls1)):
g=ls1[i]
Q.append(g[0])
for i in range(len(ls2)):
g=ls2[i]
W.append(g[0])
xmin1 = min(Q)
xmin2 = min(W)
return min(xmin1, xmin2)

def max_list(ls1, ls2):
S=[]
R=[]
for i in range(0,len(ls1)):
for z in ls1[i]:
S.append(z)
newS = S[0:len(S):2]
for i in range(0,len(ls2)):
for z in ls2[i]:
R.append(z)
newR = R[0:len(R):2]
m=max(newR)
n=max(newS)
return max(m,n)

def min_list(ls1, ls2):
S=[]
R=[]
for i in range(0,len(ls1)):
for z in ls1[i]:
S.append(z)
newS = S[0:len(S):2]
for i in range(0,len(ls2)):
for z in ls2[i]:
R.append(z)
newR = R[0:len(R):2]
m=min(newR)
n=min(newS)
return min(m,n)

#### Area and construction of Polygons

#### Line Construction and Methods

def draw_line(pt1, pt2):
"""Returns a line as a list [a,b,c] for ay+bx=c""" #Andy => I updated it with our new kind of line..
if pt2[0]-pt1[0] != 0:
m = (pt2[1] - pt1[1])/(pt2[0]- pt1[0])
b = pt1[1] -m*(pt1[0])
newline = line()
newline.a = 1
newline.b = -m
newline.c = b
return newline
else:
m = (pt2[0] - pt1[0])/(pt2[1]- pt1[1])
b = pt1[0] - m*(pt1[1])
newline = line()
newline.a = -m
newline.b = 1
newline.c = b
return newline

###We should be careful about the datatypes we allow in
### What is line here stored as!?!?! ([a,b,c] list)
### Is this a test to intersect two lines?
### A 'does_question' method should, in general, return a boolean
### This method appears to return an intersection point so that seems like it is something that MUST exist for all
### cases.

def max_list(ls1, ls2):
S=[]
R=[]
for i in range(0,len(ls1)):
for z in ls1[i]:
S.append(z)
newS = S[0:len(S):2]
for i in range(0,len(ls2)):
for z in ls2[i]:
R.append(z)
newR = R[0:len(R):2]
m=max(newR)
n=max(newS)
return max(m,n)

def min_list(ls1, ls2):
S=[]
R=[]
for i in range(0,len(ls1)):
for z in ls1[i]:
S.append(z)
newS = S[0:len(S):2]
for i in range(0,len(ls2)):
for z in ls2[i]:
R.append(z)
newR = R[0:len(R):2]
m=min(newR)
n=min(newS)
return min(m,n)

#### Area and construction of Polygons

#### Line Construction and Methods

def draw_line(pt1, pt2):
"""Returns a line as a list [a,b,c] for ax+by=c""" #Andy => I updated it with our new kind of line..
if pt2[0]-pt1[0] != 0:
m = (pt2[1] - pt1[1])/(pt2[0]- pt1[0])
b = pt1[1] -m*(pt1[0])
newline = line()
newline.a = -m
newline.b = 1
newline.c = b
return newline
else:
m = (pt2[0] - pt1[0])/(pt2[1]- pt1[1])
b = pt1[0] - m*(pt1[1])
newline = line()
newline.a = 1
newline.b = -m
newline.c = b
return newline

###We should be careful about the datatypes we allow in
### What is line here stored as!?!?! ([a,b,c] list)
### Is this a test to intersect two lines?
### A 'does_question' method should, in general, return a boolean
### This method appears to return an intersection point so that seems like it is something that MUST exist for all
### cases.

class pgon:
def __init__(self, list_vertices = []):
self.vertex_list = list_vertices
self.area = self.PArea()
self.xcent = self.cent_x()
self.ycent = self.cent_y()
self.centroid = (self.xcent, self.ycent)

##Is the 'L' a list or a line?
##And if it is a list, what is it listing?
def PArea(self):
S=[]
for i in range(0,len(self.vertex_list)):
for z in self.vertex_list[i]:
S.append(z)
A=0
for i in range(0,len(S),2):
if i < len(S)-2:
A= A + S[i]*S[i+3] - S[i+2]*S[i+1]
else:
A= A + S[i]*S[1] - S[0]*S[i+1]
return abs((1/2)*A)

### So here it appears that line is a list of form [a,b,c] ?
### The new polygon is the bottom one
def new_pgon(self, current_line):
Q=[]
if current_line.b==0:
for i in range(len(self.vertex_list)):
if i < len(self.vertex_list)-1:
#does_intersect returns intersection if it exists
r = self.does_intersect(self.vertex_list[i],self.vertex_list[i+1],current_line)
if r != None: # meaning there exists intersection
Q.append(r)
g=self.vertex_list[i+1]
if g[0] <= current_line.c:
Q.append(g)
if i == len(self.vertex_list)-1:
r= self.does_intersect(self.vertex_list[i],self.vertex_list[0],current_line)
if r!= None:
Q.append(r)
g=self.vertex_list[0]
if g[1] < current_line.c:
Q.append(g)
N = pgon(Q)
return N
else:
for i in range(len(self.vertex_list)):
if i < len(self.vertex_list)-1:
#does_intersect returns intersection if it exists
r = self.does_intersect(self.vertex_list[i+1],self.vertex_list[i],current_line)
if r != None: # meaning there exists intersection
Q.append(r)
g=self.vertex_list[i+1]
if g[1] < -(current_line.a/current_line.b)*g[0]+(current_line.c/current_line.b):
Q.append(g)
if i == len(self.vertex_list)-1:
r = self.does_intersect(self.vertex_list[0],self.vertex_list[i],current_line)
if r!= None:
Q.append(r)
g=self.vertex_list[0]
if g[1] < -(current_line.a/current_line.b)*g[0]+(current_line.c/current_line.b):
Q.append(g)
N = pgon(Q)
return N

def cent_x(self):
"""Takes a list of Pgon vertices and finds the x-centroid"""
S=[]
for i in range(0,len(self.vertex_list)):
for z in self.vertex_list[i]:
S.append(z)
B=self.area
if B==0:
return
A=0
for i in range(0,len(S),2):
if i < len(S)-2:
A= A + (S[i]+S[i+2])*(S[i]*S[i+3] - S[i+2]*S[i+1])
else:
A= A + (S[i]+S[0])*(S[i]*S[1] - S[0]*S[i+1])
return (1/(6*B))*A

def cent_y(self):
"""Takes a list of PGon vertices and finds the y-centroid"""
S=[]
for i in range(0,len(self.vertex_list)):
for z in self.vertex_list[i]:
S.append(z)
B=self.area
if B==0:
return
A=0
for i in range(0,len(S),2):
if i < len(S)-2:
A= A + (S[i+1]+S[i+3])*(S[i]*S[i+3] - S[i+2]*S[i+1])
else:
A= A + (S[i+1]+S[1])*(S[i]*S[1] - S[0]*S[i+1])
return (1/(6*B))*A

def does_intersect(self, pt1,pt2,current_line):
line1=line(pt1,pt2).give() ###.give() returns a list! [a,b,c]
line2=current_line.give()
m=matrix([line1,line2])
x0=m.solve_right(m)[0,2]
y0=m.solve_right(m)[1,2]
if pt1[0]!=pt2[0]:
if pt1[0]>=pt2[0]:
if pt2[0] <= x0 and pt1[0] >= x0:
return (x0,y0)
if pt1[0]<= pt2[0]:
if pt1[0] <= x0 and pt2[0] >= x0:
return (x0,y0)
else:
if pt1[1] <= pt2[1]:
if pt2[1] >= y0 and pt1[1] <= y0:
return (x0,y0)
if pt1[1] >= pt2[1]:
if pt1[1] >= y0 and pt2[1] <= y0:
return (x0,y0)

class line:
def __init__(self, p1 = (0,0), p2 = (0,0) ):
if p2[0]-p1[0]==0:
self.a=1
self.b=0
self.c=p1[0]
else:
d = ( (p2[1] - p1[1]) / (p2[0] - p1[0]) )
self.a = -d
self.b = 1
self.c = p1[1] - (d*p1[0])

def give(self):
ls = [self.a, self.b, self.c]
return ls

class hamline:
def __init__(self, xmin, xmax, ymin, ymax, resolution, polygons, HALF_AREA):
self.steps = int(1/resolution) #Pass in resolution as a percent (meaning less than 1)!
self.resolution = resolution
#### The xmin, etc will be important to have visible in this scope
self.xmin = xmin
self.xmax = xmax
self.ymin = ymin
self.ymax = ymax
self.polygons = polygons #pass this in a list of two pgons
self.HALF_AREA = HALF_AREA
self.found = False
self.areaLine = self.find_line()
self.topSearch = False

def __repr__(self):
#print "%sy + %sx = %s"%(self.a, self.b, self.c)
return "[%s, %s, %s]"%(self.a, self.b, self.c)

###This should just return a line!!!
def find_line(self):
steps = self.steps
xmax = self.xmax
xmin = self.xmin
ymax = self.ymax
ymin = self.ymin

###These are all just points stored as tuples
corner1 = (xmax, ymax)
corner2 = (xmax, ymin)
corner3 = (xmin, ymin)
corner4 = (xmin, ymax)
divis = [0..steps]
comps = [] #as in components
yrange = ymax-ymin
xrange = xmax-xmin

###The thing is that all of our lines that we are going to use will be generated

###Generate points along the line
### leftline_points etc are lists of 2-tuples
leftline_points = []
rightline_points = []
topline_points = []
bottomline_points = []
####Note: As implemented, we divide the region into steps by steps components.
#### In a case where xrange >> yrange or vice versa, this could be a problem
#### I suggest doing component wise entry of steps in the input box
for i in divis:
leftline_points.append( (xmin, (i*(yrange/steps) + ymin)) )
rightline_points.append( (xmax, (i*(yrange/steps) + ymin)) )
topline_points.append( ( (i*(xrange/steps) + xmin), ymax) )
bottomline_points.append( ( (i*(xrange/steps) + xmin), ymin) )
#### Post-cond: we have all our points from which to make our lovely lines!
#### Check the lines from leftline to rightline
for i in leftline_points:
for j in rightline_points:
connect = line(i,j)
#We have the set of two pgons. Those are invariant. New pgons get generated. Find areas of those:
# Each pgon has a vertex_list that can create a new pgon
# newpgon takes a vertex_list and a line
Poly1 = self.polygons[0].vertex_list
Poly2 = self.polygons[1].vertex_list
P1 = pgon(Poly1)
P2 = pgon(Poly2)
P1_intersect = P1.new_pgon(connect) ## these are new polygons, created by line
P2_intersect = P2.new_pgon(connect)
ALP1= P1_intersect.area
AP1 = P1.area
ALP2 = P2_intersect.area
AP2 = P2.area
#This is assuming that the new pgons both appear below the line
if not( (P1_intersect.area==0) or (P2_intersect.area ==0) ) and not( (P1_intersect.area == P1.area) or (P2_intersect.area == P2.area) ):
if ( (ALP1/AP1 < self.resolution + 1/2) and (ALP1/AP1 > 1/2 - self.resolution) ):
if ( (ALP2/AP2 < self.resolution + 1/2 ) and (ALP2/AP2 > - self.resolution + 1/2) ):
self.found = True
return connect
if (not self.found):
self.topSearch = True
for i in topline_points:
for j in bottomline_points:
connect = line(i,j)
#We have the set of two pgons. Those are invariant. New pgons get generated. #Find areas of those:
# Each pgon has a vertex_list that can create a new pgon
# newpgon takes a vertex_list and a line
Poly1 = self.polygons[0].vertex_list
Poly2 = self.polygons[1].vertex_list
P1 = pgon(Poly1)
P2 = pgon(Poly2)
P1_intersect = P1.new_pgon(connect)
P2_intersect = P2.new_pgon(connect)
ALP1= P1_intersect.area
AP1 = P1.area
ALP2 = P2_intersect.area
AP2 = P2.area
#This is assuming that the new pgons both appear below the line
if ( (self.HALF_AREA - (P1_intersect.area + P2_intersect.area) ) < (self.resolution*2.0*self.HALF_AREA) ):
self.found = True
return connect

def give(self):
ls = [self.a, self.b, self.c]
return ls

#### Testing new Pgon class

#####HEY THIS IS ANDY I AM JUST GOING TO TEST STUFF HERE SO IF YOU HAVE TO MAKE CHANGES PLEASE TEST IN THE ABOVE BOX K THX

@interact
def f(vertices=input_box("(225,37) (208,20) (192,13) (152,23) (122,54) (85,67) (56,65) (52,75) (45,95) (35,105) (25,112) (18,113) (3,105) (4,99) (13,101) (6,97) (18,98) (30,95) (38,80) (38,60) (40,32) (51,23) (52,8) (58,7) (60,18) (65,17) (66,4) (71,4) (74,16) (100,15) (102,3) (110, 2) (109,15) (115, 14) (117, 4) (124, 3) (125, 13) (160, 3) (194, 2) (221, 20)",
"Vertices P1 (separated by spaces)", type=str),
clr=color_selector(Color('purple'), widget='jpicker'), vertices2=input_box("(-10,135) (10,150) (20,165) (20,180) (10,190) (0,190) (-10,180) (-20,190) (-30,190) (-40,180) (-40,165) (-30,150)",
"Vertices P2 (separated by spaces)", type=str),
clr2=color_selector(Color('red'), widget='jpicker'), error=slider(0,.5,default= 0.150300601202405)):
### Are these v just vertices? Are they listed in a specific order? Are they tuples?

v = sage_eval('['+vertices.replace(')','),')+']')
v2 = sage_eval('['+vertices2.replace(')','),')+']')

### Is it possible to be a little more specific about which areas?
### We have to treat v and v2 as pgons now
shape1 = pgon(v)
shape2 = pgon(v2)

A1= shape1.area
Half_A1=A1/2
A2= shape2.area
Half_A2= A2/2

### Are these stored as point tuples?
C1= shape1.centroid
C2= shape2.centroid

### Are these the vertices located on the extrema of x and y?
m1= max_list(v,v2)
n1=min_list(v,v2)

# this is the old code: L=draw_line(C1,C2)
### So what exactly is this L? -- It is a line between centroids :)
#L=hamline(0,0,0,C1,C2).give()
CLine = line(C1, C2)
L= CLine
CP1=shape1.new_pgon(L)
CP2=shape2.new_pgon(L)

###Just to clarify, these are the two bottom areas?
CA1=float(CP1.PArea())
CA2=float(CP2.PArea())
#b1=how_big(CA1,Half_A1)
#b2=how_big(CA2, Half_A2)
bottom = CA1 + CA2
top = (A1-CA1) + (A2-CA2)
HALF_AREA = (A1 +A2)/2.0
resolution = error

#Calculate the x and y max
ymax = find_ymax(v, v2)
ymin = find_ymin(v, v2)
xmax = find_xmax(v, v2)
xmin = find_xmin(v, v2)

#Create hamline
poly_list = [shape1, shape2]
ham = hamline(xmin, xmax, ymin, ymax, resolution, poly_list, HALF_AREA)
areaCutLine = ham.areaLine
L2=areaCutLine.give()
NP1=shape1.new_pgon(ham.areaLine)
NP2=shape2.new_pgon(ham.areaLine)

###Just to clarify, these are the two bottom areas?
BA1=float(NP1.PArea())
BA2=float(NP2.PArea())
#b1=how_big(BA1,Half_A1)
#b2=how_big(BA2, Half_A2)
bottom = BA1 + BA2
top = (A1-BA1) + (A2-BA2)
HALF_AREA = (bottom + top)/2.0

#make it so it gets shown

####Print Statements

####Debug
#print 'ymax = %s , ymin = %s'%(ymax,ymin)
#print 'xmax = %s , xmin = %s'%(xmax,xmin)

if(ham.topSearch):
print "Searching top to bottom"

print "The area of P1 is %f and the area of P2 is %f" %(A1,A2)
print ''
print 'The area of the upper half of P1 is %s'%(A1-BA1)
print 'The area of the lower half of P1 is %s'%(BA1)
print ''
print 'The area of the upper half of P2 is %s'%(A2-BA2)
print 'The area of the lower half of P2 is %s'%(BA2)

#show(polygon(v, color=clr)+polygon(v2, color=clr2) + plot(-(L[0]/L[1])*x+L[2]/L[1],n1-1,m1+1, thickness=2, color='black')+plot(-#(L2[0]/L2[1])*x+L2[2]/L2[1],n1-1,m1+1, thickness=2, color='black') + polygon(NP1, color='green')+ polygon(NP2, color= 'yellow') )

###Yo Matt, I need you to add a plot of areaCutLine because theoretically that should be our line!
#show( point([C1,C2],color='black', pointsize=30) + polygon(v, color=clr)+polygon(v2, color=clr2) + plot(-(L[0]/L[1])*x+L[2]/L[1],n1-1,m1+1, #thickness=2, color='black')+plot(-(L2[0]/L2[1])*x+L2[2]/L2[1],n1-1,m1+1, thickness=2, color='black')+polygon(NP1, color='green')+ polygon(NP2, #color=# 'yellow') )

P1Plot = polygon(v, color=clr)
P2Plot = polygon(v2, color=clr2)
CLinePlot = plot(-(L.a/L.b)*x+L.c/L.a,n1-1,m1+1, thickness=.3, color='black')
HAMLINEPlot= plot(-(L2[0]/L2[1])*x+L2[2]/L2[1],n1-1,m1+1, thickness=2, color='black')
#HAMLINEPlot= line([(L2[2],ymax+1), (L2[2], ymin -1)])
P1BPLOT = polygon(NP1.vertex_list, color='green')
P2BPLOT = polygon(NP2.vertex_list, color= 'yellow')

show(P1Plot+ P2Plot +P1BPLOT + P2BPLOT + HAMLINEPlot +CLinePlot)
show('Our Hamline Approximation is: y=%s'%(-(L2[0]/L2[1])*x+L2[2]/L2[1]))
show( 'UP1/BP1= %s'%((BA1)/A1))
show('UP2/BP2 = %s' %((BA2)/A2))

THE HAM SANDWICH THEOREM

Andy Barr, Chloe Huber, Lauren Boubel, Maddie Romansic, Matt Junge

SO... What is the Ham Sandwich Theorem?

With the invention of the three ingredient sandwich by John Montagu, the Fourth Earl of Sandwich, came the desire to divide sandwiches equitably among friends.

In sandwich terms the Ham Sandwich Theorem states that:

"Given a sandwich composed of bread, cheese and meat arranged as haphazardly as you please, there exists a slice that equally divides the volume of all three ingredients."

For those of you who prefer more rigorous statements or are unaware of the existence of sandwiches we have the following very general theorem:

THEOREM: Given n Lebesgue Measurable sets A_1, A_2, ... , A_n all subsets of R^n there exists an n-1 plane that simultaneously divides the Jordan Measure of each set.

Unfortunately we lack the bravado to tackle such a broad theorem and settled on tackling a more manageable special case concerning polygons.

OUR SPECIAL CASE THEOREM:

"Given two polygons in R^2 there exists a line that simultaneously bisects the area of both polygons."

The theorem isn't constructive, so although such a line is guaranteed to exist, there is no procedure given to make the perfect cut, which we will dub the 'Hamline'.

THAT'S WHERE WE AND OUR GOOD PAL SAGE COME IN!!!

So, we broke our work up into two teams: P-Team and L-Team

P-Team Objectives:

Write an interact that plots polygons
Calculate Polygon Area
Generate the new polygons found by intersecting a line and calculate their area

L-Team Objective:

Given the P-Team code write an intelligent line generating algorithm that generates a reasonable hamline.

Here is the tender fruit of P-Team's labor.

Basically this interact lets you plot two polygons. Next it calculates the area and centroids of our two polygons. Finally, we draw a line that connects the centroids and calculate the upper and lower areas formed.

Notice that in addition to these cool plots we can:

add points
change our colors

The Lights Dim... The Room Goes Silent... All Eyes turn to the Front of the Room...

L-TEAM BREAKDOWN ENSUES

Classes are better than not having classes
Originally thought we would use rotation and translate but found a better way... rectangles.
Talk about rectangle method.

Similar to how the rain and cloud become one...

L-Team and P-Team united into one furious interact.

Don't let it's beauty distract you. This is very practical code. Notice the new things present:

We now have an error slider. This tells our line search to terminate once it finds a line that gives an area bisection with the property that | UP1/BP1 - .5 | < error.
We have a new line dubbed hamline that is found by the above method. Notice how foolish and lazy our old centroid line looks.

Some Pretty Cool Polygons that are just Begging to be bisected:

If everything went smoothly that's great, but it turns out that our method is sometimes too demanding and our line generating technique cannot find a line that meets our stringent requirements. If this is the case we wrote code that relaxes our requirements and works much more robustly.

CONCLUSION

We learned that glueing code together is way tough and that class oriented coding is great for efficiency.

We just barely scratched the surface of neat things to do with this theorem. A couple future endeavors may include:

Integrating our Polygon Class stuff (i.e Polygon.Area and Polygon.Centroid) into Sage
A Polygon Perimeter Bisection Program
General shapes in the Plane Bisection
The 3D volume bisection analogue
Turkey Sandwiches?

Select a color

Select a color

Click to the left again to hide and once more to show the dynamic interactive window

Select a color

Select a color

Click to the left again to hide and once more to show the dynamic interactive window

Select a color

Select a color

Click to the left again to hide and once more to show the dynamic interactive window

Select a color

Select a color

Click to the left again to hide and once more to show the dynamic interactive window

Select a color

Select a color

Click to the left again to hide and once more to show the dynamic interactive window