import numpy as np import random as R import math ############################## ############################## ############################## class GridException(Exception): pass class GridEnvironment: def __init__(self, size, cyclic=False): ### Initialize a square 2D grid of side length from 0 to ### 'size'. A grid position is given as the tuple (x,y) where ### x and y are integers. ### ### The grid needs to keep track of the status of each grid cell: ### whether it's empty or occupied. ### ### If cyclic is True, the edges of the grid should wrap ### around, i.e. any given position (x,y) should be accepted ### and mapped back to the initial cell. If cyclic is False, ### the boundary of the grid should be permanently marked as ### occupied. self.size = size self.cyclic = cyclic ### extend this def fill(self,pos): ### Try to fill position pos ( where pos is a tuple (x,y) ). ### If pos is full already, return False ### If pos is free, mark it as full and return True ### fix this return True def erase(self,pos): ### Try to clear grid position pos. No return value. ### If it is already empty, raise a GridException. ### fix this pass def randpos(self): ### Return a randomly chosen _free_ grid cell as an (x,y) tuple. ### fix this return (12,17) def impose_limit(self,pos): ### When cyclic is True, map 'pos' back to the main grid square. ### The input value 'pos' can lie anywhere, the output value ### (x,y) should lie in the main grid. ### fix this return pos def filled(self): ### Return the list of xs and ys for all filled positions. xs = [] ys = [] ### fix this return (xs+[1,2,3,4,6,8,12,15,19,21,22,23,24, 1,6,9,11,15,16,18,19,21, 1,2,3,6,10,15,17,19,21,22,23, 1,6,9,11,15,19,21, 1,6,8,12,15,19,21,22,23,24], ys+[5,5,5,5,5,5,5,5,5,5,5,5,5, 4,4,4,4,4,4,4,4,4, 3,3,3,3,3,3,3,3,3,3,3, 2,2,2,2,2,2,2, 1,1,1,1,1,1,1,1,1,1]) ############################## ############################## ############################## class AtomException(Exception): pass class Atom: def __init__(self,N,flexible=True): ### Initialize an Atom species with N arms (0 <= N <= 4). ### They can attach to another atom if both atoms have an arm ### free that can point in the right direction. ### ### For 2-armed atoms, 'flexible' determines if the arms can ### bend. If flexible is True, the two arms can point in any ### direction. If flexible is False, they must point opposite ### to each other. In that case the first successful ### attachment determines the direction. ### ### 'flexible' is not used for any other N != 2 ### fix this pass def canAttach(self,d): ### return True if there is a free arm that can point to direction 'd' ### return False otherwise. Do not change the status of the arms yet! # d is the direction from Atom a to b, and will be given as # 0 # ^ # 1 < a > 3 # v # 2 ### fix this return True def doAttach(self,d): ### Mark an arm as attached in direction 'd'. It is not ### necessary to repeat the 'canAttach' check here. ### fix this pass ############################## ############################## ############################## class Cluster: ### ### A cluster is a rigid collection of attached atoms. It should ### keep track of all its component atoms and their positions ### def __init__(self, env, affinity=1.0, Narms=[4], rigidPortion=0.8): ### We always initialize a cluster with a single atom. You ### can keep the following lines that create the initial atom ################################################## self.env = env self.affinity = affinity ### get a random position from env pos = env.randpos() ### choose a number of arms for the atom from the list Narms N = R.choice(Narms) a = Atom(N) ### If Narms = 2, check if they should be flexible arms if N == 2 and R.random() < rigidPortion: a = Atom(N, flexible=False) ################################################### ### extend this... ### to keep track of the atoms and their positions ### the first atom to be stored is 'a' with position 'pos' ### tell the grid environment that this cell is now busy self.env.fill(pos) def randStep(self): ### return a single step, expressed as a (dx,dy) tuple ### fix this return (0,0) def step(self): ### Move all member atoms one step in the direction chosen by ### self.randStep(). If any of them encounters an obstacle (hint: ### ask the grid about them), leave the whole cluster as it ### was before. Do _not_ move parts of the cluster and leave ### parts behind. ### ### Remember to tell the grid environment to update its ### occupied positions. ### fix this pass def pos(self): ### Return an np.array of the list of the cluster's member atoms' positions ### fix this positionlist = [(20,35),(21,35),(22,35),(23,35),(25,35), (27,35),(31,35),(34,35),(38,35),(40,35), (41,35),(42,35),(43,35),(20,34),(25,34), (28,34),(30,34),(34,34),(35,34),(37,34), (38,34),(40,34),(20,33),(21,33),(22,33), (25,33),(29,33),(34,33),(36,33),(38,33), (40,33),(41,33),(42,33),(20,32),(25,32), (28,32),(30,32),(34,32),(38,32),(40,32), (20,31),(25,31),(27,31),(31,31),(34,31), (38,31),(40,31),(41,31),(42,31),(43,31)] return np.array( positionlist ) def touches(self,pos): ### Return a list of cluster positions and atoms that touch position 'pos'. ### Format: [ (pos1,a1), (pos2,a2), (pos3,a3) ] ### ### 'touch' means atoms directly next to 'pos' along the x and y ### axes. Diagonally placed atoms do _not_ touch. ### ### Also, only a proportion of the potential touches are actual touches, ### determined by 'self.affinity'. In the extreme cases: ### If it is 1.0 all touches count. ### If it is 0.0 no touch counts. ### ### fix this return [] def merge(self,c2): ### Merge this cluster with cluster 'c2'. 'self' should take ### over c2's atoms, and c2 should be emptied. ### ### fix this pass class ClusterGroup: def __init__(self, N, env, style='o', **kwargs): ### ### Create a cluster group with N one-atom clusters. ### ### The '**kwargs' construct is used to pass all flags that we ### don't need here on to the Cluster __init__ function ### # list of member clusters self.clusters = [ Cluster(env,**kwargs) for i in range(N) ] # plotting style self.style = style def step(self): ### let all clusters in this group do one step for c in self.clusters: c.step() ### Now merge all clusters that are touching. ### Use the touches() from Cluster! ### Clusters can only merge if the two neighbouring atoms both ### still have free arms. Use the canAttach() function from ### 'Atom' for this. ### Once established that a merge is possible, mark the atoms ### using their doAttach() function, merge the clusters using merge(), ### then remove the empty cluster from this cluster group ### extend this def pos(self): poslist = [ c.pos() for c in self.clusters ] return np.concatenate(poslist) def xs(self): return self.pos()[:,0] def ys(self): return self.pos()[:,1] def fractaldim(xs,ys): ### Calculate the fractal dimension of the arrangement of points ### given by xs and ys. The procedure should follow the one shown ### in the lecture, starting from the centre point of the bounding ### box around the xs and ys. ### fix this return -1.23456 if __name__ == "__main__": import pylab as pl from Read_Type import read_float env = GridEnvironment(51, cyclic=True) # env = GridEnvironment(51, cyclic=False) Ninit = int(read_float('Number of initial one-atom clusters: ', minimum=2,maximum=300, default=100)) steps = int(read_float('Number of time steps: ', minimum=1,maximum=2000, default=100)) groups = [ ClusterGroup(Ninit,env,'b.',affinity=1.0), #ClusterGroup(50,env,'r.',affinity=0.02), ] # use something like this to test for rigid 2-armed atoms # groups = [ ClusterGroup(Ninit,env,'b.',affinity=1.0,Narms=[2],rigidPortion=1.0) ] # use something like this to test for a 3-to-1 mix of 2- and 3-armed atoms # groups = [ ClusterGroup(Ninit,env,'b.',affinity=1.0,Narms=[2,2,2,3],rigidPortion=1.0) ] pl.ion() scatter = pl.subplot('111') plots = {} for group in groups: plots[group] = scatter.plot(group.xs(), group.ys(), group.style)[0] scatter.set_aspect('equal') xs,ys = env.filled() envplot = scatter.plot(xs,ys,'ro',alpha=0.2)[0] scatter.set_xlim(0,env.size-1) scatter.set_ylim(0,env.size-1) update = 1 for i in range(steps): for group in groups: group.step() if (i+1) % update == 0: plots[group].set_data(group.xs(), group.ys()) if (i+1) % update == 0: oldxs, oldys = xs, ys xs,ys = env.filled() envplot.set_data(xs,ys) scatter.set_title( '%s of %s steps' % (i+1,steps) ) pl.draw() if any([len(group.clusters) == 1 for group in groups]): break pl.draw() pl.ioff() if env.cyclic == True: print "Fractal dimension only printed for non-cyclical boundaries" else: for group in groups: print 'Fractal dimension is %s' % fractaldim(group.xs(),group.ys()) pl.show()