Percolation and cluster growth
Remarks on completing the module
This assignement is summatively assessed. It is imperative that you submit the notebook on time.
In [15]:
Cluster growth
The cell below defines the Cluster class. Complete the methods that are only
partially implemented.
You can the use the following cells to test your implementation.
Completing the class Cluster is the main part of the assignment. Some parts of the methods need to be
completed. The cells below it take you through this, one method at a time. There are also some cells to verify
your implementation.
The main components of an instance of the class are:
ndim: the dimension of the cluster. We only examine the special case of ndim=2
periodic: if True, impose periodic boundary conditions, else, impose isolated boundary conditions
size: the extent of a cluster in any dimension
grid: a dictionary of positions currently occupied by the cluster, in the form of tuples for example ((0,1):
True) if element (0,1) is part of the cluster
boundary: a dictionary of nearest-neighbours of cluster elements, that are not occupied
offsets: a list of all neigbours of location (0,0)
for ndim=2, this is simple [(1,0), (-1,0), (0,1), (0,-1)]
import numpy as np
from scipy import special
import random as R
import math
import pylab as pl
from matplotlib import animation,rc
import matplotlib.pyplot as plt
rc('animation', html='html5')
from ipywidgets import IntProgress
from IPython.display import display
2019/11/7 Percolation
https://notebooks.dmaitre.phyip3.dur.ac.uk/miscada-sc/user/czjs88/notebooks/Percolation 19/Percolation.ipynb# 2/21
In [ ]:
class Cluster:
# initialize the cluster
def __init__(self, ndim=2, size=32, random_number_seed=-1, periodic=False):
self.size = size
self.ndim = ndim
self.periodic = periodic
self.grid = {}
self.boundary = None
# Initialize random number generator if random_number_seed is > 0
if random_number_seed > 0:
np.random.seed(random_number_seed)
else:
np.random.seed(None)
# create list of all nearest neighbours of element (0,) * ndim
self.offsets = []
for dim in range(self.ndim):
offset = np.zeros(ndim, dtype=int)
offset[dim] = 1
self.offsets.append(tuple(offset))
self.offsets.append(tuple(-offset))
self.noffsets = len(self.offsets)
# The function return True / False depending on whether position is inside / ou
def element_inside(self, position):
if self.periodic:
return True
sp = np.array(position)
return (sp >= 0).all() and (sp < self.size).all()
# This function adds the element at location "position" to the cluster.
# It returns an exception if that element is already in the cluster, or if it
def element_add(self, position):
# This routine adds the element at location "position" to the cluster
# if the element is already part of the cluster, raise an exception
# if periodic=True: wrap position periodically before inserting element
# if periodc=False: raise exception if element is outside the grid bounda
pos = np.array(position)
if self.periodic:
pos = (pos + 2 * self.size) % self.size
if(self.element_inside(pos)):
if tuple(pos) in self.grid:
print(" ++ element_add, Position = ", pos)
raise Exception("Element is already in cluster")
else:
self.grid[tuple(pos)] = True
else:
print(" ++ element_add, Position = ", pos)
raise Exception("Element is outside of grid")
# The function returns the list of all nearest-neighbours of a given position, t
def element_neighbours(self, position):
# consistency check:
if len(position) != self.ndim:
raise Exception(" Error: this position does not have the right dimension
2019/11/7 Percolation
https://notebooks.dmaitre.phyip3.dur.ac.uk/miscada-sc/user/czjs88/notebooks/Percolation 19/Percolation.ipynb# 3/21
# use self.offsets to return a list of all nearest-neighbours of the specifi
neighbours = []
# YOUR CODE HERE
raise NotImplementedError()
return neighbours
# This functon returns the distance between position1 and position2
def distance(self, position1, position2):
# consistency check
if len(position1) != self.ndim:
raise Exception("not a valid position 1 in distance")
if len(position2) != self.ndim:
raise Exception("not a valid position 2 in distance")
#
p = np.asarray(position1)
q = np.asarray(position2)
distance = 0
for i in range(len(p)):
dx = np.abs(p[i]-q[i])
if self.periodic:
if dx > self.size/2:
dx = self.size - dx
distance += dx * dx
return np.sqrt(distance)
# This function returns all the distances from each element of the clyuster, to
def element_distances(self, position):
distance = {}
for element in self.grid:
distance[element] = self.distance(element, position)
return distance
# This function periodically wraps "position" on the grid
def periodic_wrap(self, position):
if self.periodic:
return tuple(np.array(position) % self.size)
else:
return position
# Eden model for cluster growth
# The algorithm for adding an element is as follows
# - pick at random, any unoccupied site that is a nearest neighbour of an eleme
# - add it to the cluster
# The method below should add 1 element to the cluster based on this algorithm
def Eden(self):
# The first time Eden in invoked, we create a list, called self.boundary,
# of empty sites that neighbour the cluster
# When Eden is called again, we update self.boundary
if self.boundary == None:
self.boundary = {}
# iterate over all neighbours of all elements in the cluster
# if the neighbour is not in the cluster, and not in self.boundary, a
2019/11/7 Percolation
https://notebooks.dmaitre.phyip3.dur.ac.uk/miscada-sc/user/czjs88/notebooks/Percolation 19/Percolation.ipynb# 4/21
# YOUR CODE HERE
raise NotImplementedError()
# pick at random an element of self.boundary, and add this element to the cl
# YOUR CODE HERE
raise NotImplementedError()
# remove this element from self.boundary
# YOUR CODE HERE
raise NotImplementedError()
# add the neighbours of the element to self.boundary,
# provided they are not in the cluster, and not already in the boundary
for neighbour in self.element_neighbours(element):
if neighbour not in self.grid:
if neighbour not in self.boundary:
self.boundary[self.periodic_wrap(neighbour)] = True
# return current number of elements in the boundary
return len(self.boundary)
# DLA model for generating a cluster
# Start a random walk at distance "radius" away from the seed of the cluster, wh
# Terminate the walk, if either
# walker hits the cluster - then add it to the cluster
# crosses the edge, i.e. reaches a distances > edge from the centre
# DLA adds an element to cluster, and the function returns "track": the list of
# random walk
def DLA(self, centre, radius, edge):
start = np.asarray(centre)
# choose starting point in polar coordinates - then convert to cartesian coo
rad = np.asarray([1, int(radius)])
rad = rad.max()
phi = 2 * np.pi * np.random.random() # angle uniform in [0, 2pi[
x = rad * np.cos(phi)
y = rad * np.sin(phi)
dr = np.sqrt(x*x+y*y)
# round-off to integer coordinates
start[0] += int(x)
start[1] += int(y)
start = tuple(start)
cont = True
# track stores all the location that the random walk visits
track = []
track.append(start)
step = start # the random walk steps from the current step, "step", to
while cont:
# continue walking as long as cont = True
# take a random step, from step to next_step
# Hint: use RandomStep, implemented below, or write your own
# YOUR CODE HERE
raise NotImplementedError()
# If any neighbour of next_step is in the cluster, then:
# - add next_step to cluster
# - add next_step to track
2019/11/7 Percolation
https://notebooks.dmaitre.phyip3.dur.ac.uk/miscada-sc/user/czjs88/notebooks/Percolation 19/Percolation.ipynb# 5/21
# - set cont = False, to end walk
# YOUR CODE HERE
raise NotImplementedError()
# if cont is True, we need to check if walker is still in the allowed re
# if distance to centre > edge: put cont = False, add element to tra
# if distance to centre < edge: cont = True
# - add next_step to track
# - put step = next_step
# - continue walk
# YOUR CODE HERE
raise NotImplementedError()
# walk finished: return the track
# You may want to uncomment the lines below for debugging your code
# if attach:
# print("Walker attached to cluster at position = ", next_step)
# else:
# print("walker off grid")
return track
# This function randomly adds sites to a grid, based on comparing a uniform rand
# to the input (given) probability "probability"
def Random(self, probability):
# for each site: if random number < probability, add element to cluster
shape = (self.size,) * self.ndim
prob = np.random.rand(*shape)
assign = np.where(prob < probability)
assign = np.array(assign).T
for a in assign:
self.element_add(a)
# Perform a single random step, starting at location "position"
# function returns the location of the end of the step
def RandomStep(self, position):
itry = np.random.choice(range(self.noffsets))
offset = self.offsets[itry]
new = tuple(np.array(position) + np.array(offset))
return self.periodic_wrap(new)
# Starting from location "position", take nstep random walk steps
def RandomWalk(self, position, nstep):
track = []
current = position
for i in range(nstep):
current = self.RandomStep(current)
track.append(current)
return track
# Plot the current cluster
# Input is the cut to plot, and, optionally, the name of the plot, and the nam
# to save the plot
def ShowCluster(self, cut, title=None, file=None):
# extract the cooridnates of the cluster specified by cut
xs, ys = self.ExtractSlice(cut)
# create the figure, and plot all cluster sites
2019/11/7 Percolation
https://notebooks.dmaitre.phyip3.dur.ac.uk/miscada-sc/user/czjs88/notebooks/Percolation 19/Percolation.ipynb# 6/21
fig, ax = pl.subplots(figsize=(11,11))
ax.set_aspect('equal')
ax.set_xlim(-1, self.size)
ax.set_ylim(-1, self.size)
ax.set_xlabel('x')
ax.set_ylabel('y')
if not title == None:
ax.set_title(title)
ax.grid()
ax.plot(xs,ys,'ro',alpha=0.2)[0]
# save the figure
if not file == None:
print("Saving cluster to file: ", file)
plt.savefig(file)
plt.show()
# Extract a slice: this is useful if the cluster is in ndim>2 dimensions
# cut specifies the direction to cut the cluster and extract a 2D slice
def ExtractSlice(self, cut):
# cut is a tuple of length self.ndim
# the tuple should have exactly 2 elements equal to -1, which define the pla
if len(cut) != self.ndim:
raise Exception(" cut has wrong length")
plane = np.where(np.asarray(cut) == -1)[0]
if len(plane) != 2:
raise Exception(" plane has wrong dimenson ", len(plane))
for dim in range(self.ndim):
if dim not in plane:
if cut[dim] < 0 or cut[dim] >= self.size:
raise Exception(" cut is out of range ", dim, cut[dim])
#
sites = []
xs = []
ys = []
cut = np.asarray(cut)
for site in self.grid:
if np.all(np.logical_or(cut==-1,cut==np.asarray(site))):
sites.append(site)
xs.append(site[plane[0]])
ys.append(site[plane[1]])
return xs, ys
# This function plots a track, as generated by DLA for example, for a given rand
# Input is the track, and optionally, the name of the file to save the figure
def ShowTrack(self, track, file=None):
fig, ax = pl.subplots(figsize=(7,7))
ax.set_aspect('equal')
ax.set_xlim(-1, self.size)
ax.set_ylim(-1, self.size)
ax.set_xlabel("x-position")
ax.set_ylabel("y-position")
ax.set_title("Track of random")
ax.grid()
xs = []
ys = []
for el in track:
xs.append(el[0])
ys.append(el[1])
2019/11/7 Percolation
https://notebooks.dmaitre.phyip3.dur.ac.uk/miscada-sc/user/czjs88/notebooks/Percolation 19/Percolation.ipynb# 7/21
The cells below, test your implementation of methods in Cluster
If your code does not pass this cell, you should debug your implementation of element_add before continuing
# ax.plot(xs,ys,'ro',alpha=0.2)[0]
ax.plot(xs, ys)
# save the figure
if not file == None:
print("Saving track to file: ", file)
plt.savefig(file)
plt.show()
2019/11/7 Percolation
https://notebooks.dmaitre.phyip3.dur.ac.uk/miscada-sc/user/czjs88/notebooks/Percolation 19/Percolation.ipynb# 8/21
In [ ]:
# initialize a small cluster
size = 5 # size of cluster
ndim = 2 # dimensionallity of cluster
periodic = True # boundary condition
cluster = Cluster(ndim=ndim, size=size, periodic=periodic)
# add a single element at the centre
cent = int(size/2) # location of seed
seed = (cent,) * ndim # seed of cluster
cluster.element_add(seed)
# Provided your implementation is correct, cluster.grid = {{2,2)}: True}
test = {(2,2): True}
if cluster.grid == test:
print("That's right! -- continue")
else:
print("you made and error in implementing add_element")
# Using the list of offsets, add the neighbour to the right of the seed to the clust
offset = cluster.offsets[0]
element = tuple(np.asarray(seed) + np.asarray(offset))
print("We are adding element = ", element)
cluster.element_add(element)
test2 = {(2, 2): True, (3, 2): True}
if cluster.grid == test2:
print("That's right! -- continue")
else:
print("you made and error in implementing add_element")
# Test adding an element that is already in the cluster
try:
cluster.element_add(seed)
print(" this should not happen: element should not be added twice ")
except:
print("That's right! -- continue")
pass
# Test periodic wrapping
element = (2*size, int(size/2)) # this element is outside of the grid, but we are u
print("We are adding element = ", element)
cluster.element_add(element)
test3 = {(2, 2): True, (3, 2): True, (0, 2): True}
if cluster.grid == test3:
print("That's right! -- continue")
else:
print("you made and error in implementing add_element (periodic boundary conditi
# now test non-periodic cluster
cluster = Cluster(ndim=ndim, size=size, periodic=False)
try:
cluster.element_add(element)
print(" this should not happen: element should not be added")
except:
print("That's right! -- continue")
pass
2019/11/7 Percolation
https://notebooks.dmaitre.phyip3.dur.ac.uk/miscada-sc/user/czjs88/notebooks/Percolation 19/Percolation.ipynb# 9/21
In [ ]:
The Eden model of cluster growth
Use the methods you have implemented so far, to complete implementing the Eden model. The growth
algortihm is as follows (see lecture notes)
- at start-up, create self.boundary - a list of all valid, empty sites that
are nearest-neighbours of cluster elements
- pick randomly one of the elements in the list of boundary cells, self.boun
dary
- add it to the cluster
- remove it from the list of empty boundary sites
- add its neighbours to the list of boundary sites, provided theya re no
t in the cluster and not already in the list of boundary sites
# This cell test your implementation of the "element_neighbours" method"
# create a small cluster and add some elements to it
size = 7 # size of cluster
ndim = 2 # dimensionallity of cluster
periodic = False # boundary condition
cluster = Cluster(ndim=ndim, size=size, periodic=periodic)
# start with seed at centre, and add all off its neighbours
cent = int(size/2) # location of seed
seed = (cent,) * ndim # seed of cluster
cluster.element_add(seed)
for offset in cluster.offsets:
element = tuple(np.asarray(seed) + np.asarray(offset))
cluster.element_add(element)
# list all neighbours of some positions
neighbours = cluster.element_neighbours(seed)
test = [(4, 3), (2, 3), (3, 4), (3, 2)]
if neighbours == test:
print(" That's right! -- continue")
else:
print("you made an error in implementing element_neighbours")
2019/11/7 Percolation
https://notebooks.dmaitre.phyip3.dur.ac.uk/miscada-sc/user/czjs88/notebooks/Percolation 19/Percolation.ipynb# 10/21
In [ ]:
Fractal dimension of Eden cluster
In the next cell:
Generate an Eden cluster with the parameters given (size, ndim, periodic, and location of starting seed)
Compute the distance of each element to the seed
Sort the distance in ascending order
Plot the enclosed mass as a function of radius. Assume that each element has mass = 1
Compare your answer to Eden_mass.pdf
size = 21 # size of cluster
ndim = 2 # dimensionallity of cluster
periodic = False # boundary condition
cluster = Cluster(ndim=ndim, size=size, random_number_seed = 1, periodic=periodic)
cent = int(size/2)
seed = (cent,) * ndim # seed of cluster
cluster.element_add(seed)
nsites = 200 # add nsites=200 elements to the Eden cluster
for i in range(nsites):
cluster.Eden()
cut = [-1, -1]
file = 'Eden.pdf'
#file = 'Eden_solution.pdf'
title = 'Eden cluster'
cluster.ShowCluster(cut, file=file)
# Compare your result to Eden_solution.pdf
2019/11/7 Percolation
https://notebooks.dmaitre.phyip3.dur.ac.uk/miscada-sc/user/czjs88/notebooks/Percolation 19/Percolation.ipynb# 11/21
In [ ]:
Now generate 2000 realisation of an Eden cluster with the parameters of the cell above. Compute the fractal
dimension , defined as , where is the mass enclosed with radius , for .
Plot a histogram of the values of . Annotate the plot with axes labels, and save it as histogram.pdf
Compute the mean value and dispersion of d. Report your values in the cells below to the required precision.
d m = πr
d m r r = 9
d
size = 101
ndim = 2
periodic = False
cent = int(size/2)
seed = (cent,) * ndim # seed of cluster
nsites = 499
# YOUR CODE HERE
raise NotImplementedError()
# radius is an array of all the distances of each element to seed, sorted in ascenin
# you may want to exclude from it the seed itself
# provide some statistics
print(" max radius =", radius.max())
# unit volume of ndim dimensional sphere
vol = np.pi # special case for ndim = 2
# estimate fractal dimension
dimen = np.log(mass/vol) / np.log(radius)
# plot the result
fig, ax = pl.subplots(1, 2, figsize=(17,7))
ax[0].set_ylim(cluster.ndim-2, cluster.ndim+0.2)
ax[0].set_xlabel('distance to centre')
ax[0].set_ylabel('estimate of fractal dimension')
ax[0].plot(radius, dimen, '.')
ax[0].set_xscale('log')
ax[0].set_title('fractal dimension')
# enclosed mass as a function of radius
ax[1].set_xscale('log')
ax[1].set_yscale('log')
ax[1].plot(radius, mass, '.')
ax[1].set_xlabel('log distance to centre')
ax[1].set_ylabel('log enclosed mass')
ax[1].set_title('enclosed mass as a function of radius')
# overplot mass-radius relation for ndim = 2
mfit2 = vol * radius**(2)
ax[1].plot(radius, mfit2, label='ndim=2')
# overplot mass-radius relation for ndim = 1.9
d = 1.9
mfitd = vol * radius**(d)
ax[1].plot(radius, mfitd, label='ndim=1.9')
ax[1].legend()
# plt.savefig("Eden_mass.pdf")
plt.show()
2019/11/7 Percolation
https://notebooks.dmaitre.phyip3.dur.ac.uk/miscada-sc/user/czjs88/notebooks/Percolation 19/Percolation.ipynb# 12/21
p p p y q p
In [ ]:
Compute and plot a histogram of estimates of the fractal dimension.
2 marks
In [ ]:
Compute the mean (2 marks) and standard deviation (2 marks) of the fractal dimension. Enter your values in
the boxes below. Execute the hidden cells below if the boxes do not appear.
The DLA model of cluster growth
Use the methods you have implemented so far, to complete the implemention of the DLA model of cluster
growth (see course notes). In brief, the growth algorithm is as follows:
create a cluster with a single seed at the centre, using the Cluster class you've developped
choose an inner radius, , and an outer radius .
start a random walker from a location at distance close to from the seed, and let it walk
if it hits any part of the existing cluster, add it to the cluster, and start a new walk
if it crosses the outer boundary, , terminate the random walk, and start a new walk
terminate the programme, when the cluster contains a specified number of elements
You may want to base the random walk on the implementation you wrote in lecture 6 on random walks
Use the cells below to test intermediate steps of this assignment
r1 r2
r1
r > r2
# YOUR CODE HERE
raise NotImplementedError()
print("Calculation finished")
# use this cell to
# -- compute a histgram of the dimensions "dimen", estimate for each of the nEden Ed
# -- plot this histogram and save it to a file, called Eden_fractal.pdf
# Annote the histogram: add labels to the axes
# YOUR CODE HERE
raise NotImplementedError()
2019/11/7 Percolation
https://notebooks.dmaitre.phyip3.dur.ac.uk/miscada-sc/user/czjs88/notebooks/Percolation 19/Percolation.ipynb# 13/21
In [ ]:
In [ ]:
Now create a DLA cluster with nel = 400 elements in it
Use the cell below it plot the DLA cluster and saving the plot to a file.
Take care with the parameter 'radius' and 'edge'. The values in the cell below are reasonable initialize values
from where to start the walk (at a distance "radius" from the centre of the DLA cluster), and for when to abort
the random walk (if the walker ventures further than "edge' from the centre of the cluster). However, as the
cluster grows in size, you will need to increase both "radius" and "edge".
- If radius is too small, the centre of the cluster will be densely filled,
since all walkers are started too close in.
- If edge is too small, you again get a cluster that is too centrally concen
trated.
- If edge and or radius are too large, then the code becomes very inefficien
t, because the walker can run everywhere with little chance of ever hitting
the cluster.
Therefore, you may want to experiment with both these parameters.
4 marks
# create a cluster, and add a seed to the centre
size = 101 # the size here is only ever used for plotting
ndim = 2 # we only do the 2 dimensional case
periodic = False # non-periodic grid
cent = int(size/2) # put the seed in the center of the grid
seed = (cent,) * ndim # seed of cluster
# We specify a random number seed for this test, so you should always get the same s
dla = Cluster(ndim=ndim, size=size, random_number_seed=12, periodic=periodic)
dla.element_add(seed)
# as an example, perform a single walk to see what happens
centre = seed
radius = 20
edge = 40
track = dla.DLA(centre, radius, edge)
#file = "Randomwalk_solution.pdf"
file = "Randomwalk.pdf"
# the call below plots the track the random walk followed, before it terminated
# - either by travelling further than 'edge' from the seed
# - or by hitting the seed
dla.ShowTrack(track, file=file)
# you may want to compare your answer to the file "Randomwalk_solution.pdf"
# It should not look identical, but if it looks very different you may have a bug so
2019/11/7 Percolation
https://notebooks.dmaitre.phyip3.dur.ac.uk/miscada-sc/user/czjs88/notebooks/Percolation 19/Percolation.ipynb# 14/21
In [ ]:
In [ ]:
Percolation
The cell defines a class GroupFinder that identifies groups in a cluster.
The structure of the cluster required by GroupFinder is the same as what we had so far. The Eden cluster, or
DLA clusters, have the right format. Given a cluster, the groupfinder partitions the elements in a cluster by
making sure that any two elements of the cluster that are nearest-neighbours, belong to the same group.
# Use this cell to compute the DLA cluster
# The use the next cell to plot it and save it to a file called "DLA.pdf"
# You may want to use the ShowCluster method for plotting, with cut=[-1,-1]
# start with these parameters for the cluster
size = 101 # size is only relevant for plotting and setting the
cent = int(size/2)
seed = (cent,) * ndim # seed of cluster
periodic = False
dla = Cluster(ndim=ndim, size=size, random_number_seed=12, periodic=periodic)
dla.element_add(seed)
# name of file to save result, and location of cut to plot
cut = [-1, -1]
file = 'DLA.pdf'
#
centre = seed
radius = 10
edge = 20
nel = 400
# YOUR CODE HERE
raise NotImplementedError()
print("Calculation is finished")
# YOUR CODE HERE
raise NotImplementedError()
2019/11/7 Percolation
https://notebooks.dmaitre.phyip3.dur.ac.uk/miscada-sc/user/czjs88/notebooks/Percolation 19/Percolation.ipynb# 15/21
In [ ]:
# This class contains methods to
# - identify all groups in a given cluster
# a group is a set of sites that are neaarest-neighbours
# - examine whether a group percolates
# You should not have to edit any of these methods
# Use the cells below to learn how to use the method to identify groups, and examine
class GroupFinder():
# return a list of connected clusters
def __init__(self, cluster):
self.cluster = cluster
# Each element within a group will point to a common element of the group -
self.heads = {}
for element in self.cluster.grid:
self.heads[element] = element # initialize by assuming each element is
# Iteratively loop over all elements in the cluster to identify groups
def FoF(self):
# Begin partitioning
# We only loop over pairs of neighbours
maxcount = 10000
for element in self.cluster.grid:
element_head = element
count = 0
while element_head != self.heads[element_head] and count < maxcount:
element_head = self.heads[element_head]
count = 0
for neighbour in self.cluster.element_neighbours(element)[::2]:
count = 0
if neighbour in self.cluster.grid:
neighbour_head = neighbour
while neighbour_head != self.heads[neighbour_head] and count < m
count += 1
neighbour_head = self.heads[neighbour_head]
self.heads[neighbour_head] = element_head
if count == maxcount:
raise Exception("Error: partitioning does not converge")
#
count = 0
for element in self.cluster.grid:
element_head = element
while element_head != self.heads[element_head] and count < maxcount:
element_head = self.heads[element_head]
self.heads[element] = element_head
unique = {}
ngroups = 0
for head in self.heads:
if head not in unique:
unique[head] = ngroups
ngroups += 1
groupid = []
for element in self.cluster.grid:
groupid.append(unique[self.heads[element]])
return np.asarray(groupid)
2019/11/7 Percolation
https://notebooks.dmaitre.phyip3.dur.ac.uk/miscada-sc/user/czjs88/notebooks/Percolation 19/Percolation.ipynb# 16/21
# Alternative group finder
def Groups(self):
groupid = np.arange(len(self.cluster.grid), dtype=int)
groupindx = {}
for index, element in enumerate(self.cluster.grid):
groupindx[element] = index
success = False
count = 0
while not success:
success = True
for element in self.cluster.grid:
for neighbour in self.cluster.element_neighbours(element):
if neighbour in self.cluster.grid:
indx1 = groupid[groupindx[element]]
indx2 = groupid[groupindx[neighbour]]
if indx1 < indx2:
groupid[groupid==indx2] = indx1
success = False
count += 1
elif indx1 > indx2:
groupid[groupid==indx1] = indx2
success = False
count += 1
print(" count = ", count)
return groupid
# Examine the percolation of groups identified using a group finder
# Loop over all unique groups in groupdid
# for each group, loop over all dimensions
# if a group percolates, store its id and the dimension that it percoaltes i
# return this to the calling routine
def PercolateGrid(self, groupid):
# determine whether thisgroup percolates
# consistency: only works for non-periodic grids
if self.cluster.periodic:
raise Exception(" We only implement non-periodic grids ")
#
groupindx = {}
for index, element in enumerate(self.cluster.grid):
groupindx[element] = index
# extract group by group
result = []
for group in np.unique(groupid):
newcluster = Cluster(self.cluster.ndim, self.cluster.size, self.cluster.
for element in self.cluster.grid:
if groupid[groupindx[element]] == group:
newcluster.grid[element] = True
# find whether this group percolates
for dimen in range(self.cluster.ndim):
if self.PercolateGroup(newcluster, dimen):
result.append([group, dimen])
return result
# This method examines whether a given group percolates in a given dimension
def PercolateGroup(self, cluster, dimen):
# inputs:
2019/11/7 Percolation
https://notebooks.dmaitre.phyip3.dur.ac.uk/miscada-sc/user/czjs88/notebooks/Percolation 19/Percolation.ipynb# 17/21
Examine the cell to see how the methods of GroupFinder and Cluster can be
used to
create a cluster to work on
identify the groups inside the cluster
plot the result
# self: the class GroupFinder
# cluster: one of the groups identified by the group finder
# dimen: the dimension to check for percolation
# consistency check
if dimen < 0 or dimen >= cluster.ndim:
raise Exception(" too many dimensions ")
# Determine with this cluster touches both boundaries in dimension dimen of
x = [sp[dimen] for sp in cluster.grid.keys()]
return (min(x) == 0 and max(x) == cluster.size-1)
####### The routines below are ueful for plotting groups #################
def ExtractGroup(self, groupid, plotid, cut):
groupindx = {}
for index, element in enumerate(self.cluster.grid):
groupindx[element] = index
newcluster = Cluster(self.cluster.ndim, self.cluster.size, self.cluster.peri
for element in self.cluster.grid:
if groupid[groupindx[element]] == plotid:
newcluster.grid[element] = True
xs, ys = newcluster.ExtractSlice(cut)
return xs, ys
def ShowGroup(self, groupid, plotid, cut):
xs, ys = self.ExtractGroup(groupid, plotid, cut)
ax.plot(xs, ys, 'o',alpha=0.2)[0]
2019/11/7 Percolation
https://notebooks.dmaitre.phyip3.dur.ac.uk/miscada-sc/user/czjs88/notebooks/Percolation 19/Percolation.ipynb# 18/21
In [ ]:
# set-up the grid parameters
size = 512
ndim = 2
periodic = True
rseed = 12
cluster = Cluster(size=size, ndim=ndim, periodic=periodic, random_number_seed = rs
# Use the Eden model to grow some clusters
# start woth 5 seperate seeds
seed1 = (50,50)
seed2 = (50,200)
seed3 = (100,50)
seed4 = (200,400)
seed5 = (300,200)
cluster.element_add(seed1)
cluster.element_add(seed2)
cluster.element_add(seed3)
cluster.element_add(seed4)
cluster.element_add(seed5)
# run the Eden model, growing each seed
for i in range(10000):
cluster.Eden()
# Identify the groups
groupfinder = GroupFinder(cluster)
# Assign each element to its group. Elements of the same group have the same groupid
groupid = groupfinder.FoF()
grouplen = []
# Next, count the number of elements in each seperate group
for id in np.unique(groupid):
grouplen.append(len(np.where(groupid == id)[0]))
print(" Number of elements in groups = ", grouplen)
# Print some statistics
print(" Number of unique groups = ", len(np.unique(groupid)))
print(" Number of elements in largest group = ", max(grouplen))
# Now plot all groups, using different colours for each group
fig, ax = pl.subplots(1, 1, figsize=(7,7))
ax.set_aspect('equal')
ax.set_xlim(-1, cluster.size)
ax.set_ylim(-1, cluster.size)
ax.set_xlabel('x-position')
ax.set_ylabel('y-position')
ax.set_title(' connected gro
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
