首页 >
> 详细

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

联系我们

- QQ：99515681
- 邮箱：99515681@qq.com
- 工作时间：8:00-23:00
- 微信：codinghelp2

- Tsp课程作业代写、代做algorithms留学生作业、代做java，C/C 2020-06-23
- Kit107留学生作业代做、C++编程语言作业调试、Data课程作业代写、代 2020-06-23
- Sta302h1f作业代做、代写r课程设计作业、代写r编程语言作业、代做da 2020-06-22
- 代写seng 474作业、代做data Mining作业、Python，Ja 2020-06-22
- Cmpsci 187 Binary Search Trees 2020-06-21
- Comp226 Assignment 2: Strategy 2020-06-21
- Math 504 Homework 12 2020-06-21
- Math4007 Assessed Coursework 2 2020-06-21
- Optimization In Machine Learning Assig... 2020-06-21
- Homework 1 – Math 104B 2020-06-20
- Comp1000 Unix And C Programming 2020-06-20
- General Specifications Use Python In T... 2020-06-20
- Comp-206 Mini Assignment 6 2020-06-20
- Aps 105 Lab 9: Search And Link 2020-06-20
- Aps 105 Lab 9: Search And Link 2020-06-20
- Mech 203 – End-Of-Semester Project 2020-06-20
- Ms980 Business Analytics 2020-06-20
- Cs952 Database And Web Systems Develop... 2020-06-20
- Homework 4 Using Data From The China H... 2020-06-20
- Assignment 1 Build A Shopping Cart 2020-06-20