In [1]:
import numpy as np
import networkx as nx
import matplotlib.pyplot as plt
%matplotlib inline
import scipy.optimize as opt
Load the adjacency matrix (csv file to numpy array):
In [5]:
A = np.loadtxt('../data/airport/air500matrix.csv', dtype=int, delimiter=',')
Load the airport names (IATA codes) and lat/long information:
In [39]:
latlong_data=np.loadtxt('../data/airport/latlongs.dat', dtype='str', delimiter=',')
node_labels=[str(i[2:-1]).strip() for i in latlong_data[:,0]]
node_pos=[(float(i[3:-2]), float(j[3:-1])) for i,j in latlong_data[:,1:]]
Create a networkx graph from the adjacency matrix:
In [43]:
T = nx.Graph(A)
Label the nodes with the airport names:
In [44]:
N = len(node_labels)
label_mapping = dict(zip(range(N), node_labels))
T = nx.relabel_nodes(T, label_mapping)
In [45]:
nx.set_node_attributes(T, "latlon", dict(zip(node_labels, node_pos)))
In [47]:
nx.write_gpickle(T, "../data/airport/big_airportnet.gpickle")
In [20]:
T.node['FRA']['latlon']
Out[20]:
In [21]:
nx.draw(T)
In [22]:
nx.draw_networkx(T, with_labels=False, node_size=10, alpha=0.6, width=0.1)
In [23]:
nx.draw_networkx(T, with_labels=False, node_size=10, alpha=0.6, width=0.07, pos=nx.random_layout(T))
There are many other layouts available. please expore nx.layouts
In [24]:
from mpl_toolkits.basemap import Basemap
import matplotlib.colors as colors
In [25]:
bmap=Basemap(llcrnrlon=-170, llcrnrlat=-48.7,urcrnrlon=179.99, urcrnrlat=70.8, lat_0=45, lon_0=10, resolution='h')
In [26]:
def plot_over_basemap(G, node_colors=None):
fig = plt.figure(figsize=(25,25), dpi=100)
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8])
dummy=bmap.drawcoastlines(color='white')
dummy2=bmap.drawcountries(color='white')
dummy3=bmap.fillcontinents(color='0.8', alpha=0.6)
node_labels=G.nodes()
pos_on_map=[bmap(*G.node[label]['latlon']) for label in node_labels]
posdict=dict(zip(node_labels, pos_on_map))
if node_colors is not None:
nodes=nx.draw_networkx_nodes(G, with_labels=False, node_size=25, alpha=1.0, width=0.1, pos=posdict, node_color=node_colors, cmap=plt.get_cmap('cool'))
nx.draw_networkx_edges(G, pos=posdict, alpha=0.4)
plt.colorbar(nodes, orientation='horizontal', pad=0.001, aspect=40)
plt.axis('off')
else:
nx.draw_networkx(G, with_labels=False, node_size=25, alpha=1.0, width=0.1, pos=posdict)
In [27]:
plot_over_basemap(T)
This is way too big. Let's consider only the airports that fall in a lat/lon box:
In [5]:
airports_eurasia=np.loadtxt("data/eurasia_airport_labels.txt",dtype=str)
In [29]:
G=T.subgraph(airports_eurasia)
In [30]:
bmap=Basemap(llcrnrlon=-18.32, llcrnrlat=-49,urcrnrlon=179.99, urcrnrlat=71, lat_0=45, lon_0=10, resolution='h')
In [31]:
plot_over_basemap(G)
In [32]:
print G.degree()
In [33]:
G.degree()['LHR']
Out[33]:
In [34]:
deg_dist=G.degree().values()
In [35]:
freqs,bins,_=plt.hist(deg_dist)
In [36]:
freqs
Out[36]:
In [37]:
bins
Out[37]:
In [38]:
binmids=[(bins[i]+bins[i+1])/2.0 for i in range(freqs.shape[0])]
In [39]:
binmids
Out[39]:
In [40]:
powerfunc=lambda x,a,l:a*x**(-l)
In [41]:
a,l=opt.curve_fit(powerfunc,binmids[2:],freqs[2:])[0]
In [42]:
plt.plot(binmids[2:], powerfunc(binmids[2:], a,l)); plt.hist(deg_dist)
Out[42]:
In [43]:
nx.shortest_path(G, source='VIE', target='CCU')
Out[43]:
But wait! This is not the only shortest path. Fortunately, we are covered by another function
In [44]:
for path in nx.all_shortest_paths(G, source='VIE', target='CCU'):
print path
In [45]:
cents=nx.closeness_centrality(G)
In [46]:
centarray=[cents[node] for node in G.nodes()]
plot_over_basemap(G, node_colors=centarray)
In [47]:
plt.hist(centarray)
Out[47]:
Let the initial concentration of pathogen be:
$c_0=(c_1,c_2,\cdots,c_n)$
And let the spreading happen according to this law:
$c_{n+1}=(D^{-1}A)^T c_n=Sc_n$
In [48]:
A_d=np.array(nx.adjacency_matrix(G))
D=A_d.sum(axis=0)
Dinv=np.diag(np.where(D>0,1.0/D,0))
S=np.dot(Dinv,A_d).T
In [49]:
conc=np.zeros(G.number_of_nodes())
conc[G.nodes().index('FRA')]=1
In [50]:
S
Out[50]:
In [51]:
plot_over_basemap(G, node_colors=conc)
In [52]:
conc2=np.dot(S,conc)
plot_over_basemap(G, node_colors=conc2)
In [53]:
conc3=np.dot(S,conc2)
plot_over_basemap(G, node_colors=conc3)
In [54]:
plt.hist(conc3)
Out[54]:
In [55]:
plt.hist(conc2)
Out[55]:
In [56]:
plt.hist(conc)
Out[56]:
In [57]:
S50=np.linalg.matrix_power(S,50)
In [58]:
conc50=np.dot(S50, conc)
In [59]:
conc100=np.dot(np.linalg.matrix_power(S,50), conc50)
In [60]:
conc150=np.dot(np.linalg.matrix_power(S,50), conc100)
In [61]:
np.linalg.norm(conc150-conc100)
Out[61]:
In [62]:
plot_over_basemap(G, node_colors=conc50)
In [63]:
plt.hist(conc50)
Out[63]:
In [64]:
conc50[G.nodes().index('FRA')]
Out[64]:
In [65]:
H=G.copy()
In [66]:
edge_centrs=nx.edge_betweenness_centrality(H)
In [67]:
sorted_centrs=sorted(edge_centrs.keys(), cmp=lambda x,y:cmp(edge_centrs[x], edge_centrs[y]))
In [68]:
for i in xrange(10):
rmstart, rmend=sorted_centrs[-i]
H.remove_edge(rmstart, rmend)
In [69]:
A_d=np.array(nx.adjacency_matrix(H))
D=A_d.sum(axis=0)
Dinv=np.diag(np.where(D>0,1.0/D,0))
S=np.dot(Dinv,A_d).T
In [70]:
conc=np.zeros(H.number_of_nodes())
conc[G.nodes().index('FRA')]=1
In [71]:
conc50=np.dot(np.linalg.matrix_power(S,50), conc)
In [72]:
plt.hist(conc50)
Out[72]:
In [73]:
conc2=np.dot(S, conc)
In [74]:
plt.hist(conc2)
Out[74]:
In [75]:
plot_over_basemap(H, node_colors=conc2)
In [76]:
plt.hist(conc50)
Out[76]: