bgneal@1: """ Code example from Complexity and Computation, a book about
bgneal@1: exploring complexity science with Python. Available free from
bgneal@1:
bgneal@1: http://greenteapress.com/complexity
bgneal@1:
bgneal@1: Copyright 2011 Allen B. Downey.
bgneal@1: Distributed under the GNU General Public License at gnu.org/licenses/gpl.html.
bgneal@1: """
bgneal@36: import heapq
bgneal@1: import itertools
bgneal@35: from collections import deque, Counter
bgneal@1:
bgneal@1:
bgneal@36: INFINITY = float('Inf')
bgneal@36:
bgneal@36:
bgneal@5: class GraphError(Exception):
bgneal@5: """Exception for Graph errors"""
bgneal@5: pass
bgneal@5:
bgneal@5:
bgneal@1: class Vertex(object):
bgneal@1: """A Vertex is a node in a graph."""
bgneal@1:
bgneal@1: def __init__(self, label=''):
bgneal@1: self.label = label
bgneal@1:
bgneal@1: def __repr__(self):
bgneal@1: """Returns a string representation of this object that can
bgneal@1: be evaluated as a Python expression."""
bgneal@1: return 'Vertex(%s)' % repr(self.label)
bgneal@1:
bgneal@1: __str__ = __repr__
bgneal@1: """The str and repr forms of this object are the same."""
bgneal@1:
bgneal@1:
bgneal@1: class Edge(tuple):
bgneal@1: """An Edge is a list of two vertices."""
bgneal@1:
bgneal@1: def __new__(cls, *vs):
bgneal@1: """The Edge constructor takes two vertices."""
bgneal@1: if len(vs) != 2:
bgneal@1: raise ValueError, 'Edges must connect exactly two vertices.'
bgneal@1: return tuple.__new__(cls, vs)
bgneal@1:
bgneal@1: def __repr__(self):
bgneal@1: """Return a string representation of this object that can
bgneal@1: be evaluated as a Python expression."""
bgneal@1: return 'Edge(%s, %s)' % (repr(self[0]), repr(self[1]))
bgneal@1:
bgneal@1: __str__ = __repr__
bgneal@1: """The str and repr forms of this object are the same."""
bgneal@1:
bgneal@1:
bgneal@1: class Graph(dict):
bgneal@1: """A Graph is a dictionary of dictionaries. The outer
bgneal@1: dictionary maps from a vertex to an inner dictionary.
bgneal@1: The inner dictionary maps from other vertices to edges.
bgneal@1:
bgneal@1: For vertices a and b, graph[a][b] maps
bgneal@1: to the edge that connects a->b, if it exists."""
bgneal@1:
bgneal@6: def __init__(self, vs=None, es=None):
bgneal@1: """Creates a new graph.
bgneal@1: vs: list of vertices;
bgneal@1: es: list of edges.
bgneal@1: """
bgneal@6: if vs:
bgneal@6: for v in vs:
bgneal@6: self.add_vertex(v)
bgneal@1:
bgneal@6: if es:
bgneal@6: for e in es:
bgneal@6: self.add_edge(e)
bgneal@1:
bgneal@36: def set_edge_length(self, n=1):
bgneal@36: """Give each edge a length of n; this is used by the
bgneal@36: shortest_path_tree() method.
bgneal@36:
bgneal@36: """
bgneal@36: for e in self.edges():
bgneal@36: e.length = n
bgneal@36:
bgneal@1: def add_vertex(self, v):
bgneal@1: """Add a vertex to the graph."""
bgneal@1: self[v] = {}
bgneal@1:
bgneal@1: def add_edge(self, e):
bgneal@1: """Adds and edge to the graph by adding an entry in both directions.
bgneal@1:
bgneal@1: If there is already an edge connecting these Vertices, the
bgneal@1: new edge replaces it.
bgneal@1: """
bgneal@1: v, w = e
bgneal@1: self[v][w] = e
bgneal@1: self[w][v] = e
bgneal@1:
bgneal@1: def get_edge(self, v, w):
bgneal@1: """Returns the edge object that exists between the two vertices v & w,
bgneal@1: or None if no such edge exists.
bgneal@1:
bgneal@1: """
bgneal@1: try:
bgneal@1: return self[v][w]
bgneal@1: except KeyError:
bgneal@1: return None
bgneal@1:
bgneal@1: def remove_edge(self, e):
bgneal@1: """Removes the edge e from the graph."""
bgneal@1:
bgneal@4: v, w = e
bgneal@4: del self[v][w]
bgneal@4: del self[w][v]
bgneal@1:
bgneal@1: def vertices(self):
bgneal@1: """Returns a list of the vertices in the graph."""
bgneal@1:
bgneal@1: return self.keys()
bgneal@1:
bgneal@1: def edges(self):
bgneal@1: """"Returns a list of the edges in the graph."""
bgneal@1:
bgneal@1: edge_set = set()
bgneal@4: for d in self.itervalues():
bgneal@4: edge_set.update(d.itervalues())
bgneal@1:
bgneal@1: return list(edge_set)
bgneal@1:
bgneal@1: def out_vertices(self, v):
bgneal@1: """Returns a list of vertices that are adjacent to the given vertex v.
bgneal@1:
bgneal@1: """
bgneal@1: return self[v].keys()
bgneal@1:
bgneal@1: def out_edges(self, v):
bgneal@1: """Returns a list of edges connected to a given vertex v."""
bgneal@1:
bgneal@1: return self[v].values()
bgneal@1:
bgneal@5: def remove_all_edges(self):
bgneal@5: """Removes all edges in the graph."""
bgneal@5: for v in self.iterkeys():
bgneal@5: self[v] = {}
bgneal@5:
bgneal@1: def add_all_edges(self):
bgneal@1: """Makes the graph complete by adding edges between all pairs of
bgneal@1: vertices.
bgneal@1:
bgneal@1: """
bgneal@1: # Clear all edges first
bgneal@5: self.remove_all_edges()
bgneal@1:
bgneal@5: # For each combination of 2 vertices, create an edge between them:
bgneal@1: for v, w in itertools.combinations(self.iterkeys(), 2):
bgneal@6: self.add_edge(Edge(v, w))
bgneal@1:
bgneal@5: def add_regular_edges(self, k):
bgneal@5: """Makes the graph regular by making every vertex have k edges.
bgneal@5:
bgneal@5: It is not always possible to create a regular graph with a given degree.
bgneal@5: If a graph has n vertices, then a regular graph can be constructed with
bgneal@5: degree k if n >= k + 1 and n * k is even. If these conditions are not
bgneal@5: met a GraphError exception is raised.
bgneal@5:
bgneal@5: """
bgneal@5: n = len(self.vertices())
bgneal@5: if n < k + 1:
bgneal@5: raise GraphError("Can't make a regular graph with degree >= number"
bgneal@5: " of vertices")
bgneal@5: if (n * k) % 2 != 0:
bgneal@5: raise GraphError("Can't make a regular graph of degree k and"
bgneal@5: " order n where k * n is odd")
bgneal@5:
bgneal@5: # Remove all edges first
bgneal@5: self.remove_all_edges()
bgneal@5:
bgneal@5: if k % 2 != 0: # if k is odd
bgneal@5: self._add_regular_edges_even(k - 1)
bgneal@5: self._add_regular_edges_odd()
bgneal@5: else:
bgneal@5: self._add_regular_edges_even(k)
bgneal@5:
bgneal@5: def _add_regular_edges_even(self, k):
bgneal@5: """Make a regular graph with degree k. k must be even."""
bgneal@5:
bgneal@5: vs = self.vertices()
bgneal@5: vs2 = vs * 2
bgneal@5:
bgneal@5: for i, v in enumerate(vs):
bgneal@5: for j in range(1, k / 2 + 1):
bgneal@5: w = vs2[i + j]
bgneal@5: self.add_edge(Edge(v, w))
bgneal@5:
bgneal@5: def _add_regular_edges_odd(self):
bgneal@5: """Adds an extra edge across the graph to finish off a regular graph
bgneal@5: with odd degree. The number of vertices must be even.
bgneal@5:
bgneal@5: """
bgneal@5: vs = self.vertices()
bgneal@5: vs2 = vs * 2
bgneal@5: n = len(vs)
bgneal@5:
bgneal@5: for i in range(n / 2):
bgneal@5: v = vs2[i]
bgneal@5: w = vs2[i + n / 2]
bgneal@5: self.add_edge(Edge(v, w))
bgneal@5:
bgneal@7: def bfs(self, start, visit_func=None):
bgneal@7: """Perform a breadth first search starting at node start.
bgneal@7:
bgneal@7: The function visit_func, if supplied, is invoked on each node.
bgneal@7:
bgneal@8: The set of visited nodes is returned.
bgneal@8:
bgneal@7: """
bgneal@8: visited = set()
bgneal@7:
bgneal@7: # Create a work queue consisting initially of the starting node
bgneal@8: queue = deque([start])
bgneal@7:
bgneal@7: while queue:
bgneal@7: # retrieve first item from the queue
bgneal@8: v = queue.popleft()
bgneal@7:
bgneal@8: if v in visited:
bgneal@7: continue # Skip this one if we've seen it before
bgneal@7:
bgneal@7: # Mark it as visited and invoke user's function on it
bgneal@8: visited.add(v)
bgneal@7: if visit_func:
bgneal@7: visit_func(v)
bgneal@7:
bgneal@7: # Add the adjacent neigbors to the node to the queue
bgneal@22: queue.extend(c for c in self.out_vertices(v) if c not in visited)
bgneal@7:
bgneal@8: return visited
bgneal@8:
bgneal@7: def is_connected(self):
bgneal@7: """Returns True if the graph is connected (there is a path from every
bgneal@7: node to every other node) and False otherwise.
bgneal@7:
bgneal@7: """
bgneal@7: vs = self.vertices()
bgneal@7: if len(vs):
bgneal@8: visited = self.bfs(vs[0])
bgneal@7: # See if all nodes have been visited
bgneal@8: return len(vs) == len(visited)
bgneal@7:
bgneal@7: return False # Graph is empty
bgneal@7:
bgneal@35: def get_p(self):
bgneal@35: """This method returns a dictionary of probabilities where each key is
bgneal@35: the connectivity k and the value is the probability [0-1] for this
bgneal@35: graph.
bgneal@35:
bgneal@35: """
bgneal@35: # First, for each vertex, count up how many neighbors it has
bgneal@35: vs = self.vertices()
bgneal@35:
bgneal@35: c = Counter()
bgneal@35: for v in vs:
bgneal@35: n = len(self.out_vertices(v))
bgneal@35: c[n] += 1
bgneal@35:
bgneal@35: n = len(vs)
bgneal@35: if n > 0:
bgneal@35: for k in c:
bgneal@35: c[k] = float(c[k]) / n
bgneal@35:
bgneal@35: return c
bgneal@35:
bgneal@36: def clustering_coefficient(self):
bgneal@36: """Compute the clustering coefficient for this graph as defined by Watts
bgneal@36: and Strogatz.
bgneal@36:
bgneal@36: """
bgneal@36: cv = {}
bgneal@36: for v in self:
bgneal@36: # consider a node and its neighbors
bgneal@36: nodes = self.out_vertices(v)
bgneal@36: nodes.append(v)
bgneal@36:
bgneal@36: # compute the maximum number of possible edges between these nodes
bgneal@36: # if they were all connected to each other:
bgneal@36: n = len(nodes)
bgneal@36: if n == 1:
bgneal@36: # edge case of only 1 node; handle this case to avoid division
bgneal@36: # by zero in the general case
bgneal@36: cv[v] = 1.0
bgneal@36: continue
bgneal@36:
bgneal@36: possible = n * (n - 1) / 2.0
bgneal@36:
bgneal@36: # now compute how many edges actually exist between the nodes
bgneal@36: actual = 0
bgneal@36: for x, y in itertools.combinations(nodes, 2):
bgneal@36: if self.get_edge(x, y):
bgneal@36: actual += 1
bgneal@36:
bgneal@36: # the fraction of actual / possible is this nodes C sub v value
bgneal@36: cv[v] = actual / possible
bgneal@36:
bgneal@36: # The clustering coefficient is the average of all C sub v values
bgneal@36: if len(cv):
bgneal@36: return sum(cv.values()) / float(len(cv))
bgneal@36: return 0.0
bgneal@36:
bgneal@36: def shortest_path_tree(self, source, hint=None):
bgneal@36: """Finds the length of the shortest path from the source vertex to all
bgneal@36: other vertices in the graph. This length is stored on the vertices as an
bgneal@36: attribute named 'dist'. The algorithm used is Dijkstra's.
bgneal@36:
bgneal@36: hint: if provided, must be a dictionary mapping tuples to already known
bgneal@36: shortest path distances. This can be used to speed up the algorithm.
bgneal@36:
bgneal@36: """
bgneal@36: if not hint:
bgneal@36: hint = {}
bgneal@36:
bgneal@36: for v in self.vertices():
bgneal@36: v.dist = hint.get((source, v), INFINITY)
bgneal@36: source.dist = 0
bgneal@36:
bgneal@36: queue = [v for v in self.vertices() if v.dist < INFINITY]
bgneal@36: sort_flag = True
bgneal@36: while len(queue):
bgneal@36:
bgneal@36: if sort_flag:
bgneal@36: queue.sort(key=lambda v: v.dist)
bgneal@36: sort_flag = False
bgneal@36:
bgneal@36: v = queue.pop(0)
bgneal@36:
bgneal@36: # for each neighbor of v, see if we found a new shortest path
bgneal@36: for w, e in self[v].iteritems():
bgneal@36: d = v.dist + e.length
bgneal@36: if d < w.dist:
bgneal@36: w.dist = d
bgneal@36: queue.append(w)
bgneal@36: sort_flag = True
bgneal@36:
bgneal@36: def shortest_path_tree2(self, source):
bgneal@36: """Finds the length of the shortest path from the source vertex to all
bgneal@36: other vertices in the graph. This length is stored on the vertices as an
bgneal@36: attribute named 'dist'. The algorithm used is Dijkstra's with a Heap
bgneal@36: used to sort/store pending nodes to be examined.
bgneal@36:
bgneal@36: """
bgneal@36: for v in self.vertices():
bgneal@36: v.dist = INFINITY
bgneal@36: source.dist = 0
bgneal@36:
bgneal@36: queue = []
bgneal@36: heapq.heappush(queue, (0, source))
bgneal@36: while len(queue):
bgneal@36:
bgneal@36: _, v = heapq.heappop(queue)
bgneal@36:
bgneal@36: # for each neighbor of v, see if we found a new shortest path
bgneal@36: for w, e in self[v].iteritems():
bgneal@36: d = v.dist + e.length
bgneal@36: if d < w.dist:
bgneal@36: w.dist = d
bgneal@36: heapq.heappush(queue, (d, w))
bgneal@36:
bgneal@36: def all_pairs_floyd_warshall(self):
bgneal@36: """Finds the shortest paths between all pairs of vertices using the
bgneal@36: Floyd-Warshall algorithm.
bgneal@36:
bgneal@36: http://en.wikipedia.org/wiki/Floyd-Warshall_algorithm
bgneal@36:
bgneal@36: """
bgneal@36: vertices = self.vertices()
bgneal@36: dist = {}
bgneal@36: for i in vertices:
bgneal@36: for j in vertices:
bgneal@36: if i is j:
bgneal@36: dist[i, j] = 0.0
bgneal@36: else:
bgneal@36: e = self.get_edge(i, j)
bgneal@36: dist[i, j] = e.length if e else INFINITY
bgneal@36:
bgneal@36: for k in vertices:
bgneal@36: for i in vertices:
bgneal@36: for j in vertices:
bgneal@36: d_ik = dist[i, k]
bgneal@36: d_kj = dist[k, j]
bgneal@36: new_cost = d_ik + d_kj
bgneal@36: if new_cost < dist[i, j]:
bgneal@36: dist[i, j] = new_cost
bgneal@36:
bgneal@36: return dist
bgneal@36:
bgneal@36: def big_l(self):
bgneal@36: """Computes the "big-L" value for the graph as per Watts & Strogatz.
bgneal@36:
bgneal@36: L is defined as the number of edges in the shortest path between
bgneal@36: two vertices, averaged over all vertices.
bgneal@36:
bgneal@36: Uses the shortest_path_tree() method, called once for every node.
bgneal@36:
bgneal@36: """
bgneal@36: d = {}
bgneal@36: for v in self.vertices():
bgneal@36: self.shortest_path_tree(v, d)
bgneal@36: t = [((w, v), w.dist) for w in self.vertices() if v is not w]
bgneal@36: d.update(t)
bgneal@36:
bgneal@36: if len(d):
bgneal@36: return sum(d.values()) / float(len(d))
bgneal@36: return 0.0
bgneal@36:
bgneal@36: def big_l2(self):
bgneal@36: """Computes the "big-L" value for the graph as per Watts & Strogatz.
bgneal@36:
bgneal@36: L is defined as the number of edges in the shortest path between
bgneal@36: two vertices, averaged over all vertices.
bgneal@36:
bgneal@36: Uses the all_pairs_floyd_warshall() method.
bgneal@36:
bgneal@36: """
bgneal@36: dist = self.all_pairs_floyd_warshall()
bgneal@36: vertices = self.vertices()
bgneal@36: result = [dist[i, j] for i in vertices for j in vertices if i is not j]
bgneal@36:
bgneal@36: if len(result):
bgneal@36: return sum(result) / float(len(result))
bgneal@36: return 0.0
bgneal@36:
bgneal@36: def big_l3(self):
bgneal@36: """Computes the "big-L" value for the graph as per Watts & Strogatz.
bgneal@36:
bgneal@36: L is defined as the number of edges in the shortest path between
bgneal@36: two vertices, averaged over all vertices.
bgneal@36:
bgneal@36: Uses the shortest_path_tree2() method, called once for every node.
bgneal@36:
bgneal@36: """
bgneal@36: d = {}
bgneal@36: for v in self.vertices():
bgneal@36: self.shortest_path_tree2(v)
bgneal@36: t = [((v, w), w.dist) for w in self.vertices() if v is not w]
bgneal@36: d.update(t)
bgneal@36:
bgneal@36: if len(d):
bgneal@36: return sum(d.values()) / float(len(d))
bgneal@36: return 0.0
bgneal@36:
bgneal@1:
bgneal@1: def main(script, *args):
bgneal@1: import pprint
bgneal@1:
bgneal@1: v = Vertex('v')
bgneal@1: print v
bgneal@1: w = Vertex('w')
bgneal@1: print w
bgneal@1: e = Edge(v, w)
bgneal@1: print e
bgneal@1: g = Graph([v,w], [e])
bgneal@1: pprint.pprint(g)
bgneal@1:
bgneal@1: print "g.get_edge(v, w): ", g.get_edge(v, w)
bgneal@1: x = Vertex('x')
bgneal@1: print "g.get_edge(v, x): ", g.get_edge(v, x)
bgneal@1:
bgneal@1: g.remove_edge(e)
bgneal@1: pprint.pprint(g)
bgneal@1:
bgneal@1: print "vertices: ", g.vertices()
bgneal@1: print "edges: ", g.edges()
bgneal@1:
bgneal@1: g.add_edge(e)
bgneal@1: u = Vertex('u')
bgneal@1: e1 = Edge(u, v)
bgneal@1: e2 = Edge(u, w)
bgneal@1: g.add_vertex(u)
bgneal@1: g.add_edge(e1)
bgneal@1: g.add_edge(e2)
bgneal@1: print "Adding vertex u and edges:"
bgneal@1: pprint.pprint(g)
bgneal@1: print "vertices: ", g.vertices()
bgneal@1: print "edges: ", g.edges()
bgneal@1:
bgneal@1: print "Out vertices for v: ", g.out_vertices(v)
bgneal@1: print "Out edges for v: ", g.out_edges(v)
bgneal@1:
bgneal@1: x = Vertex('x')
bgneal@1: g.add_vertex(x)
bgneal@1: g.add_all_edges()
bgneal@1: pprint.pprint(g)
bgneal@1:
bgneal@7: print "g is connected?", g.is_connected()
bgneal@7: edges = g.out_edges(v)
bgneal@7: for e in edges:
bgneal@7: g.remove_edge(e)
bgneal@7: pprint.pprint(g)
bgneal@7: print "g is connected?", g.is_connected()
bgneal@7:
bgneal@7: # Isolate v and check is_connected() again
bgneal@7:
bgneal@1:
bgneal@1: if __name__ == '__main__':
bgneal@1: import sys
bgneal@1: main(*sys.argv)