first commit

README.md

@@ -0,0 +1,51 @@
+Finding all the elementary circuits of a directed graph
+-------------------------------------------------------
+
+Algorithm by D. B. Johnson
+
+    [1] Finding all the elementary circuits of a directed graph.
+    D. B. Johnson, SIAM Journal on Computing 4, no. 1, 77-84, 1975.
+    http://dx.doi.org/10.1137/0204007
+
+Using the networkx package and a modified version of its simple_cycles()
+function. The algorithm was adapted so that it would not arbitrarily order
+vertices. Three lines now contain a sorted() statement.
+
+
+Usage
+-----
+
+    python cycles.py 4 0,1 0,2 1,0 1,3 2,0 3,0 3,1 3,2
+
+First argument is the number of vertices. Subsequent arguments are ordered
+pairs of comma separated vertices that make up the directed edges of the
+graph.
+
+DOT file input
+--------------
+
+For simplicity, there is no DOT file parser included but the following allows
+to create a suitable argument string for simple DOT graphs.
+
+Given a DOT file of a simple (no labels, colors, styles, only pairs of
+vertices...) directed graph, the following line produces commandline
+arguments in the above format for that graph.
+
+	echo `sed -n -e '/^\s*[0-9]\+;$/p' graph.dot | wc -l` `sed -n -e 's/^\s*\([0-9]\) -> \([0-9]\);$/\1,\2/p' graph.dot`
+
+The above line works on DOT files like the following:
+
+    digraph G {
+      0;
+      1;
+      2;
+      0 -> 1;
+      0 -> 2;
+      1 -> 0;
+      2 -> 0;
+      2 -> 1;
+      }
+
+It would produce the following output:
+
+    3 0,1 0,2 1,0 2,0 2,1

cycles.py

@@ -0,0 +1,81 @@
+# Copyright (C) 2013 by Johannes Schauer <j.schauer@email.de>
+#
+# Copyright (C) 2010 by 
+# Aric Hagberg <hagberg@lanl.gov>
+# Dan Schult <dschult@colgate.edu>
+# Pieter Swart <swart@lanl.gov>
+# All rights reserved.
+# BSD license.
+
+
+import sys
+import networkx as nx
+from collections import defaultdict
+
+if len(sys.argv) < 3:
+    print "usage: %s num_vertices [v1,v2...]"%(sys.argv[0])
+
+def simple_cycles(G):
+    def _unblock(thisnode):
+        """Recursively unblock and remove nodes from B[thisnode]."""
+        if blocked[thisnode]:
+            blocked[thisnode] = False
+            while B[thisnode]:
+                _unblock(B[thisnode].pop())
+
+    def circuit(thisnode, startnode, component):
+        closed = False # set to True if elementary path is closed
+        path.append(thisnode)
+        blocked[thisnode] = True
+        for nextnode in sorted(component[thisnode]): # direct successors of thisnode
+            if nextnode == startnode:
+                result.append(path + [startnode])
+                closed = True
+            elif not blocked[nextnode]:
+                if circuit(nextnode, startnode, component):
+                    closed = True
+        if closed:
+            _unblock(thisnode)
+        else:
+            for nextnode in component[thisnode]:
+                if thisnode not in B[nextnode]: # TODO: use set for speedup?
+                    B[nextnode].append(thisnode)
+        path.pop() # remove thisnode from path
+        return closed
+
+    path = [] # stack of nodes in current path
+    blocked = defaultdict(bool) # vertex: blocked from search?
+    B = defaultdict(list) # graph portions that yield no elementary circuit
+    result = [] # list to accumulate the circuits found
+    # Johnson's algorithm requires some ordering of the nodes.
+    # They might not be sortable so we assign an arbitrary ordering.
+    ordering=dict(zip(sorted(G),range(len(G))))
+    for s in sorted(ordering.keys()):
+        # Build the subgraph induced by s and following nodes in the ordering
+        subgraph = G.subgraph(node for node in G 
+                              if ordering[node] >= ordering[s])
+        # Find the strongly connected component in the subgraph 
+        # that contains the least node according to the ordering
+        strongcomp = nx.strongly_connected_components(subgraph)
+        mincomp=min(strongcomp, 
+                    key=lambda nodes: min(ordering[n] for n in nodes))
+        component = G.subgraph(mincomp)
+        if component:
+            # smallest node in the component according to the ordering
+            startnode = min(component,key=ordering.__getitem__) 
+            for node in component:
+                blocked[node] = False
+                B[node][:] = []
+            dummy=circuit(startnode, startnode, component)
+
+    return result
+
+
+G = nx.DiGraph()
+
+for edge in sys.argv[2:]:
+    v1,v2 = edge.split(',', 1)
+    G.add_edge(v1,v2)
+
+for c in simple_cycles(G):
+    print " ".join(c[:-1])