Chapter 10: Sparse matrices and graphs

Robert Johansson

Source code listings for Numerical Python - A Practical Techniques Approach for Industry (ISBN 978-1-484205-54-9).

The source code listings can be downloaded from http://www.apress.com/9781484205549

In [1]:
%matplotlib inline
%config InlineBackend.figure_format='retina'
In [2]:
import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.rcParams['text.usetex'] = True
In [3]:
import scipy.sparse as sp
In [4]:
import scipy.sparse.linalg
In [5]:
import numpy as np
In [6]:
import scipy.linalg as la
In [7]:
import networkx as nx

Coordinate list format

In [8]:
values = [1, 2, 3, 4]
In [9]:
rows = [0, 1, 2, 3]
In [10]:
cols = [1, 3, 2, 0]
In [11]:
A = sp.coo_matrix((values, (rows, cols)), shape=[4, 4])
In [12]:
A.todense()
Out[12]:
matrix([[0, 1, 0, 0],
        [0, 0, 0, 2],
        [0, 0, 3, 0],
        [4, 0, 0, 0]])
In [13]:
A
Out[13]:
<4x4 sparse matrix of type '<type 'numpy.int64'>'
	with 4 stored elements in COOrdinate format>
In [14]:
A.shape, A.size, A.dtype, A.ndim
Out[14]:
((4, 4), 4, dtype('int64'), 2)
In [15]:
A.nnz, A.data
Out[15]:
(4, array([1, 2, 3, 4]))
In [16]:
A.row
Out[16]:
array([0, 1, 2, 3], dtype=int32)
In [17]:
A.col
Out[17]:
array([1, 3, 2, 0], dtype=int32)
In [18]:
A.tocsr()
Out[18]:
<4x4 sparse matrix of type '<type 'numpy.int64'>'
	with 4 stored elements in Compressed Sparse Row format>
In [19]:
A.toarray()
Out[19]:
array([[0, 1, 0, 0],
       [0, 0, 0, 2],
       [0, 0, 3, 0],
       [4, 0, 0, 0]])
In [20]:
A.todense()
Out[20]:
matrix([[0, 1, 0, 0],
        [0, 0, 0, 2],
        [0, 0, 3, 0],
        [4, 0, 0, 0]])

Not all sparse matrix formats supports indexing:

In [21]:
A[1,2]
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-21-13df355f18af> in <module>()
----> 1 A[1,2]

TypeError: 'coo_matrix' object has no attribute '__getitem__'
In [22]:
A.tobsr()[1,2]
---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
<ipython-input-22-7f4c44f6556f> in <module>()
----> 1 A.tobsr()[1,2]

/Users/rob/miniconda/envs/py27-npm/lib/python2.7/site-packages/scipy/sparse/bsr.pyc in __getitem__(self, key)
    297 
    298     def __getitem__(self,key):
--> 299         raise NotImplementedError
    300 
    301     def __setitem__(self,key,val):

NotImplementedError: 

But some do:

In [23]:
A.tocsr()[1,2]
Out[23]:
0
In [24]:
A.tolil()[1:3,3]
Out[24]:
<2x1 sparse matrix of type '<type 'numpy.int64'>'
	with 1 stored elements in LInked List format>

CSR

In [25]:
A = np.array([[1, 2, 0, 0], [0, 3, 4, 0], [0, 0, 5, 6], [7, 0, 8, 9]]); A
Out[25]:
array([[1, 2, 0, 0],
       [0, 3, 4, 0],
       [0, 0, 5, 6],
       [7, 0, 8, 9]])
In [26]:
A = sp.csr_matrix(A)
In [27]:
A.data
Out[27]:
array([1, 2, 3, 4, 5, 6, 7, 8, 9])
In [28]:
A.indices
Out[28]:
array([0, 1, 1, 2, 2, 3, 0, 2, 3], dtype=int32)
In [29]:
A.indptr
Out[29]:
array([0, 2, 4, 6, 9], dtype=int32)
In [30]:
i = 2
In [31]:
A.indptr[i], A.indptr[i+1]-1
Out[31]:
(4, 5)
In [32]:
A.indices[A.indptr[i]:A.indptr[i+1]]
Out[32]:
array([2, 3], dtype=int32)
In [33]:
A.data[A.indptr[i]:A.indptr[i+1]]
Out[33]:
array([5, 6])

Functions for constructing sparse matrices

In [34]:
N = 10
In [35]:
A = -2 * sp.eye(N) + sp.eye(N, k=1) + sp.eye(N, k=-1)
In [36]:
A
Out[36]:
<10x10 sparse matrix of type '<type 'numpy.float64'>'
	with 28 stored elements in Compressed Sparse Row format>
In [37]:
A.todense()
Out[37]:
matrix([[-2.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
        [ 1., -2.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
        [ 0.,  1., -2.,  1.,  0.,  0.,  0.,  0.,  0.,  0.],
        [ 0.,  0.,  1., -2.,  1.,  0.,  0.,  0.,  0.,  0.],
        [ 0.,  0.,  0.,  1., -2.,  1.,  0.,  0.,  0.,  0.],
        [ 0.,  0.,  0.,  0.,  1., -2.,  1.,  0.,  0.,  0.],
        [ 0.,  0.,  0.,  0.,  0.,  1., -2.,  1.,  0.,  0.],
        [ 0.,  0.,  0.,  0.,  0.,  0.,  1., -2.,  1.,  0.],
        [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  1., -2.,  1.],
        [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1., -2.]])
In [38]:
fig, ax = plt.subplots()
ax.spy(A)
fig.savefig("ch10-sparse-matrix-1.pdf");
In [39]:
A = sp.diags([1,-2,1], [1,0,-1], shape=[N, N], format='csc')
In [40]:
A
Out[40]:
<10x10 sparse matrix of type '<type 'numpy.float64'>'
	with 28 stored elements in Compressed Sparse Column format>
In [41]:
fig, ax = plt.subplots()
ax.spy(A);
In [42]:
B = sp.diags([1, 1], [-1, 1], shape=[3,3])
In [43]:
B
Out[43]:
<3x3 sparse matrix of type '<type 'numpy.float64'>'
	with 4 stored elements (2 diagonals) in DIAgonal format>
In [44]:
C = sp.kron(A, B, format='csr')
C
Out[44]:
<30x30 sparse matrix of type '<type 'numpy.float64'>'
	with 112 stored elements in Compressed Sparse Row format>
In [45]:
fig, (ax_A, ax_B, ax_C) = plt.subplots(1, 3, figsize=(12, 4))
ax_A.spy(A)
ax_B.spy(B)
ax_C.spy(C)
fig.savefig("ch10-sparse-matrix-2.pdf");

Sparse linear algebra

In [46]:
N = 10
In [47]:
A = sp.diags([1, -2, 1], [1, 0, -1], shape=[N, N], format='csc')
In [48]:
b = -np.ones(N)
In [49]:
x = sp.linalg.spsolve(A, b)
In [50]:
x
Out[50]:
array([  5.,   9.,  12.,  14.,  15.,  15.,  14.,  12.,   9.,   5.])
In [51]:
np.linalg.solve(A.todense(), b)
Out[51]:
array([  5.,   9.,  12.,  14.,  15.,  15.,  14.,  12.,   9.,   5.])
In [52]:
lu = sp.linalg.splu(A)
In [53]:
lu.L
Out[53]:
<10x10 sparse matrix of type '<type 'numpy.float64'>'
	with 20 stored elements in Compressed Sparse Column format>
In [54]:
lu.perm_r
Out[54]:
array([0, 1, 2, 3, 4, 5, 6, 8, 7, 9], dtype=int32)
In [55]:
lu.U
Out[55]:
<10x10 sparse matrix of type '<type 'numpy.float64'>'
	with 20 stored elements in Compressed Sparse Column format>
In [56]:
lu.perm_c
Out[56]:
array([0, 1, 2, 3, 4, 5, 6, 8, 7, 9], dtype=int32)
In [57]:
def sp_permute(A, perm_r, perm_c):
    """ permute rows and columns of A """
    M, N = A.shape
    # row permumation matrix
    Pr = sp.coo_matrix((np.ones(M), (perm_r, np.arange(N)))).tocsr()
    # column permutation matrix
    Pc = sp.coo_matrix((np.ones(M), (np.arange(M), perm_c))).tocsr()
    return Pr.T * A * Pc.T
In [58]:
lu.L * lu.U - A
Out[58]:
<10x10 sparse matrix of type '<type 'numpy.float64'>'
	with 8 stored elements in Compressed Sparse Column format>
In [59]:
sp_permute(lu.L * lu.U, lu.perm_r, lu.perm_c) - A
Out[59]:
<10x10 sparse matrix of type '<type 'numpy.float64'>'
	with 0 stored elements in Compressed Sparse Column format>
In [60]:
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(12, 4))
ax1.spy(lu.L)
ax2.spy(lu.U)
ax3.spy(A)
Out[60]:
<matplotlib.lines.Line2D at 0x10985a810>
In [61]:
x = lu.solve(b)
In [62]:
x
Out[62]:
array([  5.,   9.,  12.,  14.,  15.,  15.,  14.,  12.,   9.,   5.])
In [63]:
# use_umfpack=True is only effective if scikit-umfpack is installed
# (in which case UMFPACK is the default solver)
x = sp.linalg.spsolve(A, b, use_umfpack=True)
In [64]:
x
Out[64]:
array([  5.,   9.,  12.,  14.,  15.,  15.,  14.,  12.,   9.,   5.])
In [65]:
x, info = sp.linalg.cg(A, b)
In [66]:
x
Out[66]:
array([  5.,   9.,  12.,  14.,  15.,  15.,  14.,  12.,   9.,   5.])
In [67]:
x, info = sp.linalg.bicgstab(A, b)
In [68]:
x
Out[68]:
array([  5.,   9.,  12.,  14.,  15.,  15.,  14.,  12.,   9.,   5.])
In [69]:
x, info = sp.linalg.lgmres(A, b)
In [70]:
x
Out[70]:
array([  5.,   9.,  12.,  14.,  15.,  15.,  14.,  12.,   9.,   5.])
In [71]:
N = 25
In [72]:
A = sp.diags([1, -2, 1], [8, 0, -8], shape=[N, N], format='csc')

An example of a matrix reording method: Reverse Cuthil McKee

In [73]:
perm = sp.csgraph.reverse_cuthill_mckee(A)
perm
Out[73]:
array([ 7, 15, 23,  1,  9, 17,  2, 10, 18,  3, 11, 19,  4, 12, 20,  5, 13,
       21,  6, 14, 22, 24, 16,  8,  0], dtype=int32)
In [74]:
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4))
ax1.spy(A)
ax2.spy(sp_permute(A, perm, perm))
Out[74]:
<matplotlib.lines.Line2D at 0x105c2cf10>

Performance comparison sparse/dense

In [75]:
# compare performance of solving Ax=b vs system size N,
# where A is the sparse matrix for the 1d poisson problem
import time

def setup(N):
    A = sp.diags([1,-2,1], [1,0,-1], shape=[N, N], format='csr')
    b = -np.ones(N)
    return A, A.todense(), b

reps = 10
N_vec = np.arange(2, 300, 1)
t_sparse = np.empty(len(N_vec))
t_dense = np.empty(len(N_vec))
for idx, N in enumerate(N_vec):
    A, A_dense, b = setup(N)
    t = time.time()
    for r in range(reps):
        x = np.linalg.solve(A_dense, b)
    t_dense[idx] = (time.time() - t)/reps
    t = time.time()
    for r in range(reps):
        #x = la.solve(A_dense, b)
        x = sp.linalg.spsolve(A, b, use_umfpack=True)
    t_sparse[idx] = (time.time() - t)/reps
    
fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(N_vec, t_dense * 1e3, '.-', label="dense")
ax.plot(N_vec, t_sparse * 1e3, '.-', label="sparse")
ax.set_xlabel(r"$N$", fontsize=16)
ax.set_ylabel("elapsed time (ms)", fontsize=16)
ax.legend(loc=0)
fig.tight_layout()
fig.savefig("ch10-sparse-vs-dense.pdf")

Eigenvalue problems

In [76]:
N = 10
In [77]:
A = sp.diags([1, -2, 1], [1, 0, -1], shape=[N, N], format='csc')
In [78]:
evals, evecs = sp.linalg.eigs(A, k=4, which='LM')
In [79]:
evals
Out[79]:
array([-3.91898595+0.j, -3.68250707+0.j, -3.30972147+0.j, -2.83083003+0.j])
In [80]:
np.allclose(A.dot(evecs[:,0]), evals[0] * evecs[:,0])
Out[80]:
True
In [81]:
evals, evecs = sp.linalg.eigsh(A, k=4, which='LM')
In [82]:
evals
Out[82]:
array([-3.91898595, -3.68250707, -3.30972147, -2.83083003])
In [83]:
evals, evecs = sp.linalg.eigs(A, k=4, which='SR')
In [84]:
evals
Out[84]:
array([-3.91898595+0.j, -3.68250707+0.j, -3.30972147+0.j, -2.83083003+0.j])
In [85]:
np.real(evals).argsort()
Out[85]:
array([0, 1, 2, 3])
In [86]:
def sp_eigs_sorted(A, k=6, which='SR'):
    """ compute and return eigenvalues sorted by real value """
    evals, evecs = sp.linalg.eigs(A, k=k, which=which)
    idx = np.real(evals).argsort()
    return evals[idx], evecs[idx]
In [87]:
evals, evecs = sp_eigs_sorted(A, k=4, which='SM')
In [88]:
evals
Out[88]:
array([-1.16916997+0.j, -0.69027853+0.j, -0.31749293+0.j, -0.08101405+0.j])

Random matrix example

In [89]:
N = 100
In [90]:
x_vec = np.linspace(0, 1, 50)
In [91]:
# seed sp.rand with random_state to obtain a reproducible result
M1 = sp.rand(N, N, density=0.2, random_state=112312321)
M2 = sp.rand(N, N, density=0.2, random_state=984592134)
In [92]:
evals = np.array([sp_eigs_sorted((1-x)*M1 + x*M2, k=25)[0] for x in x_vec])
In [93]:
fig, ax = plt.subplots(figsize=(8, 4))

for idx in range(evals.shape[1]):
    ax.plot(x_vec, np.real(evals[:,idx]), lw=0.5)

ax.set_xlabel(r"$x$", fontsize=16)
ax.set_ylabel(r"eig.vals. of $(1-x)M_1+xM_2$", fontsize=16)

fig.tight_layout()
fig.savefig("ch10-sparse-eigs.pdf")

Graphs

In [94]:
g = nx.MultiGraph()
In [95]:
g.add_node(1)
In [96]:
g.nodes()
Out[96]:
[1]
In [97]:
g.add_nodes_from([3, 4, 5])
In [98]:
g.nodes()
Out[98]:
[1, 3, 4, 5]
In [99]:
g.add_edge(1, 2)
In [100]:
g.edges()
Out[100]:
[(1, 2)]
In [101]:
g.add_edges_from([(3, 4), (5, 6)])
In [102]:
g.edges()
Out[102]:
[(1, 2), (3, 4), (5, 6)]
In [103]:
g.add_weighted_edges_from([(1, 3, 1.5), (3, 5, 2.5)])
In [104]:
g.edges(data=True)
Out[104]:
[(1, 2, {}),
 (1, 3, {'weight': 1.5}),
 (3, 4, {}),
 (3, 5, {'weight': 2.5}),
 (5, 6, {})]
In [105]:
g.add_weighted_edges_from([(6, 7, 1.5)])
In [106]:
g.nodes()
Out[106]:
[1, 2, 3, 4, 5, 6, 7]
In [107]:
g.edges()
Out[107]:
[(1, 2), (1, 3), (3, 4), (3, 5), (5, 6), (6, 7)]
In [108]:
import numpy as np
In [109]:
import json
In [110]:
with open("tokyo-metro.json") as f:
    data = json.load(f)
In [111]:
data.keys()
Out[111]:
[u'C', u'G', u'F', u'H', u'M', u'N', u'T', u'Y', u'Z']
In [112]:
data["C"]
Out[112]:
{u'color': u'#149848',
 u'transfers': [[u'C3', u'F15'],
  [u'C4', u'Z2'],
  [u'C4', u'G2'],
  [u'C7', u'M14'],
  [u'C7', u'N6'],
  [u'C7', u'G6'],
  [u'C8', u'M15'],
  [u'C8', u'H6'],
  [u'C9', u'H7'],
  [u'C9', u'Y18'],
  [u'C11', u'T9'],
  [u'C11', u'M18'],
  [u'C11', u'Z8'],
  [u'C12', u'M19'],
  [u'C18', u'H21']],
 u'travel_times': [[u'C1', u'C2', 2],
  [u'C2', u'C3', 2],
  [u'C3', u'C4', 1],
  [u'C4', u'C5', 2],
  [u'C5', u'C6', 2],
  [u'C6', u'C7', 2],
  [u'C7', u'C8', 1],
  [u'C8', u'C9', 3],
  [u'C9', u'C10', 1],
  [u'C10', u'C11', 2],
  [u'C11', u'C12', 2],
  [u'C12', u'C13', 2],
  [u'C13', u'C14', 2],
  [u'C14', u'C15', 2],
  [u'C15', u'C16', 2],
  [u'C16', u'C17', 3],
  [u'C17', u'C18', 3],
  [u'C18', u'C19', 3]]}
In [113]:
data
Out[113]:
{u'C': {u'color': u'#149848',
  u'transfers': [[u'C3', u'F15'],
   [u'C4', u'Z2'],
   [u'C4', u'G2'],
   [u'C7', u'M14'],
   [u'C7', u'N6'],
   [u'C7', u'G6'],
   [u'C8', u'M15'],
   [u'C8', u'H6'],
   [u'C9', u'H7'],
   [u'C9', u'Y18'],
   [u'C11', u'T9'],
   [u'C11', u'M18'],
   [u'C11', u'Z8'],
   [u'C12', u'M19'],
   [u'C18', u'H21']],
  u'travel_times': [[u'C1', u'C2', 2],
   [u'C2', u'C3', 2],
   [u'C3', u'C4', 1],
   [u'C4', u'C5', 2],
   [u'C5', u'C6', 2],
   [u'C6', u'C7', 2],
   [u'C7', u'C8', 1],
   [u'C8', u'C9', 3],
   [u'C9', u'C10', 1],
   [u'C10', u'C11', 2],
   [u'C11', u'C12', 2],
   [u'C12', u'C13', 2],
   [u'C13', u'C14', 2],
   [u'C14', u'C15', 2],
   [u'C15', u'C16', 2],
   [u'C16', u'C17', 3],
   [u'C17', u'C18', 3],
   [u'C18', u'C19', 3]]},
 u'F': {u'color': u'#b96528',
  u'transfers': [[u'F1', u'Y1'],
   [u'F2', u'Y2'],
   [u'F3', u'Y3'],
   [u'F4', u'Y4'],
   [u'F5', u'Y5'],
   [u'F6', u'Y6'],
   [u'F7', u'Y7'],
   [u'F8', u'Y8'],
   [u'F9', u'Y9'],
   [u'F9', u'M25'],
   [u'F13', u'M9'],
   [u'F15', u'C3'],
   [u'F16', u'Z1'],
   [u'F16', u'G1']],
  u'travel_times': [[u'F1', u'F2', 3],
   [u'F2', u'F3', 2],
   [u'F3', u'F4', 3],
   [u'F4', u'F5', 2],
   [u'F5', u'F6', 2],
   [u'F6', u'F7', 2],
   [u'F7', u'F8', 2],
   [u'F8', u'F9', 2],
   [u'F9', u'F10', 3],
   [u'F10', u'F11', 2],
   [u'F11', u'F12', 2],
   [u'F12', u'F13', 2],
   [u'F13', u'F14', 3],
   [u'F14', u'F15', 2],
   [u'F15', u'F16', 2]]},
 u'G': {u'color': u'#f59230',
  u'transfers': [[u'G1', u'Z1'],
   [u'G1', u'F16'],
   [u'G2', u'Z2'],
   [u'G2', u'C4'],
   [u'G4', u'Z3'],
   [u'G5', u'M13'],
   [u'G5', u'Y16'],
   [u'G5', u'Z4'],
   [u'G5', u'N7'],
   [u'G6', u'N6'],
   [u'G6', u'M14'],
   [u'G6', u'C7'],
   [u'G9', u'M16'],
   [u'G9', u'H8'],
   [u'G11', u'T10'],
   [u'G12', u'Z9'],
   [u'G15', u'H16'],
   [u'G16', u'H17']],
  u'travel_times': [[u'G1', u'G2', 2],
   [u'G2', u'G3', 1],
   [u'G3', u'G4', 2],
   [u'G4', u'G5', 2],
   [u'G5', u'G6', 2],
   [u'G6', u'G7', 2],
   [u'G7', u'G8', 2],
   [u'G8', u'G9', 2],
   [u'G9', u'G10', 1],
   [u'G10', u'G11', 2],
   [u'G11', u'G12', 2],
   [u'G12', u'G13', 1],
   [u'G13', u'G14', 2],
   [u'G14', u'G15', 2],
   [u'G15', u'G16', 1],
   [u'G16', u'G17', 2],
   [u'G17', u'G18', 1],
   [u'G18', u'G19', 2]]},
 u'H': {u'color': u'#9cacb5',
  u'transfers': [[u'H6', u'M15'],
   [u'H6', u'C8'],
   [u'H7', u'Y18'],
   [u'H7', u'C9'],
   [u'H8', u'M16'],
   [u'H8', u'G9'],
   [u'H12', u'T11'],
   [u'H16', u'G15'],
   [u'H17', u'G16'],
   [u'H21', u'C18']],
  u'travel_times': [[u'H1', u'H2', 3],
   [u'H2', u'H3', 3],
   [u'H3', u'H4', 3],
   [u'H4', u'H5', 3],
   [u'H5', u'H6', 2],
   [u'H6', u'H7', 3],
   [u'H7', u'H8', 1],
   [u'H8', u'H9', 2],
   [u'H9', u'H10', 2],
   [u'H10', u'H11', 2],
   [u'H11', u'H12', 1],
   [u'H12', u'H13', 3],
   [u'H13', u'H14', 1],
   [u'H14', u'H15', 2],
   [u'H15', u'H16', 2],
   [u'H16', u'H17', 1],
   [u'H17', u'H18', 2],
   [u'H18', u'H19', 2],
   [u'H19', u'H20', 2],
   [u'H20', u'H21', 3]]},
 u'M': {u'color': u'#ff0000',
  u'transfers': [[u'M9', u'F13'],
   [u'M12', u'N8'],
   [u'M13', u'G5'],
   [u'M13', u'Y16'],
   [u'M13', u'Z4'],
   [u'M13', u'N7'],
   [u'M14', u'C7'],
   [u'M14', u'G6'],
   [u'M14', u'N6'],
   [u'M15', u'H6'],
   [u'M15', u'C8'],
   [u'M16', u'G9'],
   [u'M16', u'H8'],
   [u'M18', u'T9'],
   [u'M18', u'C11'],
   [u'M18', u'Z8'],
   [u'M19', u'C12'],
   [u'M22', u'N11'],
   [u'M25', u'Y9'],
   [u'M25', u'F9']],
  u'travel_times': [[u'M1', u'M2', 2],
   [u'M2', u'M3', 2],
   [u'M3', u'M4', 2],
   [u'M4', u'M5', 2],
   [u'M5', u'M6', 2],
   [u'M6', u'M7', 2],
   [u'M7', u'M8', 2],
   [u'M8', u'M9', 2],
   [u'M9', u'M10', 1],
   [u'M10', u'M11', 2],
   [u'M11', u'M12', 2],
   [u'M12', u'M13', 3],
   [u'M13', u'M14', 2],
   [u'M14', u'M15', 1],
   [u'M15', u'M16', 3],
   [u'M16', u'M17', 2],
   [u'M17', u'M18', 2],
   [u'M18', u'M19', 2],
   [u'M19', u'M20', 1],
   [u'M20', u'M21', 2],
   [u'M21', u'M22', 2],
   [u'M22', u'M23', 3],
   [u'M23', u'M24', 2],
   [u'M24', u'M25', 3],
   [u'm3', u'm4', 2],
   [u'm4', u'm5', 2],
   [u'm5', u'M6', 2]]},
 u'N': {u'color': u'#1aaca9',
  u'transfers': [[u'N1', u'T1'],
   [u'N2', u'T2'],
   [u'N3', u'T3'],
   [u'N6', u'G6'],
   [u'N6', u'M14'],
   [u'N6', u'C7'],
   [u'N7', u'Y16'],
   [u'N7', u'Z4'],
   [u'N7', u'G5'],
   [u'N7', u'M13'],
   [u'N8', u'M12'],
   [u'N9', u'Y14'],
   [u'N10', u'Y13'],
   [u'N10', u'T6'],
   [u'N11', u'M22']],
  u'travel_times': [[u'N1', u'N2', 2],
   [u'N2', u'N3', 2],
   [u'N3', u'N4', 2],
   [u'N4', u'N5', 2],
   [u'N5', u'N6', 2],
   [u'N6', u'N7', 2],
   [u'N7', u'N8', 2],
   [u'N8', u'N9', 2],
   [u'N9', u'N10', 2],
   [u'N10', u'N11', 2],
   [u'N11', u'N12', 3],
   [u'N12', u'N13', 2],
   [u'N13', u'N14', 2],
   [u'N14', u'N15', 3],
   [u'N15', u'N16', 1],
   [u'N16', u'N17', 3],
   [u'N17', u'N18', 2],
   [u'N18', u'N19', 2]]},
 u'T': {u'color': u'#1aa7d8',
  u'transfers': [[u'T6', u'N10'],
   [u'T6', u'Y13'],
   [u'T7', u'Z6'],
   [u'T9', u'M18'],
   [u'T9', u'C11'],
   [u'T9', u'Z8'],
   [u'T10', u'G11'],
   [u'T11', u'H12']],
  u'travel_times': [[u'T1', u'T2', 0],
   [u'T2', u'T3', 3],
   [u'T3', u'T4', 6],
   [u'T4', u'T5', 9],
   [u'T5', u'T6', 11],
   [u'T6', u'T7', 13],
   [u'T7', u'T8', 14],
   [u'T8', u'T9', 16],
   [u'T9', u'T10', 18],
   [u'T10', u'T11', 20],
   [u'T11', u'T12', 21],
   [u'T12', u'T13', 24],
   [u'T13', u'T14', 26],
   [u'T14', u'T15', 27],
   [u'T15', u'T16', 30],
   [u'T16', u'T17', 33],
   [u'T17', u'T18', 35],
   [u'T18', u'T19', 37],
   [u'T19', u'T20', 39],
   [u'T20', u'T21', 41],
   [u'T21', u'T22', 43],
   [u'T22', u'T23', 46],
   [u'T23', u'T24', 49]]},
 u'Y': {u'color': u'#ede7c3',
  u'transfers': [[u'Y1', u'F1'],
   [u'Y2', u'F2'],
   [u'Y3', u'F3'],
   [u'Y4', u'F4'],
   [u'Y5', u'F5'],
   [u'Y6', u'F6'],
   [u'Y7', u'F7'],
   [u'Y8', u'F8'],
   [u'Y9', u'F9'],
   [u'Y9', u'M25'],
   [u'Y13', u'T6'],
   [u'Y13', u'N10'],
   [u'Y14', u'N9'],
   [u'Y16', u'Z4'],
   [u'Y16', u'N7'],
   [u'Y16', u'G5'],
   [u'Y16', u'M13'],
   [u'Y18', u'H7'],
   [u'Y18', u'C9']],
  u'travel_times': [[u'Y1', u'Y2', 4],
   [u'Y2', u'Y3', 2],
   [u'Y3', u'Y4', 3],
   [u'Y4', u'Y5', 2],
   [u'Y5', u'Y6', 2],
   [u'Y6', u'Y7', 2],
   [u'Y7', u'Y8', 2],
   [u'Y8', u'Y9', 3],
   [u'Y9', u'Y10', 2],
   [u'Y10', u'Y11', 2],
   [u'Y11', u'Y12', 2],
   [u'Y12', u'Y13', 3],
   [u'Y13', u'Y14', 2],
   [u'Y14', u'Y15', 2],
   [u'Y15', u'Y16', 1],
   [u'Y16', u'Y17', 2],
   [u'Y17', u'Y18', 2],
   [u'Y18', u'Y19', 2],
   [u'Y19', u'Y20', 2],
   [u'Y20', u'Y21', 2],
   [u'Y21', u'Y22', 2],
   [u'Y22', u'Y23', 3],
   [u'Y23', u'Y24', 2]]},
 u'Z': {u'color': u'#a384bf',
  u'transfers': [[u'Z1', u'F16'],
   [u'Z1', u'G1'],
   [u'Z2', u'C4'],
   [u'Z2', u'G2'],
   [u'Z3', u'G4'],
   [u'Z4', u'Y16'],
   [u'Z4', u'N7'],
   [u'Z4', u'M13'],
   [u'Z4', u'G5'],
   [u'Z6', u'T7'],
   [u'Z8', u'M18'],
   [u'Z8', u'C11'],
   [u'Z8', u'T9'],
   [u'Z9', u'G12']],
  u'travel_times': [[u'Z1', u'Z2', 3],
   [u'Z2', u'Z3', 2],
   [u'Z3', u'Z4', 2],
   [u'Z4', u'Z5', 2],
   [u'Z5', u'Z6', 2],
   [u'Z6', u'Z7', 2],
   [u'Z7', u'Z8', 2],
   [u'Z8', u'Z9', 2],
   [u'Z9', u'Z10', 3],
   [u'Z10', u'Z11', 3],
   [u'Z11', u'Z12', 3],
   [u'Z12', u'Z13', 2],
   [u'Z13', u'Z14', 2]]}}
In [114]:
g = nx.Graph()

for line in data.values():
    g.add_weighted_edges_from(line["travel_times"])
    g.add_edges_from(line["transfers"])
In [115]:
for n1, n2 in g.edges_iter():
    g[n1][n2]["transfer"] = "weight" not in g[n1][n2]
In [116]:
g.number_of_nodes()
Out[116]:
184
In [117]:
g.nodes()[:5]
Out[117]:
[u'Z13', u'N1', u'Y20', u'Y21', u'Y22']
In [118]:
g.number_of_edges()
Out[118]:
243
In [119]:
g.edges()[:5]
Out[119]:
[(u'Z13', u'Z14'),
 (u'Z13', u'Z12'),
 (u'N1', u'N2'),
 (u'N1', u'T1'),
 (u'Y20', u'Y21')]
In [120]:
on_foot = [edge for edge in g.edges_iter() if g.get_edge_data(*edge)["transfer"]]
In [121]:
on_train = [edge for edge in g.edges_iter() if not g.get_edge_data(*edge)["transfer"]]
In [122]:
colors = [data[n[0].upper()]["color"] for n in g.nodes()]
In [123]:
fig, ax = plt.subplots(1, 1, figsize=(14, 10))

pos = nx.graphviz_layout(g, prog="neato")
nx.draw(g, pos, ax=ax, node_size=200, node_color=colors)
nx.draw_networkx_labels(g, pos=pos, ax=ax, font_size=6)
nx.draw_networkx_edges(g, pos=pos, ax=ax, edgelist=on_train, width=2)
nx.draw_networkx_edges(g, pos=pos, ax=ax, edgelist=on_foot, edge_color="blue")

# removing the default axis on all sides:
for side in ['bottom','right','top','left']:
    ax.spines[side].set_visible(False)

# removing the axis labels and ticks
ax.set_xticks([])
ax.set_yticks([])
ax.xaxis.set_ticks_position('none')
ax.yaxis.set_ticks_position('none')

fig.savefig("ch10-metro-graph.pdf")
fig.savefig("ch10-metro-graph.png")
fig.tight_layout()
In [124]:
#g.degree()
In [125]:
d_max = max(g.degree().values())
In [126]:
[(n, d) for (n, d) in g.degree().items() if d == d_max]
Out[126]:
[(u'G5', 6), (u'Y16', 6), (u'Z4', 6), (u'N7', 6), (u'M13', 6)]
In [127]:
p = nx.shortest_path(g, "Y24", "C19")
In [128]:
np.array(p)
Out[128]:
array([u'Y24', u'Y23', u'Y22', u'Y21', u'Y20', u'Y19', u'Y18', u'C9',
       u'C10', u'C11', u'C12', u'C13', u'C14', u'C15', u'C16', u'C17',
       u'C18', u'C19'], 
      dtype='<U3')
In [129]:
np.sum([g[p[n]][p[n+1]]["weight"] for n in range(len(p)-1) if "weight" in g[p[n]][p[n+1]]])
Out[129]:
35
In [130]:
h = g.copy()
In [131]:
for n1, n2 in h.edges_iter():
    if "transfer" in h[n1][n2]:
        h[n1][n2]["weight"] = 5
In [132]:
p = nx.shortest_path(h, "Y24", "C19")
In [133]:
np.array(p)
Out[133]:
array([u'Y24', u'Y23', u'Y22', u'Y21', u'Y20', u'Y19', u'Y18', u'C9',
       u'C10', u'C11', u'C12', u'C13', u'C14', u'C15', u'C16', u'C17',
       u'C18', u'C19'], 
      dtype='<U3')
In [134]:
np.sum([h[p[n]][p[n+1]]["weight"] for n in range(len(p)-1)])
Out[134]:
85
In [135]:
p = nx.shortest_path(h, "Z1", "H16")
In [136]:
np.sum([h[p[n]][p[n+1]]["weight"] for n in range(len(p)-1)])
Out[136]:
65
In [137]:
A = nx.to_scipy_sparse_matrix(g)
In [138]:
A
Out[138]:
<184x184 sparse matrix of type '<type 'numpy.int64'>'
	with 486 stored elements in Compressed Sparse Row format>
In [139]:
perm = sp.csgraph.reverse_cuthill_mckee(A)
In [140]:
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4))
ax1.spy(A, markersize=2)
ax2.spy(sp_permute(A, perm, perm), markersize=2)
fig.tight_layout()
fig.savefig("ch12-rcm-graph.pdf")

Versions

In [141]:
%reload_ext version_information
In [142]:
%version_information numpy, scipy, matplotlib, networkx, pygraphviz
Out[142]:
SoftwareVersion
Python2.7.10 64bit [GCC 4.2.1 (Apple Inc. build 5577)]
IPython3.2.1
OSDarwin 14.1.0 x86_64 i386 64bit
numpy1.9.2
scipy0.16.0
matplotlib1.4.3
networkx1.9.1
pygraphviz1.2