I should add there is more than one (quick) way to detect cycles, the end implementation definitely is nothing new but it was certainly fun to program.
Cycle Detection in the CDK
In the Chemistry Development Kit (CDK) the SpanningTree class is the usual choice for cycle detection. It is used to obtain a molecule (aka AtomContainer) that only contains cyclic atoms. This molecule is then tested as to whether each atom is pressent.
Code 1  Using the existing SpanningTree
// create the spanning tree and obtain the cyclic fragments SpanningTree tree = new SpanningTree(molecule); IAtomContainer cyclic = tree.getCyclicFragmentsContainer(); // test whether the cyclic fragment contains each 'atom' for(IAtom atom : molecule.atoms()){ cyclic.contains(atom); }
The method is relatively fast and is used in some core concepts, including, aromaticity detection and atom typing.
Depth First Target Search
I'll come back to why the spanning tree is important and what it actually is later. If we dig down in the code we find the spanning tree uses a depth first target search to test whether a vertex is in a cycle. This mirrors the core idea of a cycle in an undirected graph.
An undirected graph has a cycle if and only if a depthfirst search (DFS) finds an edge that points to an alreadyvisited vertex (a back edge)  Wikipedia
The exact piece of code which is used is a depth first target search. Lets have a look at this targeted DFS  I've restructured the code a little so it fits in the snippet, you can see the original here.
Code 2  SpanningTree internals
// see org.openscience.cdk.PathTools.depthFirstTargetSearch for original code boolean targetedDFS(IAtomContainer mol, IAtom atom, IAtom target, IAtomContainer path) { atom.setFlag(VISITED, true); // step 1 for (IBond bond : mol.getConnectedBondsList(atom)) { // step 2 IAtom connected = bond.getConnectedAtom(atom); // step 3 if (!connected.getFlag(VISITED)) { // step 4 path.addAtom(connected); // step 5 path.addBond(bond); if (connected == target) { // step 6 return true; } else { if (targetedDFS(mol, connected, target, path)) { // step 7 return true; } else { path.removeAtom(connected); // step 8 path.removeBond(bond); } } } } return false; // step 9 }
I'm going to step through each block and so it is clear what it does, I will add some commentary on various aspects after.
 indicate we are visiting atom
 access all bonds connected to the atom
 for each bond from the atom get the connected atom
 check we have not visited the connected atom  note we actually may be in the process of visiting it further up the stack
 add the connected atom to our path
 if the connected atom is our target return true  the path will contain all atoms we visited to get to the target
 connected atom was not the target, recursively visit the neighbours of the connected atom and attempt to find the target. If whilst recursively visiting we found the target, return true, the path still contains all atoms which lead to the target.
 we did not find the target whilst visiting the connected atom, connected atom cannot therefore lead to the target and we remove it from the path.
 we have visited all connected atoms and did not find the target, return false
It is important to note the atom/target can not be the same atom. For the purposes of demonstration let us presume this is allowed  this will become more clear when we look back at spanning tree. The following animation visually represents a targeted DFS on the atoms a, f and b. The visit order of the neighbours is determined by the order in which the edges are stored  the animations depict the worst possible route but this of course can occur.
A DFS from vertex a, the search terminates when we visit vertex b a second time. Vertex b was not the target, no path is returned. 

A DFS from b, the search finds the target (b) and the path of {b, h, g, e, d, c} (which are all in a cycle) is returned. 
There a couple of things to note about the complexity of the above process. Here I denote the number of vertices (atoms) as V and the number of edges (bonds) as E. This follows the standard Big O notation for graphs.
Firstly, the VISITED flag (step 1) needs to be reset for each subsequent search. This is a consequence of mixing the algorithm state with our data model. To reinvoke this method over the same set of atoms (which we need to do) we first need to reset all the flags. Reseting the flags takes time proportional to V + E. Secondly, to access the connected bonds (step 2) we need to perform a linear search (in coordinate list representation) taking at most E. Thirdly, to remove the atom and bond (step 8) we need to again do a linear search over all atoms/bonds taking at most V and E operations respectively. Phew, thats a lot of V's and E's floating about. It should be noted we need to do a targeted search for every vertex (atom). This would take a very long time, surely there is a way we can reduce the number of atoms we need to test? This is exactly where the spanning tree is used.
Spanning Tree
Okay, so what is a spanning tree? I'm not a graph theorist so my definitions may be slightly loose, the links below will provide a more accurate definition. Well a tree is a connected graph where every vertices is connected by one simple path  that is, there are no cycles. A spanning tree of a graph is the transformation of the graph (which may contain cycles) into a tree which contains only enough edges to span all connected vertices. In other words it does not contain redundant edges.
Valid spanning trees of our example graph. The dashed line indicates where an edge was removed. 
The spanning tree provides several interesting properties. Any edges that are not pressent in the tree are redundant and therefore much be in a cycle. If we obtain multiple trees from our graph, our graph is disconnected. The CDK uses an implementation of Kruskal's algorithm to calculate the spanning tree which runs in O(E log E).
Once the spanning tree is determined we now know which edges are in a cycle. We then use the depth first targeted search (above) between two vertices of any edges that were removed.
This process is efficient as we only ever check targets which must be in cycles. However, If we have two or more cycles we still traverse the graph multiple times. Can we detect all cycles and only visit each edge once...?
Changing the Representation
Firstly lets use the more friendly adjacency list for traversal. As shown previously the adjacency list allows quicker access to adjacent neighbours (step 2). When we create the adjacency list we index each vertex to an int value between 0 and V. This index allows us to quickly access the vertex value (atom) but it also allows us to decouple the algorithm visited state by storing a boolean array of size V(step 1). In fact we can do even better by storing which vertices we have visited in a BitSet or if the molecule is small enough on a single 64bit long value. Similarly as we do not need to know the order we can also store the path we visited as a BitSet (step 8).
All these changes sound great but we will still be visiting the same vertices multiple times and we still need to calculate the spanning tree. There is actually a far simpler way  which doesn't require the spanning tree at all. In truth I truncated the wikipedia quote at the start of this post. The solution should be clear from the first diagram but if not here is the full quote
An undirected graph has a cycle if and only if a depthfirst search (DFS) finds an edge that points to an alreadyvisited vertex (a back edge). Equivalently, all the back edges, which DFS skips over, are part of cycles. In the case of undirected graphs, only O(n) time is required, since at most n − 1 edges can be tree edges (where n is the number of vertices).  Wikipedia
Here O(n) means O(V + E)  Now, that's an algorithm that scales!
The animations below demonstrated what the targeted DFS was doing and what we can actually do.
With our adjacency representation we can do this very easily. The core method of the algorithm can be written in only a few lines of code.
Code 3  Is in ring
void check(int i, long parentPath, long path) { explored = path = setBit(stack[i] = path, i); for (int j : graph.get(i)) { if (isBitSet(parentPath, j)) rings = stack[j] ^ path; else if (!isBitSet(explored, j)) check(j, stack[i], path); } }
Okay I'm cheating a bit, lets make this a bit more legible and add some comments.
Code 4  Is in ring, verbose
// adjacency list representation List<List<Integer>> graph; // index of set bits indicate which vertices we have visited long explored; // index of set bits indicate the vertices in rings long rings; // keeps track of our state when we visit each vertex  each vertex is only // visited once long[] stack; boolean isInRing(int i) { // check if vertex i has been explored if (!isBitSet(explored, i)) check(i, 0, 0); // search from i to all connected components // check with the bit 'i' is set and thus in a ring return isBitSet(rings, i); } void check(int i, long parentPath, long path) { stack[i] = path; // store our state before we visit 'i' path = setBit(path, i); // set 'i' is visited in out current path explored = explored  path; // set 'i' as explored globally // for each neighbour of i for (int j : graph.get(i)) { // have we visited it before in out parent path? //  we use the parent path as one 'j' will be our // parent and this would falsely indiciate a cycle. if (isBitSet(parentPath, j)) // take the different between our current path and our state // right before we last visited 'j'. this difference is the // set of vertices we traversed to reach 'j' again. mark all // of these as in a cycle/ring. rings = stack[j] ^ path; // have we explored j? stops us revisiting parent vertices else if (!isBitSet(explored, j)) // check vertex 'j' passing the path before we // visited 'i' as the parent path and our current path check(j, stack[i], path); } }
Benchmark
To benchmark this I wrote a quick executable jar that takes an input SDF and tests how long it takes to find all the ring atoms. It prints out the answer from each method to confirm they performed correctly. It also includes a stress test which I'll come back to in a bit. The source code and cycledetection.jar can be found on github project. Here is an example usage (help will print all options):
Code 5  Runnable benchmark example
java Xmx1G jar cycledetection.jar i ~/Downloads/Compound_000000001_000025000.sdf
The times change a little but normally the BasicRingTester is about 800% to 1000% faster than the spanning tree. For 20,000 molecules the BasicRingTester usually takes about 120190 ms and the SpanningTree takes between 9001000 ms to complete.
Average (n=50) time taken to detect vertices which belong to cycles in PubChem Compound ID 0000125000 
I mentioned in the last post there is some overhead in converting our graph to the adjacency list. We can test how much this overhead is by removing the actual ring computation. This can be done by setting the simulate option.
Code 6  Runnable benchmark example 2
java Xmx1G jar cycledetection.jar i ~/Downloads/Compound_000000001_000025000.sdf simulate
If this option is set we'll see that is takes 90110 ms (on my test input) just to convert the actual graph and the true time of the computation is very quick at around 2050 ms. We can show just how quick the algorithm is if we scale up the input and perform the stress test.
The stress test loads a graphene molecule with 2,598 atoms in cycles and a single atom which is not.
Molecule loaded for stress testing 
Average (n=5) time taken todetect vertices which belong to cycles in a graphene molecule with ~2,500 atoms 
To run the stress test, yourself
Code 7  Runnable benchmark example 3
java Xmx1G jar cycledetection.jar stresstest i dummy
Final Thoughts
Although the above code is just for a boolean test (i.e. is this atom in a ring) this is still very useful. The simple calculation in the CDKHuckelAromaticityDetector use the spanning tree to simply get which atoms are in a ring system. Clearly it is possible to adapt the code to detect the independent systems. The technique also roughly preserves the concepts of a path it may well be possible to determine other cycle properties (i.e. size) just by looking at the sets of vertices as they are discovered.
Edit: added graphs to better demonstrate the performance of each method
Edit: this is now available using org.openscience.cdk.ringsearch.RingSearch.
This is a very nice approach  although I'm not a huge fan of the bittwiddling. Certainly it's hard to argue with the numbers, especially for larger systems.
ReplyDeleteAm I right that you couldn't recover the actual rings with this technique? That it only gives you a list of the atoms in a ring?
(Oh, and how do you do those animations  they're great!)
I agree, you should avoid bit twiddling. A boolean[] would be slightly faster when storing which atoms have been visited but as this search maintains a stack of the state at each vertex this might lead to excess memory usage. I use a BitSet in the 'jumbo' version which is more friendly but perhaps both should be wrapped in another class which hides the bitwise arithmetic.
ReplyDeleteI haven't implemented it yet but I believe it is possible to recover the edges at minimal cost and after you've determined what is in a cycle. It is possible to separate out the simple cycles from those which share any vertices. You can then skip the Minimum Cycle Bases (basis?) calculation unless you have a difficult case. Curiously the SpanningTree already provides this but it isn't used to bypass the MCB. It may also be possible to skip the MCB entirely for common cases in chemical graphs such as two or three cycles share an edge or vertex (presuming they are planar).
Heh, the animations weren't quick to make unfortunately. I drew the base graph and then had a separate layer for each red vertex. I then exported each frame selecting the different layers and put it all together using the first free gif maker I could find  http://www.onyxbits.de/giftedmotion