cos = math.cos
sin = math.sin
pi = math.pi
class diskObj:
"""
A disk object for use in a KWPolymer object.
Has a label, a radius, a center,
and a dictionary of the form
{adjacent:angleToAdjacent for each adjacent disk}.
"""
def __init__(self,label,radius,center=None):
"""
By default, center is None.
"""
self.label = label
self.radius = radius
self.center = center
def __str__(self):
"""
Print disk label.
"""
return str(self.label)
def branchbranch(x0,y0,a0,b0,r0,x1,y1,a1,b1,r1):
a = a0^2 + b0^2 - 2*a0*a1 + a1^2 - 2*b0*b1 + b1^2
b = 2*x0*a0 + 2*y0*b0 - 2*a0*x1 - 2*b0*y1 - 2*x0*a1 + 2*x1*a1 - 2*y0*b1 + 2*y1*b1
c = x0^2 + y0^2 - r0^2 - 2*x0*x1 + x1^2 - 2*y0*y1 + y1^2 - 2*r0*r1 - r1^2
if a == 0:
return -b/c
disc = b^2-4*a*c
if disc < 0:
return -1
r0 = (-b+sqrt(disc))/(2*a)
r1 = (-b-sqrt(disc))/(2*a)
if r1 > r0 > 0:
return r0
if r0 > r1 > 0:
return r1
return -1
class KWPolymer:
"""
A branched polymer class implementing the
Kenyon-Winkler algorithm for generating
random branched polymers in the plane.
"""
def __init__(self, radii = None, printOn=False):
"""
Constructs a new polymer consisting of a disk
at the origin. Radii must be iterable.
"""
if radii is None:
def unit():
while True:
yield 1
radii = unit()
self.nextRadius = radii.__iter__().next
self.G = DiGraph()
self.root = diskObj(label = 1,
radius = self.nextRadius(),
center = (0,0))
self.G.add_vertex(self.root)
self.printOn = printOn
def grow(self):
"""
Add a new disk to the polymer using the Kenyon-Winkler
algorithm. The algorithm works inductively, providing
a recipe for generating a uniformly random polymer of
order N+1 from one of order N. Given that an order N
polymer has already been generated, our implementation
of the algorithm consists of the following steps:
1. Place a new disk of radius zero at a uniformly
chosen point along the boundary of the polymer.
2. Find the smallest radius of the new disk less
than one at which a cycle forms in the polymer.
3. If there is no such radius, increase the new
disk's radius to 1, recalculate the centers of
the other disks while holding the new disk fixed,
and skip step 4. Otherwise, increase the new
disk's radius until a cycle forms, recalculate
the centers of the other disks while holding the
new disk fixed, and continue to step 4.
4. If L_i is the local volume change near the
polymer in the configuration space of the tree
formed by removing the i-th edge of the cycle,
with respect to a small increase in the radius of
the new disk, trees with positive L_i are those
that allow us to continue increasing the radius
of the new disk without causing an overlap.
Choose among these trees, weighting them by the
L_i. Remove the corresponding edge from the
polymer. Return to step 2.
"""
## Step 1 ##
disks = self.G.vertices()
parent = disks[ randint(0,len(self.G)-1) ]
angle_from_parent_to_new = random() * 2 * pi
targetRadius = self.nextRadius()
new = diskObj(label = len(self.G) + 1, radius = 0)
self.G.add_vertex(new)
self.addEdge(parent, new, angle_from_parent_to_new)
self.getCenterWithRespectTo(new, parent)
while True:
branches = [
(
neighbor,
self.getBranch(neighbor,new),
self.G.edge_label(new,neighbor)
)
for neighbor in self.G.neighbors(new)
]
currentRadius = new.radius
## Step 2 ##
collRadius, A, B = self.detectCollision(targetRadius, new, branches)
if self.printOn: print collRadius, A, B
## Step 3 ##
if A is None:
self.moveDisks(targetRadius - currentRadius, new, branches)
new.radius = targetRadius
if self.printOn:
plot = self.plotPolymer(True, True, True)
plot.show(axes=False, aspect_ratio=1)
break
else:
self.moveDisks(collRadius - currentRadius, new, branches)
new.radius = collRadius
## Step 4 ##
self.breakCycle(A, B, new)
if self.printOn:
plot = self.plotPolymer(True, True, True)
plot.show(axes=False, aspect_ratio=1)
def getBranch(self, inside, outside):
"""
Precondition: inside is adjacent
to outside. Returns the set of all
disks connected to inside except
through outside, i.e. the maximal
subtree containing inside but not
outside.
"""
angle = self.G.edge_label(inside, outside)
self.deleteEdge(inside, outside)
branch = self.G.connected_component_containing_vertex(inside)
self.addEdge(inside, outside, angle)
return branch
def detectCollision(self, targetRadius, new, branches):
"""
If new's radius can increase to targetRadius
without any collisions occurring between disks
(i.e., without any cycles forming), return
(0,None,None). Otherwise, return (collRadius,A,B).
collRadius is the smallest radius greater than
new's current radius but less than targetRadius
such that a collision occurs. A and B are the
disks that collide.
Let a branch be a maximal subtree of the tree
representing the polymer that does not contain new.
There are two cases to consider when checking for
collisions:
1. A disk in one branch collides with new.
2. A disk in one branch collides with a disk
in another branch.
Note that disks never collide within a branch,
since new pushes the entire branch in the same
direction. Also, new can only have as many
branches as the number of disks that can be
packed around it.
"""
collRadius, A, B = targetRadius, None, None
oldRadius = new.radius
## Case 1 ##
x0, y0 = new.center
for (n, branch, angle) in branches:
for disk in branch:
if disk is n:
continue
x1, y1 = disk.center
dx = x1 - x0
dy = y1 - y0
p = disk.radius + oldRadius
bottom = 2 * (dx*cos(angle) + dy*sin(angle) - p)
if bottom != 0:
r = (p^2 - dx^2 - dy^2) / bottom
if 0 < r < collRadius - oldRadius:
collRadius = r + oldRadius
A = disk
if self.printOn: print A, '\n', collRadius, '\n', '\n', x1, '\n', y1, '\n', x0, '\n', y0, '\n'
if A is not None:
B = new
## Case 2 ##
branchPairs = Combinations(branches, 2)
for [(n0, branch0, angle0), (n1, branch1, angle1)] in branchPairs:
for disk0 in branch0:
x0, y0 = disk0.center
r0 = disk0.radius
cos0 = cos(angle0)
sin0 = sin(angle0)
for disk1 in branch1:
x1, y1 = disk1.center
r1 = disk1.radius
cos1 = cos(angle1)
sin1 = sin(angle1)
dist = disk0.radius + disk1.radius
a = cos1 - cos0
b = sin1 - sin0
c = x1 - x0
d = y1 - y0
bottom = a^2 + b^2
r = branchbranch(x0, y0, cos0, sin0, r0, x1, y1, cos1, sin1, r1)
if 0 < r < collRadius - oldRadius:
collRadius = r + oldRadius
A, B = disk0, disk1
if self.printOn:
print A
print B
print collRadius
print dist
print
print 'a = ', a
print 'b = ', b
print 'c = ', c
print 'd = ', d
print 'discriminant = ', discriminant
return collRadius, A, B
# if bottom != 0:
# discriminant = (a*c+b*d)^2-(a^2*b^2)*(c^2+d^2-dist^2)
# if discriminant >= 0:
# r = -(a*c+b*d + sqrt(discriminant)) / bottom
# if 0 < r < collRadius - oldRadius:
# collRadius = r + oldRadius
# A, B = disk0, disk1
# if self.printOn:
# print A
# print B
# print collRadius
# print dist
# print
# print 'a = ', a
# print 'b = ', b
# print 'c = ', c
# print 'd = ', d
# print 'discriminant = ', discriminant
# return collRadius, A, B
def moveDisks(self, changeInRadius, new, branches):
"""
Keep the root branch fixed. Move new
and its other branches appropriately.
"""
parent = self.G.shortest_path(new, self.root)[1]
parentAngle = self.G.edge_label(parent, new)
dx0 = changeInRadius * cos(parentAngle)
dy0 = changeInRadius * sin(parentAngle)
x0, y0 = new.center
new.center = (x0 + dx0, y0 + dy0)
for (neighbor, branch, angle) in branches:
if neighbor is not parent:
dx1 = changeInRadius * cos(angle)
dy1 = changeInRadius * sin(angle)
for disk in branch:
x1, y1 = disk.center
disk.center = (x1 + dx0 + dx1, y1 + dy0 + dy1)
def breakCycle(self, A, B, new):
"""
Remove one of the edges in the cycle
that formed when A and B touched and
add an edge between A and B. (Since
A and B touched, the edge between
them is conveniently never the one
that is removed.) Choose among the
edges for which the local volume
change is positive as the radius of
new increases, weighting each
according to the magnitude of the
volume change.
"""
new_to_A = self.G.shortest_path(new, A)
B_to_new = self.G.shortest_path(B, new)
B_to_new.pop()
cycle = new_to_A + B_to_new
averageAngle = (self.getAngle(cycle[-1],cycle[0]) + \
self.getAngle(cycle[0],cycle[1])) / 2
L = len(cycle)
dp0 = self.dot(cycle[0], cycle[1], averageAngle)
dotproducts = [
(
(cycle[i], cycle[(i+1)%L]),
dp0 * self.dot(cycle[i], cycle[(i+1)%L], averageAngle)
)
for i in range(L)
]
goodTrees = [(edge, -dp) for (edge, dp) in dotproducts if dp < 0]
choice = self.chooseAmong(goodTrees)
self.deleteEdge(choice[0], choice[1])
self.addEdge(A, B, self.getAngle(A, B))
def dot(self, A, B, angle):
"""
Returns the dot product
of the vector from A to
B with the unit vector in
the direction angle.
"""
a,b = A.center
c,d = B.center
return (c-a)*cos(angle) + (d-b)*sin(angle)
def chooseAmong(self, choices):
"""
Given a list a tuples of the form
(choice, weight), weight each choice
with probability weight/sum(weights),
sample a choice, and return it.
Weights must be positive.
"""
sumOfWeights = sum(weight for (choice,weight) in choices)
normalizedChoices = [(choice, weight/sumOfWeights)
for (choice,weight) in choices]
rand = random()
for (choice, probability) in normalizedChoices:
if rand > probability:
rand -= probability
else:
return choice
def getCenterWithRespectTo(self, A, B):
"""
Precondition: the two disks are adjacent.
Compute disk A's center given disk B's
center and the angle of the edge from
disk B to disk A.
"""
x, y = B.center
magnitude = A.radius + B.radius
angle = self.G.edge_label(B,A)
dx = magnitude * cos(angle)
dy = magnitude * sin(angle)
A.center = (x+dx, y+dy)
def addEdge(self, A, B, angle_from_A_to_B):
"""
Adds an edge between A and B to
the graph of the polymer, weighted
with the given angle in the
direction from A to B.
"""
self.G.add_edge(A, B, angle_from_A_to_B)
self.G.add_edge(B, A, angle_from_A_to_B + pi)
def deleteEdge(self, A, B):
"""
Precondition: there exists an edge
between A and B. Removes this edge.
"""
self.G.delete_edge(A, B)
self.G.delete_edge(B, A)
def getAngle(self, A, B):
"""
Returns the angle of the vector from
the center of disk A to the center of
disk B.
"""
x = B.center[0] - A.center[0]
y = B.center[1] - A.center[1]
return CC(x,y).arg()
def plotPolymer(self, disksOn=False, labelsOn=False, linesOn=False):
"""
Return a plot of the polymer.
"""
g = Graphics()
for disk in self.G.vertices():
if disksOn:
g += circle(disk.center, disk.radius, rgbcolor='blue')
if labelsOn:
g += text(disk.label, disk.center, rgbcolor='blue', vertical_alignment='top', horizontal_alignment='left')
if linesOn:
for edge in self.G.edges():
g += line((edge[0].center, edge[1].center), rgbcolor='gray')
return g
def disks(self):
return self.G.vertices()
def geometricRadius(self):
"""
Maximal distance from center of mass realized
in the polymer.
"""
X,Y = self.centerOfMass()
D = []
for disk in self.disks():
x,y = disk.center
d = sqrt((x-X)^2 + (y-Y)^2)
D.append(d + disk.radius)
return max(D)
def geometricDiameter(self):
"""
Maximal distance between pairs of points in
the polymer.
"""
m = 0
for disk0,disk1 in Combinations(self.disks(),2):
x0,y0 = disk0.center
x1,y1 = disk1.center
d = disk0.radius + disk1.radius + sqrt((x1-x0)^2 + (y1-y0)^2)
if m < d:
m = d
return m
def totalArea(self):
"""
Total area = sum of area of disks.
"""
return sum(pi*(d.radius^2) for d in self.disks())
def centerOfMass(self):
"""
Center of mass = weighted average of centers of disks,
X = 1/A sum( pi * x_i * r_i^2 : i = 1 .. n )
Y = 1/A sum( pi * x_i * r_i^2 : i = 1 .. n )
"""
A = self.totalArea()
X = sum(pi*(d.center[0])*(d.radius^2) for d in self.disks()) / A
Y = sum(pi*(d.center[1])*(d.radius^2) for d in self.disks()) / A
return (X,Y)
