Recently I had to quickly come up with Python code that found the longest path through a weighted DAG (directed acylic graph). This particular DAG had a single start node and a single end node, but there were multiple weighted paths from the start to the end. "Longest path" in this context meant the path with the highest sum of weights.
To begin with, I coded up an exhuastive algorithm for trying all possible paths. This wouldn't be fast enough in the general case, but it didn't take long to write and is useful for comparison.
Internet rumour (contradicted by Wikipedia I now see) had it that the well-known Dijkstra's algorithm for the shortest path could be adapted for longest paths (but only in a DAG) by using the negative of the weights. However, a test of a couple of libraries for doing so indicated that negative weights were not supported (segfault with R for example), and I didn't want to code Dijkstra myself, so I decided to hunt around the internet for a direct solution to the longest path problem.
Geeksforgeeks (appropriately enough I guess) had an lucid description of a 3-step algorithm which I implemented. The only snag was that Step 1 involved a topological sort. Again, I didn't want to spend any time figuring this out so I looked for an implementation in Python and found one at ActiveState's Python recipes. I note that another popular Python implementation is available at PyPI.
Here's the code (if you are using an RSS reader to read this, you may have to click through to the original blog post):
Saturday, 30 November 2013
Thursday, 28 November 2013
Got CPUs to burn? Put 'em to work with GNU parallel
I've just been using GNU parallel for the first time. It makes running jobs over multiple CPUs trivial.
In the past, if I had a large number of single-CPU computationally intensive jobs and multiple CPUs to run them over, I would create separate bash scripts for each CPU with a line for each calculation, e.g. ./runthis input1.smi > output1.txt. This is not super-ideal as different jobs take different lengths of time and so any CPU that finishes its bash script ahead of schedule just sits there idle. It also involves making N separate bash scripts.
Enter GNU parallel. This comes with several Linux distributions but on Centos I just quickly installed from source. Once done, you just need to put all of the jobs in a single script and pipe it through parallel:
In the past, if I had a large number of single-CPU computationally intensive jobs and multiple CPUs to run them over, I would create separate bash scripts for each CPU with a line for each calculation, e.g. ./runthis input1.smi > output1.txt. This is not super-ideal as different jobs take different lengths of time and so any CPU that finishes its bash script ahead of schedule just sits there idle. It also involves making N separate bash scripts.
Enter GNU parallel. This comes with several Linux distributions but on Centos I just quickly installed from source. Once done, you just need to put all of the jobs in a single script and pipe it through parallel:
cat myjobs.sh | parallel -j7 # e.g. for 7 CPUsThere are a whole load of complicated ways of using parallel. I'm very happy with this one simple way.
Wednesday, 20 November 2013
Anatomy of a bug fix
There's been a longstanding problem with the 2d layout in Open Babel. When support for correct layout of cis/trans stereo was added, it was done as a correction to the initial layout rather than as an intrinsic part of the layout. This means that after nicely laying out a molecule to avoid overlaps, it may then flip two groups at the end of a double bond, crashing one part of a molecule into another.
At the time, given that the depiction code had just moved from no support for stereo to supporting it correctly (at a small expense of layout), sorting this out was way down my priority list. Anyhoo, a couple of days ago I decided to revisit this issue. Here's an example before and after:
I thought it might be interesting for some (or perhaps not) to see behind-the-scenes how one might go about sorting out something like this. In this case, I already knew the technical fix, but it still took some time to ensure it was done right. It's not quite #realtimechem, but probably as close I'll get.
Day one
1. Update local git repo by pulling from the main branch
2. Create (and switch to) branch for bugfix
3. Build Open Babel under MSVC (debug and release)
4. Some build issues with recent code (presumably only tested on Linux) for chemicaljson and chemdoodle plugins. Ignore for now as it doesn't impact this work (note to self: fix).
5. Need to generate problematic test cases. Convert several thousand ChEMBL molecules to SMIs and use grep to remove those with "." and keep those with either "\" or "/". Write Python script to sort by size and added line number as title (these are useful for debugging; smaller test cases are easier to handle).
6. Add print statements to the code to print out the measured "overlap error" before and after calling the code to fix the 2D layout
7. Use obabel to convert my test cases to SVG and pipe the output to a file
Night came. Morning came. The second day.
8. Coding time: move the double-bond correction to immediately after the initial layout.
9. When doing the somewhat exhaustive rotation around bonds for layout, do not alter bonds with defined double bond stereo.
10. Add print statements to the code to print out the measured "overlap error" afterwards
11. Use obabel to convert the test cases to SVG and pipe the output to a file
12. Write a Python script to compare the overlap errors for the old and new code
13. Visually inspect some of the cases where an improvement is supposed to have occured.
14. Oh oh! I see an example where the layout has improved but the stereo is now wrong.
Night came. Morning came. The third day.
15. I run the problem case in the debugger and step through the code inspecting the molecule object in the depiction code. Funny thing is, my code appears to be correct but none of the bonds are marked with stereo refs and so it doesn't do anything.
Night came. Morning came. The fourth day.
16. I run the problem case in the debugger with the original code. The stereo refs are there.
17. I update to the new code and run it again. The stereo refs are there originally but must be disappearing later. I step further along. A copy is made of the molecule; it seems like it uses its own special clone method. Looking further along, the bonds are cloned but the stereo refs are not being copied at that point. I fix this and all is well.
18. I run the test cases again, run my comparison Python script and start visually inspecting problem cases.
19. Seems like everything is working better than before.
20. I write a blog post before I forget everything.
21. Still need to remove the print statements, inspect the diff, commit the code, and send a pull request.
At the time, given that the depiction code had just moved from no support for stereo to supporting it correctly (at a small expense of layout), sorting this out was way down my priority list. Anyhoo, a couple of days ago I decided to revisit this issue. Here's an example before and after:
Day one
1. Update local git repo by pulling from the main branch
2. Create (and switch to) branch for bugfix
3. Build Open Babel under MSVC (debug and release)
4. Some build issues with recent code (presumably only tested on Linux) for chemicaljson and chemdoodle plugins. Ignore for now as it doesn't impact this work (note to self: fix).
5. Need to generate problematic test cases. Convert several thousand ChEMBL molecules to SMIs and use grep to remove those with "." and keep those with either "\" or "/". Write Python script to sort by size and added line number as title (these are useful for debugging; smaller test cases are easier to handle).
6. Add print statements to the code to print out the measured "overlap error" before and after calling the code to fix the 2D layout
7. Use obabel to convert my test cases to SVG and pipe the output to a file
Night came. Morning came. The second day.
8. Coding time: move the double-bond correction to immediately after the initial layout.
9. When doing the somewhat exhaustive rotation around bonds for layout, do not alter bonds with defined double bond stereo.
10. Add print statements to the code to print out the measured "overlap error" afterwards
11. Use obabel to convert the test cases to SVG and pipe the output to a file
12. Write a Python script to compare the overlap errors for the old and new code
13. Visually inspect some of the cases where an improvement is supposed to have occured.
14. Oh oh! I see an example where the layout has improved but the stereo is now wrong.
Night came. Morning came. The third day.
15. I run the problem case in the debugger and step through the code inspecting the molecule object in the depiction code. Funny thing is, my code appears to be correct but none of the bonds are marked with stereo refs and so it doesn't do anything.
Night came. Morning came. The fourth day.
16. I run the problem case in the debugger with the original code. The stereo refs are there.
17. I update to the new code and run it again. The stereo refs are there originally but must be disappearing later. I step further along. A copy is made of the molecule; it seems like it uses its own special clone method. Looking further along, the bonds are cloned but the stereo refs are not being copied at that point. I fix this and all is well.
18. I run the test cases again, run my comparison Python script and start visually inspecting problem cases.
19. Seems like everything is working better than before.
20. I write a blog post before I forget everything.
21. Still need to remove the print statements, inspect the diff, commit the code, and send a pull request.
Saturday, 16 November 2013
GaussSum 3.0 now out
I've just made a release of GaussSum 3.0. This marks the first time I've released code that uses Python 3 so I'm interested to see whether this will affect take-up.
The files are available at SourceForge, where there are also install instructions.
Here are the release notes:
This release has no new functionality apart from updates to the parser. The main reason for this new 3.0 release is that the codebase has been modernised so that GaussSum will continue to work in future. Specifically, GaussSum is now a Python 3 application, and it uses Python's Matplotlib for all graphs instead of Gnuplot. Also the dependency on Python Imaging Library has been removed.
Unfortunately, because of the move to Python 3, Linux users that are not on fairly recent distributions may not find either python3-numpy or python3-matplotlib in their package manager. If this is the case, you'll either have to update your distro or try compiling them yourself.
The Windows release is now a self-installing .msi courtesy of Python's cx_Freeze.
Since I'm now using Matplotlib instead of Gnuplot, if you want to spruce up your staid scientific plots à la xkcd, uncomment "# pyplot.xkcd()" in mpl.py. Then you can make graphs like the following:
May the best drug discovery tutorial win
Like a perfect alignment of stars, the Teach-Discover-Treat initiative touches on several areas of interest of mine. Its aim is to encourage the development of tutorials for drug discovery that involve freely-available software. The winners get a cash prize plus partial reimbursement of expenses to attend the ACS National Meeting where they will present the results. The competition is now open over at the TDT website.
I was fortunate to be present at the recent Spring ACS in New Orleans where the winners of the inaugural competition were presenting: Paulo Tosco, David Koes, George Nicola and Khaled Elokely. Unfortunately the tutorials themselves do not seem to be directly available on the web (and so will not be found by Google) but after agreeing not to be unethical you can download them from this page.
The results of the current competition will be presented at the ACS in San Francisco next year. Maybe I'll see you there.
I was fortunate to be present at the recent Spring ACS in New Orleans where the winners of the inaugural competition were presenting: Paulo Tosco, David Koes, George Nicola and Khaled Elokely. Unfortunately the tutorials themselves do not seem to be directly available on the web (and so will not be found by Google) but after agreeing not to be unethical you can download them from this page.
The results of the current competition will be presented at the ACS in San Francisco next year. Maybe I'll see you there.
Wednesday, 23 October 2013
Summary slides should be summarily skipped
This is a summary of my blog post. It will begin with an introduction. After this, I will present the meat of it. I will then surprise you by discussing this. Finally, guess what, I will end with a pithy conclusion.
Ok - so you think I'm overdoing it. But this is exactly how summary slides in presentations now appear to me. Once I thought they were essential; now I think they take away from the presentation.
Think about it. The audience is sitting there all ready to hear about the work you're doing; they are at their most alert and ready to engage with your topic and try to understand your results. And what do you do? You show them a slide that sends them to sleep. Either it's that pro forma slide that we've all seen before, or you can't resist talking through every single aspect of the entire talk so that there's nothing to look forward to; no-one wants to hear the same talk twice.
Perhaps there's an argument to be made for a summary slide for longer talks, especially if it's the last 30 years of some professor's work and it's going to be rambling about a bit so a bit of steering could be useful. However (thankfully) the majority of talks these days are of the 25-30 minute variety and that time spent on a zero-information-content summary slide may translate into cramming the last five slides of results into 30 seconds.
Anyhoo, this is on my mind at the moment as I'm putting together a talk for the GCC in Fulda in a few weeks time. Thing is, I'm one slide short at the moment. Hmmmm...
Ok - so you think I'm overdoing it. But this is exactly how summary slides in presentations now appear to me. Once I thought they were essential; now I think they take away from the presentation.
Think about it. The audience is sitting there all ready to hear about the work you're doing; they are at their most alert and ready to engage with your topic and try to understand your results. And what do you do? You show them a slide that sends them to sleep. Either it's that pro forma slide that we've all seen before, or you can't resist talking through every single aspect of the entire talk so that there's nothing to look forward to; no-one wants to hear the same talk twice.
Perhaps there's an argument to be made for a summary slide for longer talks, especially if it's the last 30 years of some professor's work and it's going to be rambling about a bit so a bit of steering could be useful. However (thankfully) the majority of talks these days are of the 25-30 minute variety and that time spent on a zero-information-content summary slide may translate into cramming the last five slides of results into 30 seconds.
Anyhoo, this is on my mind at the moment as I'm putting together a talk for the GCC in Fulda in a few weeks time. Thing is, I'm one slide short at the moment. Hmmmm...
Thursday, 10 October 2013
QM up for testing - The Quantum Chemistry Speed Test
Ever wondered which quantum chemistry package is fastest? No? Well you're not alone - I can't find anywhere on the webs a comparison of the speeds of different QM packages. Enter the Quantum Chemistry Speed Test...Over a series of weeks (weeks which may be spaced months or years apart depending on the ebb and flow of life), I will carry out the same calculation using a variety of packages. The calculation will be a geo-opt of a small size organic molecule on a single CPU and without any attempt to tune.
And you can play along at home. I'll be making all of the input and output files available for your viewing pleasure. Commercial software is unlikely to feature in this comparison (as I don't have access to any) but don't let that stop you. Note that for the usual reason, you should avoid publishing Gaussian timings.
The reason I'm interested in this is that it seems that the focus of many QM packages these days is towards carrying out massively-parallel super accurate calculations. But what I'd really like (and I think most users would be with me in this) is faster speeds for standard calculations. For example, in a project with Geoff Hutchison a few years back I carried out 1000s of single-CPU calculations using Gaussian on a 8-cpu per node cluster. If I had run these in parallel it would have been much slower (see Amdahl's Law) and those CPU hours would not have stretched so far.
Maybe if there were more of these speed tests, it would encourage developers to bring more performance to single-CPU calculations. Well, probably not, but it'll be fun to find out.
Image credit: Paul Townsend on Flickr
Subscribe to:
Posts (Atom)