Created
May 18, 2012 15:22
-
-
Save tpoisot/2725839 to your computer and use it in GitHub Desktop.
Metapop model using networkx
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
#!/usr/bin/env python | |
# encoding: utf-8 | |
import networkx as nx | |
import numpy as np | |
import matplotlib.pyplot as plt | |
Patches = 100 # Number of patches | |
P_ext = 0.01 # Probability of extinction (e) | |
P_col = 0.014 # Probability of colonization (c) | |
P_init = 0.02 # Probability that a patch will be occupied at the beginning | |
Distance = 1.4 # An arbitrary parameter to determine which patches are connected | |
class patch: | |
def __init__(self,status=0,pos=(0,0)): | |
self.status = status | |
self.pos = pos | |
def __str__(self): | |
return(str(self.status)) | |
G = nx.Graph() | |
for i in xrange(Patches): | |
Stat = 1 if np.random.uniform() < P_init else 0 | |
Pos = (np.random.uniform()*10-5,np.random.uniform()*10-5) | |
G.add_node(patch(Stat,Pos)) | |
for p1 in G.nodes(): | |
for p2 in G.nodes(): | |
Dist = np.sqrt((p1.pos[1]-p2.pos[1])**2+(p1.pos[0]-p2.pos[0])**2) | |
if Dist <= Distance: | |
G.add_edge(p1,p2) | |
pos = {} | |
for n in G.nodes(): | |
pos[n] = n.pos | |
occup = [n.status for n in G] | |
nx.draw(G,node_color=occup,with_labels=False,cmap=plt.cm.Greys,pos=pos,vmin=0,vmax=1) | |
plt.show() | |
Time = [0] | |
Occupancy = [np.sum([n.status for n in G])/float(Patches)] | |
for timestep in xrange(2000): | |
## Check for extinctions | |
for n in G.nodes(): | |
if (n.status == 1 and np.random.uniform() < P_ext): | |
n.status = 0 | |
## Check for invasions | |
for n in G.nodes(): | |
if n.status == 1: | |
neighb = G[n] # That's it, a list of the neighbors | |
for nei in neighb: | |
if nei.status == 0: | |
if np.random.uniform() < P_col: | |
nei.status = 1 | |
break | |
Time.append(timestep+1) | |
Occupancy.append(np.sum([n.status for n in G])/float(Patches)) | |
nx.draw(G,node_color=[n.status for n in G],with_labels=False,cmap=plt.cm.Greys,pos=pos,vmin=0,vmax=1) | |
plt.show() | |
plt.plot(Time,Occupancy,'g-') | |
plt.show() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment