Mercurial > public > think_complexity

annotate SmallWorldGraph.py @ 24:5c2c4ce095ef

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
A stab at the L(p)/L(0) plot. I still don't quite get how the graphs in the Watts and Strogatz paper were generated. My results have basically the same shape, but don't converge to 0. I'm not sure how this is possible if the rewire function does not remove edges.
author Brian Neal <bgneal@gmail.com>
date Thu, 03 Jan 2013 18:41:13 -0600
parents 74c9d126bd05
children a46783561538
rev   line source
bgneal@23 1 """This module contains the SmallWorldGraph class for 4.4, exercise 4.
bgneal@23 2
bgneal@23 3 """
bgneal@23 4 import itertools
bgneal@23 5 import random
bgneal@23 6
bgneal@23 7 from Graph import Edge
bgneal@23 8 from RandomGraph import RandomGraph
bgneal@23 9
bgneal@23 10
bgneal@24 11 INFINITY = float('Inf')
bgneal@24 12
bgneal@23 13 class SmallWorldGraph(RandomGraph):
bgneal@23 14 """A small world graph for 4.4, exercise 4 in Think Complexity."""
bgneal@23 15
bgneal@23 16 def __init__(self, vs, k, p):
bgneal@23 17 """Create a small world graph. Parameters:
bgneal@23 18 vs - a list of vertices
bgneal@23 19 k - regular graph parameter k: the number of edges between vertices
bgneal@23 20 p - probability for rewiring.
bgneal@23 21
bgneal@23 22 First a regular graph is created from the list of vertices and parameter
bgneal@23 23 k, which specifies how many edges each vertex should have. Then the
bgneal@23 24 rewire() method is called with parameter p to rewire the graph according
bgneal@23 25 to the algorithm by Watts and Strogatz.
bgneal@23 26
bgneal@23 27 """
bgneal@23 28 super(SmallWorldGraph, self).__init__(vs=vs)
bgneal@23 29 self.add_regular_edges(k)
bgneal@23 30 self.rewire(p)
bgneal@23 31
bgneal@24 32 # give each edge a default length of 1; this is used by the
bgneal@24 33 # shortest_path_tree() method
bgneal@24 34 for e in self.edges():
bgneal@24 35 e.length = 1
bgneal@24 36
bgneal@23 37 def rewire(self, p):
bgneal@23 38 """Assuming the graph is initially a regular graph, rewires the graph
bgneal@23 39 using the Watts and Strogatz algorithm.
bgneal@23 40
bgneal@23 41 """
bgneal@23 42 # We iterate over the edges in a random order:
bgneal@23 43 edges = self.edges()
bgneal@23 44 random.shuffle(edges)
bgneal@23 45 vertices = self.vertices()
bgneal@23 46
bgneal@23 47 for e in edges:
bgneal@23 48 if random.random() < p:
bgneal@23 49 v, w = e # remember the 2 vertices this edge connected
bgneal@23 50 self.remove_edge(e) # remove from graph
bgneal@23 51
bgneal@23 52 # pick another vertex to connect to v; duplicate edges are
bgneal@23 53 # forbidden
bgneal@23 54 while True:
bgneal@23 55 x = random.choice(vertices)
bgneal@23 56 if v is not x and not self.get_edge(v, x):
bgneal@23 57 self.add_edge(Edge(v, x))
bgneal@23 58 break
bgneal@23 59
bgneal@23 60 def clustering_coefficient(self):
bgneal@23 61 """Compute the clustering coefficient for this graph as defined by Watts
bgneal@23 62 and Strogatz.
bgneal@23 63
bgneal@23 64 """
bgneal@23 65 cv = {}
bgneal@23 66 for v in self:
bgneal@23 67 # consider a node and its neighbors
bgneal@23 68 nodes = self.out_vertices(v)
bgneal@23 69 nodes.append(v)
bgneal@23 70
bgneal@23 71 # compute the maximum number of possible edges between these nodes
bgneal@23 72 # if they were all connected to each other:
bgneal@23 73 n = len(nodes)
bgneal@23 74 if n == 1:
bgneal@23 75 # edge case of only 1 node; handle this case to avoid division
bgneal@23 76 # by zero in the general case
bgneal@23 77 cv[v] = 1.0
bgneal@23 78 continue
bgneal@23 79
bgneal@23 80 possible = n * (n - 1) / 2.0
bgneal@23 81
bgneal@23 82 # now compute how many edges actually exist between the nodes
bgneal@23 83 actual = 0
bgneal@23 84 for x, y in itertools.combinations(nodes, 2):
bgneal@23 85 if self.get_edge(x, y):
bgneal@23 86 actual += 1
bgneal@23 87
bgneal@23 88 # the fraction of actual / possible is this nodes C sub v value
bgneal@23 89 cv[v] = actual / possible
bgneal@23 90
bgneal@23 91 # The clustering coefficient is the average of all C sub v values
bgneal@23 92 if len(cv):
bgneal@23 93 return sum(cv.values()) / float(len(cv))
bgneal@23 94 return 0.0
bgneal@24 95
bgneal@24 96 def shortest_path_tree(self, source, hint=None):
bgneal@24 97 """Finds the length of the shortest path from the source vertex to all
bgneal@24 98 other vertices in the graph. This length is stored on the vertices as an
bgneal@24 99 attribute named 'dist'. The algorithm used is Dijkstra's.
bgneal@24 100
bgneal@24 101 hint: if provided, must be a dictionary mapping tuples to already known
bgneal@24 102 shortest path distances. This can be used to speed up the algorithm.
bgneal@24 103
bgneal@24 104 """
bgneal@24 105 if not hint:
bgneal@24 106 hint = {}
bgneal@24 107
bgneal@24 108 for v in self.vertices():
bgneal@24 109 v.dist = hint.get((source, v), INFINITY)
bgneal@24 110 source.dist = 0
bgneal@24 111
bgneal@24 112 queue = [v for v in self.vertices() if v.dist < INFINITY]
bgneal@24 113 sort_flag = True
bgneal@24 114 while len(queue):
bgneal@24 115
bgneal@24 116 if sort_flag:
bgneal@24 117 queue.sort(key=lambda v: v.dist)
bgneal@24 118 sort_flag = False
bgneal@24 119
bgneal@24 120 v = queue.pop(0)
bgneal@24 121
bgneal@24 122 # for each neighbor of v, see if we found a new shortest path
bgneal@24 123 for w, e in self[v].iteritems():
bgneal@24 124 d = v.dist + e.length
bgneal@24 125 if d < w.dist:
bgneal@24 126 w.dist = d
bgneal@24 127 queue.append(w)
bgneal@24 128 sort_flag = True
bgneal@24 129
bgneal@24 130 def big_l(self):
bgneal@24 131 """Computes the "big-L" value for the graph as per Watts & Strogatz.
bgneal@24 132
bgneal@24 133 big_l() is defined as the number of edges in the shortest path between
bgneal@24 134 two vertices, averaged over all vertices.
bgneal@24 135
bgneal@24 136 """
bgneal@24 137 d = {}
bgneal@24 138 for v in self.vertices():
bgneal@24 139 self.shortest_path_tree(v, d)
bgneal@24 140 t = [((w, v), w.dist) for w in self.vertices() if v is not w]
bgneal@24 141 d.update(t)
bgneal@24 142
bgneal@24 143 if len(d):
bgneal@24 144 return sum(d.values()) / float(len(d))
bgneal@24 145 return 0.0