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:
L-Team Objective:
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.
{{{id=1| %hide def PArea(L): S=[] for i in range(0,len(L)): for z in L[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) 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) def cent_x(L): """Takes a list of Pgon vertices and finds the x-centroid""" S=[] for i in range(0,len(L)): for z in L[i]: S.append(z) B=PArea(L) 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(L): """Takes a list of PGon vertices and finds the y-centroid""" S=[] for i in range(0,len(L)): for z in L[i]: S.append(z) B=PArea(L) 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 centroid(L): M=[] M.append(cent_x(L)) M.append(cent_y(L)) return tuple(M) def draw_line(pt1, pt2): """Returns a line as a list [a,b,c] for ax+by=c""" if pt2[0]-pt1[0] != 0: m = (pt2[1] - pt1[1])/(pt2[0]- pt1[0]) b = pt1[1] -m*(pt1[0]) return [-m,1 ,b] else: m = (pt2[0] - pt1[0])/(pt2[1]- pt1[1]) b = pt1[0] - m*(pt1[1]) return [1,-m,b] def does_intersect(pt1,pt2,line): line1=draw_line(pt1,pt2) line2=line 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) def new_pgon(L,line): Q=[] for i in range(len(L)): if i < len(L)-1: r=does_intersect(L[i],L[i+1],line) if r != None: Q.append(r) g=L[i+1] if g[1] < -(line[0]/line[1])*g[0]+(line[2]/line[1]): Q.append(g) if i == len(L)-1: r= does_intersect(L[i],L[0],line) if r!= None: Q.append(r) g=L[0] if g[1] < -(line[0]/line[1])*g[0]+(line[2]/line[1]): Q.append(g) return Q def how_big(x,y): if xNotice that in addition to these cool plots we can:
The Lights Dim... The Room Goes Silent... All Eyes turn to the Front of the Room...
L-TEAM BREAKDOWN ENSUES
Similar to how the rain and cloud become one...
L-Team and P-Team united into one furious interact.
{{{id=3| %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 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 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("(1,-2) (2,2) (4,-1) (6,5) (-1,2)", "Vertices P1 (separated by spaces)", type=str), clr=color_selector(Color('purple'), widget='jpicker'), vertices2=input_box("(-1,-1) (-2,-2) (-5,-1/2) (-5,-8) (-4,-3) (-3,-7) (-5/2,-2) (-2,-6) (0,-1)", "Vertices P2 (separated by spaces)", type=str), clr2=color_selector(Color('red'), widget='jpicker'), error=slider(0,.5,default= 0.0360721442885772 )): ### 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)) /// }}}Don't let it's beauty distract you. This is very practical code. Notice the new things present:
Some Pretty Cool Polygons that are just Begging to be bisected:
{{{id=8| %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 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 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("(2,-5) (3, -5) (3,6) (5,1) (7,6) (7,-5) (8,-5) (8,8) (7,8) (5,3) (3,8) (2,8)", "Vertices P1 (separated by spaces)", type=str), clr=color_selector(Color('purple'), widget='jpicker'), vertices2=input_box("(-10,-5) (-9,-5) (-9,1) (-6,1) (-6,-5) (-4,-5) (-3,1) (-1,1) (-7/6,2) (-17/6,2) (-2,7) (0,-5) (1,-5) (-1,8) (-3,8) (-5,-5) (-5,8) (-6,8) (-6,2) (-9,2) (-9,8) (-10,8)", "Vertices P2 (separated by spaces)", type=str), clr2=color_selector(Color('red'), widget='jpicker'), error=slider(0,.5,default=.0400801603206413)): ### 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( 'BP1/P1= %s'%((BA1)/A1)) show('BP2/P2 = %s' %((BA2)/A2)) /// }}} {{{id=11| %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 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 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("(15,10) (60,10) (60,45) (65,50) (65,55) (60,60) (45,60) (35,55) (30,60) (15,60) (10,55) (10,50) (15,45)", "Vertices P1 (separated by spaces)", type=str), clr=color_selector(Color('purple'), widget='jpicker'), vertices2=input_box("(-15,-10) (-60,-10) (-60,-45) (-65,-50) (-65,-55) (-60,-60) (-45,-60) (-35,-55) (-30,-60) (-15,-60) (-10,-55) (-10,-50) (-15,-45)", "Vertices P2 (separated by spaces)", type=str), clr2=color_selector(Color('red'), widget='jpicker'), error=slider(0,.5,default=.45)): ### 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( 'BP1/P1= %s'%((BA1)/A1)) show('BP2/P2 = %s' %((BA2)/A2)) /// }}}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.
{{{id=12| %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)) /// }}}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:
