Wednesday, 26 June 2013

AllRingsFinder, Sport Edition

Recently I've been looking at efficient perception of cycles in graphs. The goal is to improve and unify ring perception methods in the Chemistry Development Kit (CDK). Slowly working through, optimising and adding algorithms there was one I was dreading. The AllRingsFinder in CDK provides all simple cycles (every possible ring) in a graph by progressively collapsing edges (Hanser et al, 1996). The AllRingsFinder is used (sometimes incorrectly) throughout the library in many different modules.

To perceive all simple cycles can take a long time and I was pessimistic about making any improvement with it's speed. Fortunately, I was pleasantly surprised.

I'll go through the technical aspects of the implementation in the next post. This post will primarily focus on the performance. It should be noted that the algorithm is exactly the same and only the data structures and operations were modified.

Preliminaries

We will be testing the new implementation named AllCycles against the existing AllRingsFinder and HanserRingFinder. The HanserRingFinder was written by Rich Apodaca (depth-first) and is included in the CDK SMSD module. The AllCycles takes as input an adjacency list graph with the vertices described by integers. This is standard for graph algorithms - the time to convert the CDK molecule (IAtomContainer) is included in all results.

The benchmarks were run with Java 7 (21) on a 2.66 GHz Intel Core i7.

Complete Graphs Benchmark

When each vertex in the graph is adjacent to every other vertex we have the worst case for the algorithm. The number of simple cycles in these graphs grows very rapidly. The graphs are known as complete graphs and are named Kn, n being the number of vertices. The table below summarises the results of each implementation in milliseconds.

AllRingsFinder' is a modified version of AllRingsFinder that does not convert the paths to IRing objects. This conversion is currently inefficient and can impact the performance when there is a large number of cycles. The AllCycles and HanserRingFinder implementations already return the paths as a sequence of integers and atoms respectively.

Edit - I realise now the reason for the difference here is that I accidentally used the default and not silent domain objects. The default objects fire events on updates which can cause massive performance hits - particularly when adding atoms. Using the silent objects (no events) yields similar times to AllRingsFinder' without the IRing conversion. In short - don't use the default objects, ever... you don't need it, trust me. I'll leave the timings as reminded not to do this - the dataset results presented were already using the silent objects.

Graph Number of Cycles Implementation Time (ms)
K3 1 AllCycles 0.00 ms
HanserRingFinder 0.00 ms
AllRingsFinder 0.15 ms
AllRingsFinder' 0.00 ms
K4 7 AllCycles 0.00 ms
HanserRingFinder 0.00 ms
AllRingsFinder 0.95 ms
AllRingsFinder' 0.00 ms
K5 37 AllCycles 0.00 ms
HanserRingFinder 0.15 ms
AllRingsFinder 2.20 ms
AllRingsFinder' 0.75 ms
K6 197 AllCycles 0.10 ms
HanserRingFinder 0.85 ms
AllRingsFinder 7.80 ms
AllRingsFinder' 2.00 ms
K7 1172 AllCycles 0.55 ms
HanserRingFinder 26.10 ms
AllRingsFinder 193.80 ms
AllRingsFinder' 39.80 ms
K8 8018 AllCycles 12.0 ms
HanserRingFinder 817.00 ms
AllRingsFinder 10078.75 ms
AllRingsFinder' 1360.75 ms
K9 62814 AllCycles 994.70 ms
HanserRingFinder 54524.85 ms
AllRingsFinder 389597.00 ms
AllRingsFinder' 81699.53 ms

Performance of each implementation on complete graphs - log scale
It's quite amazing that using the same algorithm and simply changing the data structures (and not firing object events) can produce more than a 300000% improvement in performance (K9).

Dataset Benchmark

Examining the worst case gives us a concise view of how the implementation scales. To judge the real performance several datasets were used; ChEBI, HMDB and a subset of PubChem Compound (CID 0 - 25000).

There are molecules in the datasets with far more cycles than we can reasonably handle. Such molecules are considered impractical and depends on some predefined escape clause. 
  • AllRingsFinder measures the time when the algorithm starts and terminates if a set time limit is exceeded.
  • AllCycles terminates when a preset maximum degree (connections) on a vertex is exceeded.
  • HanserRingFinder does not have a timeout and so was not tested on the datasets.
The table below summarises the results for each dataset. To allow a reasonable number of repeats the timeout for AllRingsFinder was reduced from the default 5000 ms to 50 ms. Running on ChEBI with the default timeout took 389597 ms (29 impractical) and clogged the memory so much the timing was not stable. The maximum degree used for AllCycles was arbitrarily chosen as 500. Additionally the findAllRingsInIsolatedRingSystem() was used instead of findAllRings(). This allowed nearly identical prepossessing for both implementations using the fast RingSearch. Invoking the findAllRings() directly will yield slower results than those listed below. 


Data Molecules Implementation Time (ms) Impractical Cycles
ChEBI 29954 AllCycles 256 ± 56 ms 36 155476
AllRingsFinder 10324 ± 189 ms 86.86 131124.72
HMDB 37495 AllCycles 390 ± 48 ms 2 93148
AllRingsFinder 5410 ± 180 ms 40.83 86230.38
PubChem Compound Subset 22676 AllCycles 195 ± 55 ms 1 61229
AllRingsFinder 1942 ± 173 ms 5.05 60587.58

We can see that for ChEBI we were able to perceive 20000+ more cycles in ~2% of the time. As mentioned the times include the preprocessing and conversion, were molecules loaded direct as adjacency lists the time for PubChem Compound subset (195 ms) would be even faster! The counts for AllRingsFinder are fractional due to the fluctuations between repeats.

Choosing a Maximum Degree Limit

The performance of new implementation depends on choice of the maximum degree. The limit doesn't really have a meaning and may be chosen incorrectly. Set the limit too low and we lose information, too high and we may generate more information than we can handle . 

Running the implementation on all of PubChem Compound (Dec '12) with a limit of 5000 we can retrospectively calculate the limit for a given percentile. In PubChem Compound there were 17842848 ring systems with more than one ring,  987 (0.005%) of which exceeded our hard maximum degree of 5000.

Plotting the cumulated percent (T) against the maximum degree reached indicates that 99% of ring systems can be completed with a limit of just 9. Due to the sheer number of systems, if we did choose that value there would still be 166131 systems which were considered impractical and did not finish


To reduced the number of impractical molecules we need to go to fractions of a percent.

Maximum DegreePercent (%)Completed
(ring systems)
Impractical
(ring systems)
7299.95178340138835
8499.96178358766972
12699.97178376925156
21699.98178392933555
68499.99178410651783
88299.991178413421506
106299.992178414291419
144099.993178416021246
307299.994178417891059

This breakdown provides a meaningful insight for an appropriate value. I haven't decided yet where to include a Run Forever option - due to the combinatorial explosion there will always be molecules which we can't compute. I'll leave you with this entertaining video which explains the amazing power of combinatorial explosion.



Results