Created
July 18, 2013 21:39
-
-
Save cameronneylon/6033364 to your computer and use it in GitHub Desktop.
Code for percolation network simulation used in my talk at #BOSC July 19 2013 in Berlin. Based on code from http://dragly.org/2013/03/25/working-with-percolation-clusters-in-python/. The license for the source code is not entirely clear so this code is made available under a CC BY-SA license as per the blog it was derived from. For reasons not e…
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
from pylab import * | |
from scipy.ndimage import measurements | |
import numpy as np | |
import matplotlib.pyplot as plt | |
import time | |
import sys | |
import signal | |
def signal_handler(signal, frame): | |
sys.exit(0) | |
signal.signal(signal.SIGINT, signal_handler) | |
# Set up the figures | |
plt.figure(0, figsize=(9,9)) | |
plt.title("Clusters by area") | |
plt.figure(1, figsize=(7,4)) | |
plt.title("Number of clusters") | |
plt.figure(2, figsize = (7,4)) | |
plt.title("Largest Cluster") | |
plt.ion() | |
time.sleep(20) | |
for L in [50, 100, 200, 500, 1000]: | |
for marker in ['ko', 'bo', 'ro']: | |
# Randomise the percolation network | |
r = rand(L,L) | |
for p in np.arange(0., 1., 0.05): | |
# Determine the connectivity for probability p | |
z = r<p | |
# Define the clusters using the ndimage.measurements routines | |
lw, num = measurements.label(z) | |
# Plot the clusters for given p | |
plt.figure(0) | |
area = measurements.sum(z, lw, index=arange(lw.max() + 1)) | |
areaImg = area[lw] | |
plt.imshow(areaImg, origin='lower', interpolation='nearest', vmin=1, vmax=L*5) | |
plt.draw() | |
# Plot the number of clusters | |
plt.figure(1) | |
plt.plot(p, lw.max(), marker) | |
plt.draw() | |
# Plot the size of the largest cluster | |
plt.figure(2) | |
plt.plot(p, area.max(), marker) | |
plt.draw() | |
plt.ioff() | |
plt.show() | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment