Some Stuff

McKay's canonical augmentation method explained for simple graphs

2011-08-31T17:38:00.005+01:00

The previous post talked about generating one type of combinatorial object (chessboards) using a method similar to that outlined by Brendan McKay in a paper called "Isomorph-free exhaustive generation" (J Algorithms, 26 (1998) 306-324.). This one will focus instead on simple graphs, which requires both parts of the method.

The canonical construction (or canonical augmentation) method has two components. Firstly, only one 'expansion' of a graph is tried at each step from the set of equivalent expansions. Secondly, the expansions are checked to see if they are the inverse of a 'canonical deletion' for that graph.

For an example of the first rule, consider this set of expansions of a 4-vertex graph on the left:

Each of the 5-vertex graphs on the right are shown with the newly added vertex and edges in red; the arrows are labelled by the added edge set - so {1:4, 3:4} means edges added from 1 to 4 and 3 to 4. The sets of vertices to add to - {{0}, {1}, {1,3}} - are representatives of the orbit of these vertices. For example, the orbit of {1} in G(4) is {1, 2} as these two vertices are equivalent in G(4) on the left.

This is now quite similar to the situation with chessboards : trying only minimal orbit representatives for extending an object. In McKay's paper, the process of generating child objects is split into 'upper' and 'lower' objects. An upper object is a pair where X is (say) a graph, and W is a set of vertices to connect to a new vertex. A lower object is a pair where v is a vertex to delete. This is illustrated here:

Click for bigger, as usual. There is a function shown between a lower object for X' and an upper object for X. This is the 'deletion' function, and its inverse is the important one : f^-1, the function that adds a new vertex by connecting it to all the vertices in W.

This process will generate isomorphic graphs, so there has to be a way to reject children that are not canonical. This is where the second part comes in ... unfortunately it is harder to describe.

Roughly, we need to check that the newly added vertex is the one that should have been added if it was canonical. To verify this, the child graph is canonically labelled (eg : see this post, or possibly this one) and then the code checks if the added vertex (under the canonical labelling) is in the same orbit as the last one. Kind of.

The upshot is that this code now produces results very similar to nauty (geng) for graphs up to 8-12 vertices. For the larger numbers, I started to restrict the maximum degree, to shorten the runtime. It's definitely not as fast as nauty, but not too bad. I still have the lingering suspicion that I might start missing graphs for larger spaces, but it's not bad, not bad at all...

Generating Chessboards With K-Independant Vertex Sets

2011-08-04T11:02:00.007+01:00

After looking at Julio Peironcely's poster on generating chemical structures, where he describes using Canonical Path Augmentation (I think due to Brendan McKay) I went looking for more about it. One thing I found was this talk/slideshow by Mathieu Dutour Sikirić - incidentally coauthor of a nice book on Chemical Graphs.

Anyway; that's the context. Now : about chessboards? Well one of the examples given for augmentation (or 'orderly') schemes of independent vertex sets. I'm not sure what made me think of chessboards for this, but I think it's a fairly standard simple toy example. Let me show an example for 3x3 boards:

So, these are all the three by three boards where no two black squares share an edge. In other words, if the board is considered as a grid-shaped graph, then these are the k-independent vertex sets. So the question is : how to generate these?

The simple way, of course, is just to fill in every square and eliminate those boards that have pairs of black squares across an edge. This is the 'brute force' approach, and scales badly : there are 2⁹ boards but we only want 20 of these, or just 3%. For 4x4 boards, this fraction is even smaller - 131 of 2¹⁶ boards which is only 0.19%. This is 0.0026% for 5x5 boards.

So, the number of boards increases rapidly as the size increases and any way of decreasing this large search space is essential. The approach outlined in Mathieu's talk is quite simple : only try sets (boards) that are the minimal representative in their orbit. Ok, so maybe that doesn't sound so simple :) but look at this:

Here are 3x3 boards, with a numbering (any one will do), the equivalence classes of the cells, and a set of orbits. Lets say we've just generated {0, 4, 6} : how do we check if it is the minimal representative? Well, so long as we have the automorphism group of the board we just apply each permutation in the group to the set of numbers and check to see if any are smaller. In this case, {0, 2, 4} is smaller, so we don't consider {0, 4, 6}.

It's really as simple as that. Of course, we have to check each newly added cell to make sure it is not adjacent, but that can be done without consideration of previously generated boards. This means we don't have to store solutions, which means it can be done in parallel without communication. See how many there are just for 5x5 boards:

I suppose that one last question about all this is : so what? Does this have anything to do with chemistry? Well, actually, it does. Consider the bond order assignment problem : the second image shows all the possible assignments, which you may notice has duplicates (3, 5, 10, and 12 for example). Also consider the line graph approach to double bond systems where again the second image shows a pair of colorings of line graphs. In fact, these are k-independent sets...

Oh, and there is code here

Ring Plate Visualisation

2011-06-12T21:48:00.004+01:00

One diagram I've often wanted to make was filled-in rings in molecules :

Mainly for the purposes of highlighting rings without highlighting the atoms involved. This image was made using the CDK renderbasic module, and a small toy AWTRenderingVisitor that fills in paths. I'm not sure if the current one does this...

The gist for the main drawing method is here but is really just stuff seen before. The custom generator for rings is probably more interesting - if very simple - and is here. Note that it doesn't look too nice with inner-ring double bonds:

and would look nicer if the double bonds were symmetric.

EDIT : Coloring by ring equivalence classes didn't do what I expected...

Shouldn't all the outer rings be in the same class? Steran is how I expect, though:

External Symmetry Numbers and Graph Automorphism Groups

2011-05-29T15:01:00.005+01:00

So, there was a question on BioStar about calculating the 'external symmetry number' of a molecule - something I hadn't heard of, but turns out to be something like the subgroup of rotations and reflections of the automorphism group of a graph. Since I have some code to calculate the automorphism group, I naïvely thought it would be simple...

The questioner - Nick Vandewiele - kindly provided some test cases, which ended up as this code. Although many of these tests now pass, they only do so because I commented out the hydrogen adding! :)

On the one hand, there are some recent improvements that try to handle vertex and edge 'colors' - in other words, element symbols and bond orders. For example, consider the improbable molecule C1OCO1 :

These are the three permutations that leave the carbons and the oxygens in the same positions; when you include the identity, that makes 4. Cyclobutane (without hydrogens!) has a symmetry group of order 8. Similarly, cyclobutadiene now gives 4 instead of 8.

So what goes wrong when there are hydrogens? Well, it's a deeper problem than just hydrogens, but it starts there. Consider methane : it has an external symmetry number of 12, but my code gives 24 - why? Well the main answer is 'inversion', look:

The permutation (0)(1)(2, 3)(4) just swaps hydrogens 2 and 3. This effectively changes the chirality of the molecule ... sortof. It's not actually chiral, but its a reasonable description of the transformation. Apparently, this does happen (another thing I didn't know; there are lots more :) according to this document, but quite slowly compared to rotations - "slower than 1 cycle s^-1".

This kind of pseudo-chirality will happen at any tetrahedral center. Or at any atom with 4 neighbours, I think - like XeF4, which is square planar. As an example, take this spira-fused ring system:

with a transform that swaps 7 and 9 but not the pairs (0, 5)(1, 4)(2, 3). Effectively this changes the parity at carbon 6. Somehow I doubt that this kind of 'movement' actually occurs in solution, but I could well be wrong. In any case, it seems likely that the external symmetry number is 2, and not 4.

In summary, it is probably not possible to calculate the external symmetry number correctly without 3D coordinates, or symmetry axes, or point groups. I have a feeling that the positional info could be recorded as a 3D combinatorial map which would give explicit orientations for atoms with four neighbours.

Colored Tree Paths to Represent Bond Order Assignements

2011-04-10T16:45:00.004+01:00

A couple of different sources pointed me to a paper by Dehof et al (Bioinformatics, 2011; doi: 10.1093/bioinformatics/btq718) about bond order assignment by various methods. One of them is the A*-algorithm, which finds paths. Not in the molecule graph, but in a kind of tree of possible bond combinations. Precisely, the tree has a layer for each decision, and a leaf for each combination of decisions.

For example, here is the usual example of a square four-cycle (G):

with the tree shown below. If we pick a particular path through the tree - say to leaf 5 - we get the sequence of colors in the box in the upper center of the image. This corresponds to the assignment on the upper right. Clearly, it has too many double bonds to be carbon skeleton, but still.

So, it's a nice, simple way to represent the set of solutions. Well, it could be even simpler but representing it as a tree lets the paper authors assign weights to the edges based on the chemical rules outlined in work they reference (Wang et al, 2006). What happens if you don't have any weights? Here are all the assignments from the tree:

where the assignments are labelled by the leaf number (0-15) and have a letter label (A-F). These letters stand for 'equivalence classes' - really just isomorphism classes. Oh, and the grey vertices have unlikely atom types for carbons.

Final post on Combinatorial Maps of Cuneane

2011-03-24T18:39:00.005Z

EDIT : Updated with new image.
Well, unless I can make a more horrifying diagram than this:

..but I think that's unlikely. An explanation - if such a thing is possible - is that the red/purple arrows correspond to a kind of inside-out operation. So the triangular map (M3) at the bottom can be converted to the square (M4) by choosing the purple darts (2, 14, 16, 5) as the boundary, and putting everything else on the inside. The cycle (1, 4, 9, 20, 22) is red because it is the boundary of the pentagonal map (M5).

So, it seems like ϕ(M4) is α(ϕ(M5)) - so the cycle (0,8,23,21,5) in ϕ(M4) is (1,9,22,20,4) in ϕ(M5). That is, because (α[0] = 1, α[8] = 9, α[23] = 22 ...). However ϕ(M3) = ϕ(M5), which is confusing. Oh well

So I have updated the image so that all three 'maps' have the same cycles. Which means they have the same ϕ and therefore the same σ. I guess this means I misunderstood what a CM actually is : you can get 'different' embeddings with the same σ. Here is a summary diagram:

Which is just the same diagram, without all the numbers. Colored arrows in the faces show cycles of darts chosen as the bounding faces according to the same colored arrows joining maps. Note that the 4-cycle in CM3 chosen as the bounds for CM4 has to be flipped. So does the 5-cycle in CM3 chosen as the bounds of CM5. Perhaps there is still a transformation wrong here somewhere.

ANOTHER EDIT : After flipping CM5, and labelling the cycles (A-F/a-f : clockwise is uppercase, anticlockwise is lowercase) it is better.

The transitions are quite simple, really. Take a clockwise face, and make it anticlockwise (eg: B->b, red arrow from CM3 to CM5) and make the anticlockwise outer face clockwise (d -> D). These changes are on the red/purple arrows between maps. Under each map is a summary of the cycles, which really just shows that the outer face is anticlockwise and the others are clockwise.

Cuneane Maps

2011-03-23T14:36:00.007Z

So, to continue about combinatorial maps, here are some more intricate diagrams. Firstly, an embedding of cuneane with a 4-cycle as the outer face:

The permutations below are just for reference. Anyway, with this embedding, the cycles of ϕ are (0,8,23,21,5)(1,2,11,7)(3,4,17,15)(6,12,9)(10,14,18,22,13)(16,20,19) which does indeed have one for each face, including the boundary. If you use a different embedding, naturally you get a different map:

Which is quite different, and has ϕ of (0,6,10,3)(1,4,20,22,9)(2,14,16,5)(7,8,13)(11,12,23,19,15)(17,18,21) - again, 6 cycles for the 6 rings. And one to rule them all, and in the darkness bind them, of course. The cycles of darts are shown in this composite image:

Some of the triangular faces are missing the dart labels, as they were getting too crowded. Well, the whole thing is too crowded, but still. Highlighted in red in each are the darts corresponding to the outer face of the other embedding. Not sure what it means, though.

Finally, some references:

1) doi:10.1.1.149.3832 citeseer link - it's a dissertation, not a paper, but interesting.

2) Signatures of combinatorial maps (direct link to pdf) S Gosselin, G Damiand, and Christine Solnon. Damiand is the author of some of the images on wikipedia about combinatorial maps...

Squashing Molecules Flat and Combinatorial Maps

2011-03-20T19:02:00.009Z

So the traditional way to say it is probably 'planar embedding', but I've gone for the alternate terminology of 'squashing flat'. For many molecules, this is a fairly easy process : some decisions have to be made about rotatable bonds, but quite a few drugs are just things sticking off a benzene ring or two.

For fully 3-dimensional molecular graphs (for example) it is trickier. I don't yet know a simple algorithm for choosing different possible embeddings ... er ... squashings. Perhaps they are all complicated; they seem to involve tree data structures with names like "SPQR Tree" (see "Optimizing over all combinatorial embeddings of a planar graph").

Far easier, though, is the intermediate data structure between a graph and some 2D coordinates - known as a combinatorial map. These mathematical objects store the detail of the embedding by recording the order of 'half bonds' called flags (or darts?) around a vertex. Maybe a flag is the vertex and the half-bond. Anyway, here is a picture of K4 (or a squashed tetrahedrane skeleton):

Now the left hand side just shows the graph with labelled vertices, while the one on the right shows the flags (I'm going to call them that now, sorry). At the bottom is the permutation σ applied to the set of flags (f). To understand σ consider just the central vertex (1) : if we travel clockwise round from flag 1, we get the flags (1, 8, 6). In fact, this is how I constructed the permutation : clockwise round all the vertices. It can be helpful to consider σ in cycle notation as (0, 2, 4)(1, 6, 8)(5, 11, 9)(7, 10, 3) which has a cycle for each vertex.

This image shows two different embeddings of the same graph. Well, I suppose technically it is different embeddings of the same labelling. In any case, L and R are flipped : if there wasn't so much symmetry it might be better. The same procedure as before is used to get σL and σR (clockwise both times). Now check this out:

What is this crazy thing? Well, the other part of the combinatorial map is an 'involution' called α, which is really just another permutation that stores the pair of flags for each edge. In the images above, I have just ordered the edges using the vertices (as 0:1,0:2,0:3,1:2,1:3,2:3) then ordered the flags by vertex index : that is, α = [1, 0, 3, 2, 5, 4, 7, 6, 9, 8, 11, 10].

Anyway, using the formula φ = σ ⋅ α we get a permutation of what I actually am going to call 'darts'. This is all on the wikipedia page, but it should be clear that the labelled arrows on the image go in cycles. So there is (0, 6, 3) for example. Indeed φ = (0, 6, 3)(1, 4, 9)(2, 10, 5)(7, 8, 11) which is - of course - the four faces of the embedding, including the boundary face.

Phew! Fun stuff, and there is code here as well as some tests here. Well more like just System.out statements, but the beginning of tests...

Further work on PDB hetdict

2011-03-06T21:05:00.002Z

Previously... After using setFormalCharge instead of setCharge, some of the null atom types disappeared. Specifically, the quaternary (SP3?) nitrogens, that have to have a +1 formal charge.

What is left? Mostly FeS clusters, and other metals. Turns out that the height-1 signature is clearer for 'clustering' the atoms into types. This is because a signature of this height records only the immediate neighbours of the atom.

One fairly frequent example is "[Co]([N][N][N][N])", but fortunately there is a patch for this based on a bug report. The atom is found in cobalamine and other cobalt-haemes.

Still, there are quite a few 'odd' ligands - either due to unusual elements, like Iridium; or unusual coordinations, like this iron with 6 oxygen neighbours. Hmmm, PDBeChem's image (at the bottom of the page) doesn't show all six bonds. It's clearer in Jmol from PDBSum, where the C3 symmetry of the iron ligands makes it seem (to me) that it is 6-coordinate.

Wow, that's a lot of links.

Null Atomtypes in PDB Het Dictionary now in some sort of RDF thing

2011-02-08T20:12:00.003Z

This blog post won "Best Title" in the award for "Random Combinations of Words" category.

Anyway, the atomtypes that were not found by the CDKAtomMatcher are now down to just 1,400. There were various errors on my part (such as using element symbols like "CL" rather than "Cl') and CIF files are trickier to parse than I thought (an atom name can be "HOH'" or 'P"' - look closely!).

Just to make life even harder, and because Egon suggested it, I tried to get the output in some kind of RDF using Jena. In practice, this seems to mean a sort of vague tree of nodes, some of which are data and some are types.

Well, its in the same repository as before, with source code. Although looking at the N3 file it still looks like it would be impossible to work out which atom belongs to which hetgroup...

Atom-typing the Hetgroup Dictionary

2011-02-07T18:36:00.003Z

So I posted this to CDK-devel, but probably this is the better place...

I've been trying to make a map between the atom IDs used in the HET dictionary (which is in CIF format) and atom types of some sort. To see what this looks like, here is a tail of the file:

ZZZ.O6A:O.sp2
ZZZ.H7C1:H
ZZZ.H7C2:H
ZZZ.H8:H
ZZZ.H2N1:H
ZZZ.H2N2:H
ZZZ.H3:H
ZZZ.H5:H
ZZZ.H6:H
ZZZ.H6A:H

'Zzzz', you may be thinking, but although many atom ids are quite obvious (like H8 is a hydrogen), some are probably not. One annoying aspect of this process was that the CIF file format is not especially friendly, and particularly, the file has 'loop_'s that don't terminate in octothorpes ('#'), as I thought they would.

Probably the parser (an IteratingCIFReader) could be much better written - in fact, it will probably only parse this one CIF! So my initial estimates of 13,000 typing failures is now down to only 3,508. What are the atoms that fail? There are some that are bound to like TBR, which is decidedly a non-organic molecule.

Quite a few, however, are aromatic chlorines:

Null type for 00A CL4A [CL]([C]([C]=[C]))
Null type for 014 CL [CL]([C]([C]=[C]))
Null type for 01A CL4A [CL]([C]([C]=[C]))
Null type for 01W N [N]([C]([C][C][H])[H][H][H])
Null type for 024 BR19 [BR]([C]([C]=[C]))
Null type for 032 CL13 [CL]([C](=[C][C]))
Null type for 039 CL [CL]([C](=[C][C]))
Null type for 055 CL1 [CL]([C](=[C][C]))
Null type for 062 CL1 [CL]([C]([C]=[C]))
Null type for 064 CL13 [CL]([C]([C]=[C]))
Null type for 064 CL15 [CL]([C](=[C][C]))
Null type for 064 CL25 [CL]([C](=[C][C]))
Null type for 088 CL32 [CL]([C](=[C][C]))
Null type for 088 CL37 [CL]([C]([C]=[C]))

The last part is the height-2 signature, just to give a quick idea of the environment of the atom. Aha! Some quick "cut | sort | uniq -c | sort -n" gives me:

24 [BR]([C]([C][H][H]))
29 [CL]([C](=[C][N]))
29 [CL]([C](=[C][S]))
38 [N]([C]([C][C][H])[H][H][H])
42 [N]([C]([C][H][H])[H][H][H])
48 [CL]([C]([C][H][H]))
75 [N]([C]([C][H][H])[C]([H][H][H])[C]([H][H][H])[C]([H][H][H]))
337 [BR]([C]([C]=[C]))
1176 [CL]([C]([C]=[C]))

gives the 'top-10' worst offenders. That nitrogen one is N(CH3)4 - perhaps charge is a problem here?

Final example for today : bowtieane

2010-12-20T20:52:00.004Z

The simplest-yet-most-complex example I could find is this one from another previous post (I like to recycle :)

I should mention again that I don't know if this is an actual chemical. I call it 'bowtieane', but maybe it has a proper name. In any case, it makes a nice small test case. The ring equivalence classes shown (A, B, C) are defined by the signatures as usual.

What's really strange though is SSSRFinder's partition of the rings. It divides the two A-rings into separate classes. This seems like the opposite behaviour to the fullerene-26 case where there were fewer equivalence classes than for the signature method. I really, really should read the papers referenced in the code.

Fullerene-26 revisited : more vertices, less complexity

2010-12-20T18:48:00.003Z

There is an example graph that makes for a better comparison between my code and the SSSRFinder. It is the 'fullerene-26' molecule that I used in a previous post. Unfortunately the diagram is wrong in that post; the one below has the correct vertex equivalence classes (colors):

The difference is just a pair of yellows that should have been pink (and v.v.) in the topmost C ring. Anyway, the graph is complex enough that there are different vertex symmetry classes without multiple bonds. Dodecahedrane, on the other hand, is so symmetric that without double bonds all vertices are the same.

So how does the SSSRfinder do? Well it considers all the rings to be different apart from two of the B rings. It is not reporting the final B ring, which could be considered as the bounding face of the map. In short, there are two equivalence classes (5-ring and 6-ring) which makes some kind of sense but isn't terribly informative.

Dodecahedrane has 12 faces, right?

2010-12-20T16:43:00.004Z

So there was this guy called Euler, and he had a formula that goes something like F = E - V + 2. Well, actually it is χ = V - E + F, where χ is the Euler characteristic, and this is equal to 2 for polyhedra. Anyway, the point is that dodecahedrane has 12 faces (cycles).

For the SSSRFinder, however, it has only 11; which is annoying. Moreover the ring equivalence class method only distinguishes based on the underlying simple graph - in other words it ignores bond order. In some applications this might be exactly what is needed, but I'm glad that my method gives a more detailed result:

So, apart from being a ridiculously detailed image, the above shows the face (ring, cycle) equivalence classes for dodecahedrane with a particular double bond network. Clearly any face could be 'glued' to another along one of the edges, following the vertex classes. All possible combinations of faces are shown in the 'face quotient graph' at the bottom right.

SSSRFinder and Ring Equivalence Classes

2010-12-19T22:21:00.005Z

With hubris and arrogance, I implemented my own ring equivalence class finder. With humility and grace, the SSSRFinder method getRingEquivalenceClasses does better for at least one case:

This is the test case used in SimpleCycleBasisTest, and probably in the paper on the subject. I guess I should read that instead of other papers on cycle bases. I do get the same answer for cuneane as when I do it by hand. That is to say, 3 equivalence classes with [1, 2, 2] elements.

For the 'Bauer' graph above (as I am currently calling it - Ulrich Bauer is the author of the SSSRFinder) the signature method puts almost every vertex in its own equivalence class. Now my method makes equivalence classes from cycles with the vertices in the same vertex equivalence class. Naturally enough, this makes a lot of ring equivalence classes.

The graph on the left has vertices colored by degree, to show how dissimilar the rings are. Is the SSSRFinder really getting the right answer?

Herschel graph with different planar embeddings

2010-12-18T21:30:00.004Z

Another example from a paper that I have mentioned before. This time the Herschel graph which is another of these crazy graphs thought up to prove or disprove some conjecture. The spring-layout paper gives two different layouts:

These really are the same graph :) One way to see this is to look at the vertex colors that I have applied. Each face is given a label (A or B) based on the key shown below the two graphs. The 'A' face is {Black, White, Grey, White} for example. The left-hand embedding (or layout) has such an A-face for its border, while the right-hand one has a B-face for a border.

Presumably, it is possible to convert a vertex partition (into equivalence classes) into a partition of faces. It seems easy for examples - like this - that have a planar embedding. More difficult for graphs that don't have one. On the other hand, not all planar layouts look very informative:

This is twistane (again) but not looking as symmetric as it can. However, the faces show the regularity - they are all the same, even the boundary. The colors used are the same as in the previous post.

The many faces of fused cycles

2010-12-16T19:29:00.004Z

Although many ring systems in molecules are quite easy to layout as 2D diagrams, there are some that are inherently 3D. Bridged rings are usually in this class; consider my favourite example molecule, cuneane:

Each of (A, B, C) is a particular layout of the same molecule, but with a different boundary (hexagon, pentagon, er...kind of fused squares). It would be nice to have a layout method that picked the same choice each time - regardless of the permutation of atoms and bonds. Even better if it could allow enumeration of the alternative possibilities.

As another example, consider a series based on twistane (which is a molecule) to two other graphs that may well not be actual molecules:

Twistane itself is in the middle, surrounded by five- and seven- ring equivalents. The upper layouts emphasise one ring in the graph while the lower ones emphasise the dual rings in each.

Why modular decomposition is not very useful for chemical graphs

2010-12-15T14:30:00.004Z

It is difficult to publish negative results in a journal, but a blog post seems like a good place to record the experience. Especially situations like this, where it probably should have been obvious not to try in the first place...

So; what is modular decomposition? Briefly, a module is a little like a connected component in a graph - indeed, a connected component is made up from one or more modules, but modules can overlap. Decomposition of a graph into its modules is, therefore, like finding the connected components of the graph. An example is shown here:

Two modules in the graph are circled, there may be others. The definition is a set of vertices that have the same neighbours outside the set. So, there was no need for me to make them complete graphs, but it looked nicer. Anyway, already looking at this example it is clear that these are not very 'chemical' graphs. They look more like networks (for example, see : [1]).

Indeed, I tried out some code made by the authors of a paper on a linear-time algorithm for modular decomposition [2] that was in java (makes a change from c). It was also nicely commented, unlike the mathematician - whose name I won't mention - that said he doesn't comment or document his code because "no program has ever improved through comments" which is just a lazy excuse, frankly.

The results for the molecules in the CDK MoleculeFactory were that almost all of them are prime modules; which means that they are elementary, or unbreakdownable. Notable exceptions are cyclobutane and a propellane-like graph (see image, modules are circled).

In fact, I suspect that chemical graphs with non-prime modular decomposition trees are rare. Partly because most graphs are irregular, but mainly due to the low degree (valence) of atomic vertices. Anyway, modules are not a solution to structure diagram layout [3].

[1] : J. Gagneur, R. Krause, T. Bouwmeester, and G. Casari. Modular decomposition of protein-protein interaction networks. Genome Biology 5:R57 (2004). doi:10.1186/gb-2004-5-8-r57

[2] : Marc Tedder, Derek Corneil, Michel Habib, Christophe Paul Simple. Simpler linear-time modular decomposition via recursive factorizing permutations. DOI:10.1007/978-3-540-70575-8_52

[3] Drawing Graphs Using Modular Decomposition DOI: 10.1007/11618058_31

Muhahahaha! Things can always be more complex...

2010-09-29T18:43:00.003+01:00

A good measure of how right a model or an implementation is can be how quickly it extends to more complex situations:

This is essentially the same except that 'SSE' has been added (no big deal) but also the leaf list has been generalised to AbstractLeafCollection. :)

I'm sure there are better ways to do this, but it fits neatly with some existing ideas I had on searching through lists vs searching through sets.

Rose Forests

2010-09-29T17:42:00.004+01:00

Carl Masak blogged about tree data structures, which caught my interest because of a pet-project of mine (tailor; a structure description and measurement tool) where I found myself using trees a lot. An awful lot. Perhaps ... too much.

Anyway, a related tweet by AudreyT mentioned an article called "Origami Programming" by Jeremy Gibbons. Which is in haskell (perhaps not surprisingly), a language I don't speak very well. However, while reading - and not understanding it - I did get one thing which was the idea of having a tree datatype where the node (called a 'rose') references a forest (a list of roses). I think that's right.

In any case, it solves a object-modelling problem for me that I had. The difficulty was that protein structures are hierarchical, yes, but have a strange mixed hierarchy of types. Perhaps this is obvious to haskell programmers and compiler-code writers, but this makes it very difficult to use the 'simple' tree datatype, where a Node class has a List of Node children.

Specifically, I mean situations like: a Chain composed of Atoms or a Chain of SSEs of Residues of Atoms. The 'rose tree' way of doing things makes this possible, at the price of a more complex model. Now here is a picture of a sketch of it:

So, for example, you can have a Protein:ChainList(Chain:AtomList(Atom),Chain:ResidueList(Residue...)) or several other possibilities. Also a visitor to the hierarchy can do separate things to a Tree than to a LeafList. Neat! Oh, and the code for the implementation (just the bare bones, not usable) is here.

Molecule Layouts

2010-09-20T20:15:00.005+01:00

I've been doing experimental work on layouts for the CDK. Not for atoms, exactly, for which the StructureDiagramGenerator is doing a pretty good job - could be better, of course, but what couldn't?

No, layout of MoleculeSets, and Reactions. Well actually IMoleculeSets and IReactions. With an ILayout class - my apologies to anyone who doesn't like generics, but it can be quite useful. Anyway, here is an example of what it is looking like at the moment:

Hmmm. Well, it is a grid I suppose. The problems with the ring bonds are known to me, please do not mention them >:|

The code for this is quite short:

IMoleculeSet moleculeSet = makeMolSet();
ILayout<IMoleculeSet> gridLayout = new GridMoleculeSetLayout(3, 3);
makeImage(moleculeSet, gridLayout, "three_by_three", 500, 500);

where the methods 'makeMolSet' and 'makeImage' do what you might expect (I hope :). Similarly:

IMoleculeSet molSet = makeMolSet();
AxisOrientation o = AxisOrientation.PLUS_X_PLUS_Y;
ILayout<IMoleculeSet> layout = new LinearMoleculeSetLayout(new StandardMoleculeLayout(), 3, o);
makeImage(molSet, layout, "xypos", 500, 500);

and this produces an image like this:

which is also ... alright. Getting there.

Consistent Zoom with Models of Different Scales

2010-09-12T16:36:00.004+01:00

So there is a way to get the zoom to work:

(to zoom on the picture, click for bigger :)

The approach taken here is to create graphical objects (LineElement, RectangleElement, etc) that are scaled at the origin, but not zoomed or translated to the center of the draw area. These last two parts of the transform are then added to the graphics transform.

One downside of altering the transform in the graphics is that if we want to draw extra stuff on the panel (like the detail string "Zoom = x, Scale = y" in the picture above) the original transform has to be captured before drawing, then restored after.

For example, see this commit:

http://github.com/gilleain/toyrenderer/commit/c3bad966cd37c604f2ab4eb0e177603e88bee2f8

Scaling and Text

2010-09-11T15:21:00.005+01:00

An obvious question about the CDK rendering code is : "Why not scale text with AffineTransform?" So, of course this is possible, and works quite nicely - but there is a cost.

One of the goals of the rendering code is to start from models of any scale, and render them as consistently sized diagrams on screen. By 'scale' here I simply mean the average distance between points. So the CDK layout code might use a distance of (say) 1 between two carbon atoms, but a file with a structure made in some other chemical editor might have an average atom-atom distance of 100. These are unitless values, by the way.

Now, what we could have done was transform the coordinates in the model to a consistent scale, then rendered these transformed coordinates. What we chose to do, however, was to calculate a single transform for the model and draw with this. If you use this transform to scale the graphics object before drawing you get this for a model scale of 10:

and this for a model scale of 100 (in other words a 'bond length' of 100):

the images are at slightly different zooms, but the point is clear I think. The choice is between scaling fonts independently, and altering the model directly.

Generic Rendering

2010-09-05T19:31:00.004+01:00

Egon++ is continuing the process of merging the CDK-JCP rendering core into CDK master. Some proposed generification of the classes was made on the mailing list, and here is a sketch of some of the classes and interfaces:
I realise that this looks horribly complex, but the question is : "Is it just complex enough, or too complex?". One of the things missing from the diagram is layout - there may be a need for classes like LinearMoleculeSetLayout or GridMoleculeSetLayout. Oh, and yes (you guessed it!) an IMoleculeSetLayout and ChemObjectLayout classes :)

The goal here is not to make convoluted code, but to avoid repeating stuff. A reaction renderer should know how to layout and paint molecule sets, and then pass on the task to the molecule set renderer, and so on. Some key things to avoid will be a) not to relayout on each paint and b) generate the diagrams in the right places, at the correct scales. I think that this will be possible.

How (not) to remove items from a (CDK) list

2010-09-03T11:45:00.005+01:00

So there I was, trying to remove all mappings from a Reaction like this:

for (int i = 0; i < reaction.getMappingCount(); i++) { reaction.removeMapping(i); }

and found that only half the mappings were being removed ... can you see why? :)

In fact, this is not some obscure CDK bug, but a logic error on my part. Equivalent code is this:

List list = getListSomehow();
for (int i = 0; i < list.size(); i++) { list.remove(i) }

using for example a java.util.ArrayList. The problem is that the index (i) is being tested against a changing number (the size). Once half the items have been removed, i is at the half-way point, so on the next pass it stops.

One way to 'solve' this is to go backwards:

for (int i = list.size(); i > 0; i--) { reaction.removeMapping(i); }

but this is slightly less clear than just :

List list = getListSomehow();
int size = list.size();
for (int i = 0; i < size; i++) { list.remove(i) }

which is clearer. Of course, even better is the List method removeAll(). It would be nice if Reaction had a similar method...

Combinations and Filters

2010-08-15T20:15:00.003+01:00

So there is now the beginning of a possible re-write of the DBST that uses basically the same approach, but is a bit more flexible. The code is here, but it's still a bit rough.

The original idea seems to have been to encode arrangements of double bonds for different ring sizes as a kind of 'library'. For each ring, a particular arrangement is picked until all possible combinations are generated. As a concrete example, see this example for a napthalene skeleton:

Here, the arrangements (1, 2) are applied to each ring (A, B) and then these are combined. Of the four combinations (A1B1, A2B1, A1B2, A2B2) only three are valid. The A1B2 combination has two atoms highlighted in red that have two double bonds and one single bond.

So one way to filter the combinations is to try and type the atoms, and reject any structure that has untypeable atoms. Another possible filter rejects structures that don't have atoms that are SP2 hybridized. Both of these are from the original code, but implemented as instances of a ChemicalFilter interface.

This is quite similar - uncoincidentally- to the approach in SMSD where graph-theoretical tools are used to generate possible subgraph matches, and then a chemical filter is used to rank the results. Ranking and filtering are not quite the same, so perhaps there should be a ChemicalRanker interface? It would be a little like an Enumeration, except that it might not be a total ordering, but a partial order.

1,4-Benzoquinone and the DeduceBondSystemsTool

2010-08-14T16:23:00.003+01:00

Once upon a time, there was a DeduceBondSystemsTool, and...

Er, anyway. Further to a patch made on the tool (patch ID : 3040138), there is a failing test for 1,4-benzoquinone:

The tool generates A, and the test wants B. Now, the problem is not that the tool is not trying B as a possibility, but that it generates A first and the final step doesn't remove it or rank it as better than A.

Understanding this requires an understanding of the algorithm. This is (roughly):

For each ring, generate a list of possible positions for all numbers of double bonds.
Generate a set of molecules by combining these positions together.
Remove 'bad' solutions and pick a solution with the least number of 'bad' N/S atoms.

where the definition of 'bad' is based on chemical rules like atom types.

Now, neither A nor B are bad solutions, and they don't contain N or S atoms, so they both have a rank of zero, and the first one generated will be returned. So, there is really no particular reason that the test should pass.

In general, it might be good to separate the generation of possible solutions from the ranking/filtering process. So that the computational or mathematical problem of generation is done by one class, while other classes determine which is the optimal solution (or set of solutions).

Line Graphs and Double Bonding Systems

2010-08-09T22:31:00.006+01:00

After looking at a CDK tool for fixing bond orders for aromatic systems (DeduceBondSystemTool in the smiles package) I wondered if there was a more general approach. That is, the problem is to take a molecular graph with no double bonds and generate all possible double bonded systems.

One possibility might be to first convert the graph (G) into a form known as a line graph (lg(G)) where every vertex in lg(G) is an edge in G. If these vertices are labelled to represent the bond order, then an aromatic system has a particular line graph. For example, here is benzene:

The dashed lines show the construction of the line graph, and the labels '-' and '=' mean single and double. Now obviously, the two resulting graphs are essentially the same, so it would be nice to remove this redundancy. An example of two different bonding systems comes from phenanthrene:

Which is great, but how to generate all non-redundant colorings of the line graphs? Since a line graph is just a graph, it can have a signature, and a signature quotient graph. This can then be colored:

However, in this example it is necessary to 'half-color' some of the vertices of the quotient graph ... which doesn't quite seem to work. The numbers in between the colored quotient graphs show how many line graph vertices are in each of the symmetry classes.

In any case, this is a bit of a toy problem, with only a partial solution, but here is a code repository for a sketch of the code. Note that the algorithm is missing!

CDK Export

2010-07-28T15:35:00.002+01:00

A short post just to get a picture in place. Someone else's code, my diagram..

Click for bigger, as usual.

Formula Debugger

2010-06-29T13:22:00.003+01:00

(Please click for bigger - there is quite a bit of detail :)

The left panel shows the search tree of structures; each circle is a structure, the red ones at the tips of some branches are fully connected and saturated. The blue branch has been selected, and is shown in the middle, with the structure at the tip at the top.

This solution is also highlighted (in red - colors are not very consistent) and shown in detail on the right. The right/middle panel is a traditional molecule layout, but with numbers in place of atom symbols. The lower right panel is the spanning tree of the atom highlighted in red in the upper right panel.

Clearly, there are still a lot of unnecessary structures generated since most branches are dead ends. However, this version does at least avoid duplicates - sadly it also misses a few structures :( Clearly refining partitions based on element symbols doesn't totally work...

CDK Signature implementation now in review

2010-06-20T18:28:00.002+01:00

What is it?

I have written about them quite a bit, but here is a quick summary : signatures are a little like SMILES, but also somewhat like HOSE codes. They are a description of the connectivity of a molecule, or an atom in the molecule. A more detailed description can be found in these papers by Faulon et al: [1], [3] or in this blog post by me.

The java implementation of this algorithm is a collaboration between Lars Carlsson (who wrote a C++ version) and me (who ported this version to java). However, I was also influenced by my previous attempt at a port from the c implementation by Faulon's group. There is an online service for using their program called "sscan" here. It also deals with stereochemistry.

What is it used for?

So, what can be done with all this new code? Here are some possibilities:

Smiles-like canonical strings that represent molecules. Note that signatures are considerably longer than smiles, but are guaranteed to work for cuneane, and indeed a broad range of graphs.
As with HOSE-codes (which can describe molecule connectivity up to different 'spheres') signatures can vary in height. Practically, this means an atom's environment can be described with different levels of detail.
Due to the canonisation of the structure, the core algorithm can be used to give a canonical labelling of the structure, which can be useful for atom-atom mapping of isomorphic structures.
Calculating signatures for all atoms of a molecule produces a partition of the atoms into sets of equivalent positions. This is useful for a variety of analyses of a molecule's graph structure.

Obviously, there are advantages and disadvantages of using signatures for these applications, compared to existing techniques. I am not sure, for example, of speed differences in using signatures to get a canonical representation.

How do you use it?

The MoleculeSignature class is a wrapper around an instance of an IMolecule and provides several useful methods, many of them from the base class AbstractGraphSignature. For example:

IMolecule thiazole = MoleculeFactory.makeThiazole();
MoleculeSignature moleculeSignature = new MoleculeSignature(thiazole);
System.out.println(moleculeSignature.toCanonicalString());
// Result = "[C](=[C]([N](=[C,0]))[S]([C,0]))"

This is the canonical signature for the whole molecule. To get this, canonical signatures are made for each atom, and the canonical one from the list is returned. To get all the signatures - rather, the equivalance classes (or 'orbits') - use the calculateOrbits method like this:


MoleculeSignature moleculeSignature = new MoleculeSignature(MoleculeFactory.makeQuinone());
for (Orbit orbit : moleculeSignature.calculateOrbits()) {
System.out.println(orbit);
}

which gives this output (the 'makeQuinone' method makes 1,4-benzoquinone:


[O](=[C]([C](=[C]([C,0](=[O])))[C](=[C]([C,0])))) [0, 7]
[C]([C](=[C]([C,0](=[O])))[C](=[C]([C,0]))=[O]) [1, 4]
[C](=[C]([C]([C,0]=[O]))[C]([C](=[C,0])=[O])) [2, 3, 5, 6]

which tells us that the two oxygen atoms ([0, 7]) are in the same orbit, as are the carbons attached to them, and that the other four are in another orbit. I have written about more complex examples of orbits : in C60 or in other fullerenes or in some other regular graphs. In practice, most chemicals will have automorphism partitions that are (nearly) discrete.

So, finally, an example of how to get the canonical labelling of a graph:

MoleculeSignature moleculeSignature =
        new MoleculeSignature(MoleculeFactory.makeCyclobutadiene());
System.out.println(Arrays.toString(moleculeSignature.getCanonicalLabels()));

which gives "[0, 3, 2, 1]" - essentially this is the permutation which gives a canonical arrangement of atoms.

Non-CDK implementations?

There are other chemistry projects other than the CDK, and it should be fairly easy to make a mychemlib.MoleculeSignature by subclassing signature.AbstractGraphSignature (and similarly for AtomSignature/AbstractVertexSignature). All the concrete classes need do is tell its superclass about the underlying molecule graph - getVertexCount, getConnected - and the MoleculeSignature has to act as a factory for the concrete AtomSignature instances via getSignatureForVertex.

The signature project is on github and has some of the maven machinery for building/testing/packaging. There are a couple of 'toy' implementations for chemicals and simple (mathematical) graphs.

Any feedback, suggestions, and so on are welcome. I am also happy to help with other people's implementations in the form of code or just hints. Enjoy!

Non-chemical example : Grinberg's graph

2010-05-16T13:43:00.003+01:00

While looking for an algorithm to lay out fullerenes as 2D graphs, I came across this paper (PDF). It describes an annealing/spring layout method for Schlegel diagrams. I don't think I can spare the time to implement it at the moment, but one of the graphs in the paper is this:

known as the Grinberg graph.

C60 double bonding networks

2010-05-14T16:56:00.005+01:00

C60 (or buckminsterfullerene) has no hydrogens, so it must have quite a few double bonds. I am beginning to understand that bond orders are in some sense a simplification, and I suppose that the bonding is in some way delocalised across the sphere. However, here is a picture of two different bonding patterns:

(click for bigger, as usual). The 'ChEBI' bonding pattern on the left is from the molfile in a ChEBI entry while the 'radial' bonding one on the right is bonded according to schemes from this site which has an interesting graph-theory perspective on fullerenes.

The radial bonding version has a simpler, layered structure like this:

Ok, so that's a slightly comical picture of the slices. What is also nice is the quotient graph for the ChEBI structure:

I didn't color this, but what is great is that it looks like a subgraph of a fullerene! Apart from the loops along the top and bottom, of course.

Fullerene symmetries

2010-05-13T22:08:00.005+01:00

Continuing the theme of colored graphs, some of the more interesting examples are fused ring structures, especially those with some symmetries, but not completely symmetrical. Fullerenes fit this description, for example this 26 atom example:

the distribution of colors might look a little odd, but the dark blue atom surrounded by three cyan atoms is actually repeated at the top - which is really the other side of the sphere. These kind of 3D molecules don't lay out very well with the CDK's layout code, so I used JChempaint instead to make a more symmetric 28-fullerene:

This is actually a screengrab of a crude viewer I put together (commit) that takes a molfile and calculates the signatures. Selecting a signature from the right hand list highlights it on the graph. Anyway, it's easier than making images like this by hand:

Except that the little graph at the bottom (which is a quotient graph) can't be drawn by the renderer as it has loop-edges.

Chemicals as colored graphs

2010-05-13T17:04:00.004+01:00

The interface between maths and chemistry can be tricky when it comes to terminology - sets (maths) have elements, chemistry has a different kind of element; graphs have colors which are usually just numbers, diagrams of chemicals have colors which usually relate to the element type of the atom, and so on.

So, for maximum confusion, here are two pictures of graphs (that could represent chemical connectivity) colored by equivalence class (determined by signature). The signature trees are also drawn with graphical colors, but these represent the integer colors in the signature, which are not the same as the colors used to indicate equivalence class. Firstly, a structure that the smiles algorithm is meant to have trouble with (but may not exist):

It looks quite strained, so I expect that it may not be possible to synthesise. Another multi-ring system is this one:

I don't even know what this one would be called, even if it did exist. Annoyingly, this structure triggers a bug if the two dark blue atoms are connected. This makes the graph 3-regular, but the yellow equivalence class is split, which shouldn't happen.

Orienting a Pyrene diagram

2010-05-12T14:45:00.004+01:00

Another example from Symmmetry in Chemistry, of drawing pyrene:

the standard (IUPAC?) orientation on the left doesn't show the symmetries of the molecule as well as the rotated version on the right. The letters indicate equivalent atom positions - oh, and aromatic indicators are missing :)

It seems like there should be a way to discover symmetry axis from the graph - without coordinates. In a similar way that some ring perception algorithms work by reducing the graph:

However, it is not obvious that this would be better than laying out the structure, then using the 2D coordinates to determine symmetries. Also not clear to me is how to choose which vertices to merge, and which to duplicate. On the left hand side of the image above, the fragments have duplicate (a), (d), and (e) vertices, but (e) has also been merged before being duplicated.

1,2-dichlorocyclopropane and a spiran

2010-05-12T11:25:00.005+01:00

As I am reading a book called "Symmetry in Chemistry" (H. H. Jaffé and M. Orchin) I thought I would try out a couple of examples that they use. One is 1,2-dichlorocylopropane :

which is, apparently, dissymmetric because it has a symmetry element (a C2 axis) but is optically active. Incidentally, wedges can look horrible in small structures - this is why:

The box around the hydrogen is shaded in grey, to show the effect of overlap. A possible fix might be to shorten the wedge, but sadly this would require working out the bounds of the text when calculating the wedge, which has to be done at render time. Oh well.

Another interesting example is this 'spiran', which I can't find on ChEBI or ChemSpider:

Image again courtesy of JChempaint. I guess the problem marker (the red line) on the N suggests that it is not a real compound? In any case, some simple code to determine potential chiral centres (using signatures) finds 2 in the cyclopropane structure, and 4 in the spiran. Since the code is not using a 3D structure, only a connected graph, it can't work out the spiran's S4 axis.

Stuck : Detailed Description

2010-05-11T21:29:00.004+01:00

Ok, so this is the detailed version of the previous post.

To recap; the structure generation code is still missing a vital piece - the canonical checking. I have been implementing Jean-Loup Faulon's algorithm for generation, but there is no precise algorithm given for canonical checking. Here the relevant paragraph from the enumeration paper:

"Checking for canonicity is a common procedure of orderly enumeration algorithms, the procedure guarantees that the graphs generated are nonisomorphic. ...To verify that a graph is canonical, one labels the vertices of the graph in all possible ways. The graph is canonical if the initial labeling leads to a list of edges that is lexicographically smaller than the lists obtained with all other labelings. In the present paper, we have implemented two algorithms to verify canonicity, Tarjan tree canonization algorithm if the tested graph is acyclic and McKay’s Nauty technique otherwise"

Okay, so I don't understand this for a couple of reasons : firstly, 'labelling all possible ways' sounds like an n-factorial (n!) operation; secondly, nauty does not lexicographically compare edge strings (as far as I know). Sadly, I know that I have misunderstood something, since my code doesn't work, and theirs does. :(

The technique used by nauty is also used by bliss and saucy and is known as iterative refinement of partitions. A 'partition' is a division of a set into subsets called 'cells', and refinement of a partition roughly means making a finer partition that has cells at least as small, if not smaller that the original. The end of the refinement process are 'discrete' partitions that are the same as permutations, since each cell has only one member. This is sketched in the image below:

A simple implementation of a partition refiner is in my repository. It is a slightly modified version of the algorithm from a book I've mentioned before (CAGES). It tries to deal with vertex and edge colours, although I don't think that it does so all that well. However, for simple graphs, it does indeed produce the automorphism group, as well as check for a canonical graph.

For example, the simple 5-cycle graph in the image might be canonical (or 'canonically labelled') if the refinement process produces the identity permutation as the first discrete partition. This will depend on the particular cell selection algorithm, and the choice made of how to arrange newly split cells.

The modifications I made to the CAGES refiner were done to align the canonical checking with the graph enumeration process. This is important, as edges are added to a graph in a characteristic order, and this has to produce graphs that will be accepted as canonical. This is illustrated in the following image:

The cycle graphs are laid out in a standard way on the left, and in a linear view on the right. The linear view emphasises the order in which the edges might be added. In the canonical scheme I have chosen, vertices are connected to the next possible partner. To put it another way, the resulting graphs will have a minimal edge 'length' when laid out as a linear graph.

What the partition refinement process is doing is the same as labelling all possible ways. The leaves of the refinement tree are permutations. Assuming the refinement process gives automorphic permutations, this is almost the same as labelling within the orbits of the atoms. It seems like partitioning the atoms by signature, then refining this partition should work, but it doesn't.

Any better suggestions or clarifications are very welcome :)

Stuck : The Summary

2010-05-11T19:03:00.002+01:00

A couple of people have asked how the structure generation stuff is going, and the short answer is that I am stuck. This post will give a short summary of the problem, and the next will give a much more detailed description.

So the problem is this : given an elemental formula (like C6H12) or a list of fragments plus a formula (like {2 * CH, 2 * CH2, 2 * CH3 : C6H12}) return all possible connected structures.

There are simple ways to do this, such as connecting every atom to every other atom, and removing duplicates. The downside is that this takes forever, because this procedure will make many, many isomorphic copies of each solution. At the final filtering step, an all-v-all comparison would have to be done on these many copies.

A better solution is to check each structure each time a bond is made, to see if it is canonical. Although I know how to do this in theory, it turns out to be more difficult in practice. For simple graphs, I have a solution that seems to work. Chemicals are not simple graphs, however, as they have elements and bond orders.

This is not code optimisation - without checking for canonical graphs, the running times for even quite simple problems are far too long. For more realistic, reasonable problems the code would be too slow to be useful.

Debugger

2010-05-09T15:21:00.003+01:00

Compatibility Table

2009-07-30T10:25:00.004+01:00

With the canonical code in place (thanks to open source :) the structure generation goal is much nearer. The first thing to improve is the bond compatibility.

This image shows cuneane (again!) and the bond compatibility table. The table is tricky to calculate, but relatively easy to understand; there is only one bond between atoms of type A - so there is a one in the cell (A, A).

Conversely, no atom of type C is connected to an atom of type A (see this for more detail), so there is a O in both (C, A) and (A, C). Note that the table is not symmetric, as can be seen with (A, B) = 1 and (B, A) = 2. This makes sense, in that an atom of type A is connected to two of type B, and yet an atom of type B is only connected to one of type A.

Crude port of Faulon's C-implementation

2009-07-22T13:55:00.002+01:00

I finally did the sensible thing and just ported over the c code. I don't particularly like doing this, as the resulting code is probably quite fragile, and unreadable - java spoken with a c 'accent' usually is.

However, it does reproduce the same signatures as the c code, so that's good. It only took 1 day to port, but 2 days to debug...

Check it out here : http://github.com/gilleain/generation/tree/master

Nicer tree pictures

2009-07-09T22:21:00.004+01:00

Here are some pictures of proper trees, or rather proper signatures, to celebrate getting the translator program to compile (it involved replacing some 'local' keywords with 'define's - who knew?).

Also, thought I would post the code for tree layout I'm using:


    public int layout(Node node) {
        node.y  = node.depth * ySep;
        if (node.isLeaf()) {
            leafCount += 1;
            node.x = leafCount * xSep;
            return node.x;
        } else {
            int min = 0;
            int max = 0;
            for (Node child : node.children) {
                int childCenter = layout(child);
                if (min == 0) {
                    min = childCenter;
                }
                max = childCenter;
            }
            if (min == max) {
                node.x = min;
            } else {
                node.x = min + (max - min) / 2;
            }
            return node.x;
        }
    }

basically, it lays out the leaves first, then returns their centers. Each non-leaf node uses the min/max values of its children to position itself. If the min and max are the same, it only has one child, so it is laid above. Otherwise, it is put at the center of the range.

Simple, probably nothing new, but does the job.

edit: Heh. Although, with large trees...

Tree Canonization Simplified

2009-07-09T12:27:00.009+01:00

While debugging the methods to make canonical signatures, I learned something about tree isomorphism from various sources, including Prof. Valiente's excellent looking book on trees.

One way of checking isomorphism is canonisation, since two trees are only isomorphic if they have the same canonical form. For simple labelled trees, it looks like there is an almost trivial way to get a canonical string representation. Say we have two trees:

The are rooted, labelled trees. So the conversion to a canonised string proceeds as follows; for each node, lexicographically sort the string form of the labels of its children, and return the concatenated string. In python this looks like:

which...is unreadable. hmmm. Wish there was a better way to get marked-up code into blog posts. Perhaps there is one, and I don't know of it. Anyway, the point is that it is very short.

edit : the code is...


def printSorted(node):
 if len(node.children) > 0:
   childStrings = [printSorted(child) for child in node.children]
   return node.label + "("+ "".join(sorted(childStrings)) + ")"
 else:
   return node.label

Square Grids, Cylinders, Spheres, and Toruses

2009-07-06T21:01:00.006+01:00

This is straying from the point; but if any graph can be described by canonized trees made from its subgraphs, then what are the properties of very large (regular) graphs? A grid, for example?

Starting with a square, and fusing squares together results in this situation:

two and three fused squares are similar, at the height-two tree level. With four rings, a new type c appears (b becomes b' with three rings). Beyond this point, any number of rings fused in a row like this has the same structure.

Moving to a second dimension of growth, a square grid (G_nn) has the following structure:

On the left are 'snapshots' of the advancing wave of the tree. Looks a lot like a breadth-first search, I suppose. On the right are the trees for each snapshot, with increasing heights. Although this is not shown, 8 of the 20 leaves of the third tree are duplicates. This should be clear from the third snapshot, which has only 12 (filled) circles on the 'wavefront'.

These trees can even be extended to grids wrapped around three-dimensional objects. Or, to put it another way, grids wrapped up as surfaces:

The diamond shapes with lines radiating out are the expanding trees. To the left of each surface is a very rough sketch of what the spanning tree would be like. The dashed lines indicate cutpoints where the wavefronts meet - these are duplicated on the trees. This is similar to the idea of 'gluing' surfaces in topology.

Warning : Abstraction!

2009-07-02T00:11:00.003+01:00

This is a throwaway mathematical point, that I am not qualified to make, but it looks like three of the previous examples (diamantane, twistane, and cuneane) have a very abstract connection when colored by signature:

what I mean by this diagram is that diamantane has atoms colored by (a) connected to both other (a) atoms, and to (b) atoms. Its (c) atoms are only connected to (b)s; the arrows could well be double-headed, by the way.

The most complex situation is cuneane, where each 'type' of atom is connected to another in its type and to two in another type. Adamantane would just look like : (a)-(b).

Interesting, but it doesn't get the signature canonization methods debugged any faster...

Adamantane, Diamantane, Twistane

2009-07-01T19:36:00.003+01:00

After cubane, the thought occurred to look at other regular hydrocarbons. If only there was some sort of classification of chemicals that I could use look up similar structures. Oh wate, there is.

Anyway, adamantane is not as regular as cubane, but it is highly symmetrical, looking like three cyclohexanes fused together. The vertices fall into two different types when colored by signature:

The carbons with three carbon neighbours (degree-3, in the simple graph) have signature (a) and the degree-2 carbons have signature (b). Atoms of one type are only connected to atoms of another - the graph is bipartite.

Adamantane connects together to form diamondoids (or, rather, this class have adamantane as a repeating subunit). One such is diamantane, which is no longer bipartite when colored by signature:

It has three classes of vertex in the simple graph (a and b), as the set with degree-3 has been split in two. The tree for signature (c) is not shown. The graph is still bipartite according to the degree of each vertex.

A different interesting case is a structure with the excellent name of twistane.

This structure is also 3-colored by signatures (if that's the right terminology...) but is not bipartite with respect to degree. Or, to put it more simply, there are carbons with two carbon neigbours connected together. It seems quite clear from looking at the colored structure on the right that atoms with the same color are similar in the context of the structure.

Cuneane

2009-06-28T15:41:00.005+01:00

While taking a look at the InChI discussion mail archives, I came across a discussion on graphs that are difficult to canonize (so called 'isospectral' graphs - I've heard the term before, but I haven't worked with eigenvalues of adjacency matrices, so did not pay them much attention).

Anyway, one such was this structure, cuneane:

which shows 3D and 2D representation of the molecule (reproduced from the wiki page), and an arbitrary numbering in the center. The three signatures A, B, and C are shown labelled by this numbering.

What's interesting about this particular example is the number of times that the same atom is represented in the signatures. For the signature called 'B', which is rooted at atom 2, the atom numbered 5 appears four times in the lowest layer of the tree. This naturally follows from the large number of rings of different sizes that make up cuneane - two 5-membered rings, two 4-membered rings, and two three-membered rings.

Two pass rendering

2009-06-26T09:42:00.007+01:00

So, there was a question on the cdk-devel mailing list about bounding boxes, reactions, and text. An unfortunate consequence of the new design is that the renderer will not calculate bounding boxes that can fully contain the text. Concretely, this is what it would look like (not made in JCP!)

The blue box is the bounds that would be created, which is minimal with respect to the atom centers. The black box is the bounds that should be created, if we respected the text size. The problem is, the size of text is not known until the point it is drawn. Or, more precisely, until we have some sort of GraphicsContext to ask about the width in a particular font.

So, a two-pass system was suggested. When this was mentioned before, I was dismissive - perhaps even rude. Sorry about that Egon, Sam. I still think it is better avoided; in the case of transparency, I don't know why alpha values can't be used for fill colours. I understand there was some SWT problems..?

Anyway, here is a sketch of a possible two-pass system, that would allow some of these adjustments to be made:

That's unreadably small in thumbnail - click for bigger, as usual. The basic idea would be to have one element tree with model-space values, and one with screen space values. I've made the distinction between double and integer, but Java2D will draw with doubles, so that is not important.

More chemical signature example

2009-06-23T14:58:00.006+01:00

A neat and tidy example of generating a structure from atoms and signatures occurred to me : adenine. It is known that this molecule forms readily from HCN under conditions similar to those thought to exist on the early Earth.

Anyway, the point is that the structure is almost exactly 6 C-N units (ignoring hydrogens and multiple bonds). A C-N is very similar to a height-1 signature, and the height-1 signatures for adenine are shown in the figure below:

Each signature (a-e) is mapped to an atom in the original structure, and then a table is shown of the maximum number of compatible bonds possible between each pair of signatures. Note that the table is not symmetric, as the compatibility operation is not commutative. Also, notice that the diagonal is not all 0 - bonds can form between atoms that both have b as a target signature.

Although these target signatures are smaller than the diameter of the adenine molecular graph, they still seem to define it completely. However, they are not guaranteed to produce only adenine as a solution - there may be other molecules with the same signature of (3a + 2b+3c+d+e).

Edit: Hmmm. That's not actually adenine, is it? In fact, it doesn't seem to be in ChEBI either. I probably should check my examples before posting them to the world.

Signature bond compatibility

2009-06-23T10:37:00.004+01:00

So, given my previous posts on what Faulon's signatures are, here is an explanation of how they are used in the structure enumeration algorithm that I am almost finished implementing.

The core test in this algorithm is for compatible bonds. Two atoms are only joined if : a) they have compatible target signatures and b) there are less than the target number of bonds already. A target signature here is just a signature that is set on the atom for it to match, like a pattern.

The first of these tests is illustrated here:

Another (overly) complex diagram! But the formula here is a bit difficult to interpret otherwise. In the top left corner is a graph G (slightly resembling hexane without the hydrogens) which is, by convention, composed of vertices (V) and edges (E).

The equation to the right of the graph defines part of the condition for a compatible bond. The tau terms are just target signatures, as shown on the upper right. The tricky term is ^h-1σ_τ(y)(z) which means 'a signature starting from the neighbour z of y in the subgraph defined by τ(y)'. This requires using a target signature (τ(y)) as if it was a subgraph - shown on the bottom left for the target b and the neighbour n.

The same process is repeated on the bottom right of the figure for b and m - which matches the height - 1 target signature for c. This should make sense, since the atoms labelled with b in the graph are attached to both a and c - so the signature for b must be compatible with both. It is easy to check that a and c are not compatible, and cannot therefore be bonded.

The Voynich blog?

2009-06-18T01:24:00.002+01:00

I apologise to any readers. I do tend to make complex diagrams. Perhaps some distant historian will find them in an archive and a new voynich-style mystery will result..

This does mean something, but you would have to read the previous post to discover what, exactly.

Faulon's Signatures : A Possible Interpretation

2009-06-17T23:11:00.012+01:00

Several recent papers by Faulon concern an idea he calls 'signatures'. This post is just a record of what I understood them to be.

Firstly, a signature is a subgraph of a molecular graph. There is a distinction between atomic signatures - which is a tree rooted at a particular atom - and a molecular signature, which is the set of atomic signatures for each atom in a molecule.

A tree is a graph with no cycles, so an atomic signature is not just a subgraph. Like a path, a signature has a length - or rather a height. Here is a picture of signatures of heights 1-4 for a fused ring structure:

The graph G on the left has one of its atoms labelled (a), and each of the trees in the center is a signature rooted at that atom. On the right, is the simple string form of the tree, as a nested list. I should point out that the signatures in these images may not be canonical, as I worked them out by hand (as I have not yet fully implemented the canonization algorithm).

Signatures of the same height may be different for the atoms in a molecule. At a height of zero, it is simply the atoms. A signature of height one is each atom, plus its neighbours. For G, above, there are two distinct height-1 signatures. For greater heights in G, there are more:

These are three subgraphs (S_G) of G, rooted at three different atoms (a, b, c). Each one corresponds to a signature tree, which also correspond to different signature strings (not shown). The trees have been given square nodes, instead of circular ones, just to make them look different. From the symmetry of G, it may be clear that the other atoms also have one of these same height-2 signatures.

Finally, there are some odd properties of the trees created from the subgraphs, that become noticable in height-3 signatures of G. As mentioned above, a tree cannot have cycles, so when the paths radiating out from the root atom meet on the same atom, it will appear in the tree twice. Further, when paths cross the same bond - at the same time - both atoms in the bond will appear in both orders across two layers:

The subgraph S_Gshows the former case, by putting two new atoms corresponding to the duplicate visit to the bridging atom in G. For the subgraph S_H of the pentagon H the whole of the last bond visited is duplicated, and the signature tree has a pair of duplicate bonds at the leaves. The tree construction process forbids duplication of bonds except in these two ways.

Automorphism groups and fragment graphs

2009-06-16T18:08:00.007+01:00

Structure generation involves not just graph theory, but group theory. Or, I should say, it does in some of the papers I have read. For example, in this paper by J.L.Faulon, there is the sentence:

"The two main steps are to compute the orbits of the automorphism group of G and to saturate all the atoms of a chosen orbit

which may well be incomprehensible to many readers, except if the reader is a mathematician.

I am no mathematician, but thanks to some books on groups, I now understand both what an automorphism group is and what an orbit is. On the other hand, I also believe that this definition of how the algorithm works is overly complex. A more simple term might just be "fragment sets" - as it is fairly clear, if not mathematically exact. So, for the fragment graph [CH₃, CH₃, CH₂, CH₂, CH, CH] the fragment set is [CH₃, CH₂, CH].

Anyway, here is a short analysis of the automorphism group of the fragment graph [CH₂, CH₂]. This first image shows the tiny group of permutations that swaps the two fragments:

The notation is taken from an excellent book called "Visual Group Theory" that is also associated with some software called group explorer on sourceforge. It might be quite general, I suppose (and I hope I'm using it right), but it shows the permutation that swaps the fragments as a circled s. This is an automorphism with respect to the edges - in other words, after the swap, there are still bonds between [1, 2], [2, 3], [4, 5], and [5-6].

Another part of the automorphism group is a 'flip' like:

which is a little more complex, but shows how 'flipping' each fragment separately combines to form four possible permutations. If this does not seem particularly tricky, consider what happens if you take the direct product of these two groups:

Assuming I have done it right, this should show most (all?) of the automorphisms of the fragment graph. It does look pretty cool, but I don't think that it gets me any closer to implementing the cursed algorithm :)

Simplest example : 1-methylcyclobutane

2009-05-23T11:05:00.003+01:00

This is the simplest possible hydrocarbon example (there being none with four carbons) of multiple isomorphic solutions:

The only difference in numbering is the swap of 2 and 4. Oh well.

Doesn't quite work

2009-05-23T00:25:00.002+01:00

Sadly, using the degree sequence as a kind of 'mask' on canonically numbered matrices isn't enough. I should have guessed, it was too simple to be true :(

The smiles for all matrices in the (2, 2, 2, 2, 2, 3, 3) partition are:

C1CC2CC(C1)C2

C1C2CC12.C1CC1

C1CC2CCC1C2

C1CC1CC2CC2

C1CCC2CC2(C1)

and (ignoring the third one, which is disconnected) there are only 3 or 4 different simple graphs in the list.

Perhaps I would have known all this if I knew what automorphisms were, or orbits of elements in sets, or any of that...

Bridged Cyclohexanes - same partition, different PVR sequences

2009-05-22T18:16:00.003+01:00

I wanted to find an example of two non-isomorphic molecules with the same degree sequence, but different PVR sequences. The simplest I could think of was this pair:

although there might be a similar pair of bridged 5-membered rings, I suppose.

The meaning of this is quite simple - for a particular degree sequence (partition), there are multiple (different) molecules with a valid PVR sequence. I thought that this must be true, but there are none in the G4 set.

Rows and Columns

2009-05-21T17:12:00.003+01:00

Aha! So, a comment by Egon (on my last post) showed the benefits of showing people what you do. He suggested summing the columns - not converting from binary to integers, as with the rows - to remove the last traces of redundancy. So this seems to work for graphs with 4 vertices:

Sorry for the extremely detailed diagram, but it is necessary to show my point. These matrix/graph pairs are all the PVR numbered adjacency matrices for n=4. There are isomorphic structures here, but note the column sums along the bottom of each matrix.

These column sums form another sequence - which can be used to select only one of the isomorphs. Arbitrarily, we choose sequences that are partially ordered - i.e. no number in the sequence is less than the previous number in the sequence. This seems to work...

PVR numbering scheme not solution to all woes : film at 11

2009-05-21T14:02:00.004+01:00

On a whim, I decided to try generating all adjacency matrices with the property that they are PVR numbered. The short summary of that link is that a matrix can be expressed as a sequence of positive integers by considering each row of the matrix as a binary number.

The point of doing this (I thought) was that you can number a molecule in such a way that the adjacency matrix is PVR-numbered, and that this is canonical. So my cunning plan was to generate all sequences of n numbers that are partially ordered, choosing them from [1, 2ⁿ] to give all non-redundant (simple) graphs with n vertices.

Unfortunately, it seems like this can't work:

This image shows all adjacency matrices for n = 3 which are PVR-numbered. They are made by backtracking through all sequences of integers with a partial order, pruning the solutions using the symmetry of the matrix as a constraint.

Anyway, the point is that the first two graphs are clearly isomorphic! More simply, they both represent propane. Maybe this is well known, but it's a surprise to me...

Generation : Overview

2009-05-18T17:07:00.004+01:00

To sum up the previous post flood; generation of constitutional isomers from the elemental formula can be done by generating all partitions of the total 'free' valence of the heavy atoms. The overall scheme is shown here:

(click for bigger, as usual). So, for each formula, multiple partitions can be made, and each of these makes multiple sub-partitions, and each of these correspond to one or more molecules.

Now, I won't pretend that any of this is particularly novel. I am no doubt re-expressing the problem of generating all possible molecules in a slightly different way. Having tried (and failed) to implement published methods, this was the best I could come up with.

I suspect that there are many improvements that could be made to the algorithm, and the implementation of it. Getting something that works, even in a limited way, seems like progress, however :)

Final step : nested partition thingies to actual molecules

2009-05-18T13:28:00.004+01:00

So, the last step:

One thing to note about this is that the algorithm again has to backtrack to get all the molecules for any sub-partition list (another name for the things like [[3, 1], [3], [1, 1], [1]]).

Anyway, on the right hand side is the final (only) molecule made for this nested partition. It has the [4, 3, 2, 1] partition structure, naturally, and the correct constitutional formula.

Now what would be nice, would be to combine the second and third steps, so that only those nested partitions that produce valid molecules were tried. However, as the saying goes : "First, make it work, then make it work fast".

Partitions to er..'nested partitions'

2009-05-18T12:48:00.003+01:00

So I don't have a good name for the objects that I create half-way through this process : the code uses 'solution', which is confusing. Anyway, here is the process:

The left hand side is clear enough, I think, and follows on from the image in the previous post. Conceptually, it is similar to 'gathering' the attachment points into half-bonds in all possible ways. So, the 4 attachment points on the bare carbon fragment can become a triple (half) bond, and a single half bond. This is shown, as [3, 1].

Of course, there are other possibilities, and each combination at each atom fragment has to be paired with each other possibility at each other fragment! If this sounds like a backtracking problem, then you might understand why I did exactly this in the code.

What would be nice, would be to prune the solutions - for example, [[3, 1], [2, 1], [1, 1], [1]] is generated, but is clearly impossible. Neither the triple half-bond, nor the double half-bond have partners. Pruning these would be easy enough; but looking more than one position ahead seems tricky.

Formula to Partitions

2009-05-18T11:51:00.004+01:00

For the benefit of my (2) readers, here is the process of going from the chemical formula to the partitions, and what this actually means.

So, the list of numbers at the bottom (the partition) is the simplest possible representation of the attachment points for each atom fragment.

Oh, and the python code for generating partitions for is here - it is basically a straight copy of the code from the book "Combinatorial Algorithms : Generation, Enumeration, and Search", which I highly recommend - a good balance of maths and computer science.

Uniqified

2009-05-17T21:01:00.003+01:00

Hmm. Not my favourite solution, but I think that the isomorphism checks can be done in batches, resulting in actual isomorph spaces, like this (for C5H10):

I should probably check that these are right...

Step one : catch them all

2009-05-17T17:11:00.003+01:00

Well, this is the C4H6 space again. Or <10,4> as I have taken to calling it, for no obvious reason.

There are duplicates, yes. But there are nine (unique) structures here, which is good.

Exploring the wild beasts of the layout jungle

2009-04-30T14:37:00.007+01:00

There was a bug submitted to the CDK sourceforge tracker (bug number 2783741) with a list of molecules that are laid out badly. I had a look at some of them with the help of bioclipse. For example, this calix[4]arene:

or this:

which is a clearer case of something going wrong. More difficult is structures that are fully 3D, like:

Can you guess what it is? :) Try the 3D version (also made with bioclipse, using the CDK 3D layout):

It's a paracyclophane! The phenyl rings are lost in the 2D layout because there is a bigger 'ring'. Perhaps a chemist would look at the 3D structure and think that those chains are linkers, not parts of a ring, but the algorithm doesn't know this.

I think that it is difficult to have general rules for this. Of course, any fully 3D structure will be difficult to lay out in 2D (if it is not embeddable in the plane then it is impossible) so things like this:

are truly awful.

Bioclipse : Safe When Used As Directed

2009-04-27T23:51:00.004+01:00

Finally used bioclipse for a real purpose, and to good effect, too:

what this shows (the images do get larger if you click on them! :) is the following basic workflow:

1) Exploring manager functions in the js console (bottom).

2) Writing a script in the js editor (top left).

3) Running and getting feedback in the rhino console (far right).

4) Viewing the results in the sdf viewer (top right).

What I was doing was searching through an sdf file (C10H16_filtered.sdf) for all structures with a cyclohexane ring as a substructure, then writing those out to a file.

Probably could have been done 5 other ways, but, well, it was more fun this way.

Oh, and it is a gist here.

ChEBI in Bioclipse

2009-04-23T12:08:00.003+01:00

Says it all, really.

The scrolling can hang if you move too fast, which might be due to garbage collection? It is a very large file....

Portable whiteboard : deployed

2009-04-17T01:07:00.002+01:00

Proposed CDK changes related to PDBReader and BioPolymer

2009-04-15T17:25:00.003+01:00

This is expanding on one of the points that Rajarshi made in his blog (which he followed up here) on the PDB file handling capabilities of the CDK. There are two related topics : reading of PDB format files (the ancient, fixed column-width ATOM files) and the model that these are read into.

PDBReader

The old PDB format is being replaced with mmCIF and/or PDBML formats. Only there are lots of programs that write out this format, so it makes sense to still support it for a while at least.

However, it is a quite nasty format, in some ways. Not so much the fixed column width, but the fact that crystallographers abuse the file format in all sorts of ways. Even simple things like expecting that atom numbers will always increase, may not be true.

So it is not easy making a good reader for PDB files. The current CDK one won't read a file with just ATOM records, for example. Think that's reasonable? Well, tough luck for people that made programs that produce simple files like this.

A more serious problem is the fact that you can't get properly connected ligands from a PDB file. Or easily get at the disordered regions. Or get at the waters. Well, sort of - I suppose that many of these things can be done after reading, with CDK classes.

BioPolymer

In some ways - so long as the atoms are read in - anything can be done to the model post-reading. However, the point of having data model classes like Polymer and BioPolymer is to capture some of the complexity of the macromolecule's organisation.

BioPolymer and the PDBReader do a good job of reading and storing the information in the header files (except that it is not always right!). Apart from calling chains 'Strands', some things are done reasonably well. I don't think that ligand atoms should be stored 'loose' as they are, but probably in referenced atom containers.

The real difficulty with modelling proteins lies with representing the hierarchy. One way is 'PMCA' - protein, model, chain, atom. This misses out secondary structure, but it may be too literal to have objects for every concept; the CDK stores the secondary structure given in the header file as IPDBStructure objects - with insertion codes, which is good to see.

However, there are more secondary structures than helix, turn, and strand. I'm never quite sure what the best way to model the more flexible situation, though.

Integration with biojava

It seems a shame that several open source java projects have very little in the way of integration (that I can see). At least on this topic. Of course, Egon (among others) has done work on both both CDK and Jmol, but it concerns me that incompatible ways of doing things lead to projects drifting apart, that should work together.

For example, here is a class I wrote today (for someone else) that uses both biojava and the CDK. It is very much a hack, but the key point is the method makeLigandsFromGroups that takes both a List of biojava.Group objects and a List of IMolecules, along the way converting biojava Atoms to CDK Atoms.

Clearly, biojava has a better interface to its model as you can get a List of the hetatm groups. The CDK, on the other hand, has better support for determining atom types and setting properties on them.

Generated Wallpaper, anyone?

2009-04-09T16:47:00.001+01:00

Heh.

This is exhaustive generation of all C4Hn (where n=any number of hydrogens), badly smashed together onto one canvas.

Structure zoo

2009-04-08T17:43:00.004+01:00

So, I'm getting better at generating structures :
but, not quite there yet. Oh, and the [1.1.1] propellane is drawn oddly (with one of its carbons at 0,0) due to a mac-specific bug in the layout. Irritating, but difficult to fix.

edit: Oh, and for interests sake, here is an auto-generated tree of graphs

Even more complex looking is this! :

which is a set of structures produced by descending to depth 6 for 5-carbon graphs.

Numbering atoms, numbering vertices

2009-04-03T14:42:00.003+01:00

Further to the similarities between numbering atoms in a structure, and generating unique graphs here is this:

which shows the same molecule with two different numberings on the left, and the resulting graphs on the right. The double bond is not shown on the graphs; but it would probably have to be a labelling of the edge, rather than an actual multiple edge, to still be a simple graph.

So, this quickly shows how - if you start with the vertices of the graph and connect 'all possible ways' - you get molecules that are isomorphic, but numbered differently. Therefore (perhaps) the numbering of the vertices and edges is one of he keys to not creating all the isomorphs and then having to expensively check them all.

Step 3 : Extend all possible ways

2009-04-03T14:05:00.003+01:00

The title of this post refers to the tendency of algorithms in papers to have detailed explanation of every step except the most crucial one. Right now I am rediscovering this peculiar pleasure in structure generation.

As an example - or more as visual decoration - have this image:

which looks nice, but needs some explanation. The diagrams in boxes that look like parachutes (as one of my colleagues put it :) are simple representations of atoms connected by bonds. Each point is an atom, and a curved line connecting them is a bond.

These are grouped together by what structure they correspond to; which is shown on the right of each set of diagrams. What this shows, then, is the redundancy you get from a simple generator. If you connect all atoms like this:

  for atomA in atoms:
  for atomB in atoms greater than atomA:
    connect(atomA, atomB)

You quickly get a very large number of isomorphic structures. And then your process runs out of memory, in my experience.

Oh, and the numbers below each graph are the sum of the indegree/outdegree (the degree?) of the that vertex - or the number of bonds the atom is in, in other words. This seems similar to the ideas behind Morgan numbers, which I now understand a bit better.

Symmetric Generations

2009-03-12T14:50:00.003Z

I've been trying to work on an implementation of a structure generator algorithm due to Faulon (its a paper in this list somewhere). One problem I have difficulty with is reduction of the number of isomers generated. For example:

It may be hard to read, but the idea is that the full tree of possible ways to attach {Br, Cl, F, I} to c1ccc1 is 4 * 3 * 2 * 1 (you can attach Bromine to all four of the carbons, leaving 3 places to attach chlorine, and so on).

Of course c1(Br)c(Cl)cc1 is the same as c1(Cl)c(Br)cc1 [never mind that C1(Br)=C(Cl)C=C1 is not the same as C1=C(Br)C(Cl)=C1]. Or, in other words, there is a high degree of symmetry in the tree.

There is a way to solve this, perhaps even described in the paper - if I could just understand them...

Multi-view

2009-02-25T18:43:00.005Z

A quick example of a mini-application made using the new JChempaint renderer:

the picture is of an isomer space (C₃H₇NO), with compact mode on and atoms rendered as circles. The code is here, for now:

http://gist.github.com/70342

which is a kind of dumb way to achieve this, as it creates an instance of a renderer for each molecule, instead of adding them to a molecule set, and then laying that out...

edit: just realised; it's rendering the hydrogens as compact black circles :(

edit2: Ahhh. that's better:

edit3: Ha! This was a mistake, but it looks kind of cool:

lone pair rendering

2009-02-14T19:06:00.007Z

JCP now has minimal lone pair display. I would prefer the layout to be at the corners of a square, rather than on the edges. Videlicet, they are currently only at N, W, S, E; I think that NW, SW, SW, and NE would be better.

Strangely lone pairs don't seem to appear in CML files when written out, but Egon says he will look at this. Perhaps I should file a bug report...

Oh, and radicals are implemented too, ~~but I don't have a picture of that~~ (see right). They are in different generators, but I guess a single 'DotGenerator' could do both :)

Spectral Full House

2009-01-29T20:40:00.005Z

So, all of the isomers of C4H11O are in NMRShiftDB and here are all the experimental and predicted carbon spectra:

It's not obvious from this picture, but not all of the predicted spectra are unique matches for their experimental partners. In other words, you could not pick out the right molecule by comparing the predicted and experimental spectra.

The situation is more difficult still for larger isomer spaces, where the predicted spectra may be exactly the same for sub-sets of the isomers. There are still many with unique predictions, but the rest follow a sort of power-law distribution of spectral-equivalent sets.

EDIT: As per a suggestion by egon, here is a table of top hits (a yellow square indicates the top match):

C4H11N network, labelled with pictures

2009-01-16T18:38:00.004Z

A really tiny isomer set, with a particularly regular network:

Nicely enough, all 8 of these structures are in NMRShiftDB, so it will be possible to compare all against all.

One with smiles strings on

2009-01-14T15:12:00.004Z

Slightly more useful/comprehensible.
Also, this is a more interesting isomer space (C3H7ON) which is more fragmented at a 50% similarity cutoff.

just one more

2009-01-12T20:12:00.003Z

C8H16 this time. The central bridge has split into two parts, it seems. The clustering is somewhat artificial, I suppose...
An example from one of these is C=C(C(C)C)C(C)C, and from the other is C=CC(CC)CCC. Hmmm.

another network

2009-01-12T19:30:00.003Z

Another one (C7H14 this time):
It's a bit cumbersome, but I managed to get smiles strings to show as tooltip text. This tells me that the three vertices in the center of this picture (the bridging ones) are C=C(C)CCCC, C=C(CC)C(C)C, and C=C(CC)CCC. Not sure why.

Oh, and the 'lobe' on the left seem to be (all?) non-cyclic, while the ones on the right seem to be cyclic, which makes sense.

isomer networks

2009-01-12T17:05:00.004Z

hmmm. Been a month. Oh well. Here is a picure of an 'isomer network':
for the isomers of C6H12. I generated them using molgen, compared the structures all-v-all with the CDK, and visualized the connectivity graph with the help of JUNG. So hardly any work actually by me...

Each vertex in the graph is a structure, and each edge is a tanimoto similarity between fingerprints of greater than 0.5 - fairly arbitrary, but I just wanted to see if it worked. The next step is to use predicted spectrum similarity instead of molecular-fingerprint similarity.

Why you can't make custom progress monitors in eclipse

2008-12-09T17:47:00.003Z

I wanted to make a custom progress monitor with a temperature graph in Eclipse, but it does not seem to be possible. Here's why:

There is a ProgressProvider class that has an abstract createMonitor(Job job) method that is called in the JobManager when creating jobs. The standard concrete implementing class is the ProgressManager, shown above in a dashed box labelled "Discouraged Access". So you can't just override this class.

The JobMonitor is an inner class of ProgressManager, and it is not possible to provide an alternative. Once it has been created and passed to the Job, it has already been displayed as a Dialog, so it can't be wrapped with another Dialog.

Seems impossible, really.

Intermediate vs direct rendering

2008-12-07T17:56:00.003Z

Just a small, informal diagram of how the new rendering system works, to balance out all the speculative UML sketches I made...
I have not shown the elements inside the element groups, but they are things like line elements, oval elements, and so on.

Creating new molecule editing programs with cdk

2008-12-05T14:24:00.004Z

After all the work that has gone into making a flexible core for JChemPaint (the controller and renderer modules) I finally have the opportunity to use it for what I wanted : a small custom program for editing molecules from scratch and displaying the spectrum as you go.
This may not look like much, but compare to this next screenshot:
Obviously this is a fairly unlikely molecule, but (if you believe me :), the new ring has caused the program to re-predict a spectrum (using NMRShiftDb). The point of this was to teach me more about the connection between molecules and spectra.

CDK-City

2008-11-25T17:51:00.003Z

On the mean streets of CDK-city, one massive tower block dominates the skyline in the renderer district - the Renderer2DModel tower...

Masak pointed me to this site due to Richard Wettel which is about a static-analysis tool for making so called 'code cities'.

Basically, the parking lots (like SMARTSParserConstants) have lots of attributes (in fact, this class only has attributes) while skyscrapers (like in the debug package) have many more methods than attributes (most of the debug classes have no attributes except inherited ones).

There's an eclipse plugin from here called MooseBrewer, and I think there are other tools there I haven't explored. It's great to have a map to this busy code metropolis, as I wondered what my workplace looked like :)

chalky

2008-11-23T15:33:00.004Z

I wanted to show something that hints at the things that the new architecture can afford us:
This is using a Java2D graphics Paint object to make it look like chalk...kindof. It's a very simplistic way of doing it by making a small image with a random number of white, gray, lightgray, and black pixels.

edit: it doesn't look so good at small scales

some tweaking of stroke widths and so on is essential.

Seeing Double?

2008-11-20T23:27:00.005Z

Comes from too much screen-time:
although refactoring often seems like running on the spot (you get nowhere fast), things happen behind the scenes..
A short code snippet.

More Annealing

2008-11-08T10:59:00.004Z

Some more detail here:

Here is a larger molecule, cembrane (pubchem link) which has 20 carbons. The run is longer, too, with 2,000 steps. It finds a spectrum with 100% match within 300 steps, in this case.

An even larger example is lanosterol (pubchem link) which has 30 carbons. This actually seems to be too large for C-NMR prediction using NMRShiftDB. It doesn't get the answer within 2,000 steps, and does not look like it is on course to do so:

The highlighted step (1928) is shown as a molecule in the box marked "final", but the score graph has levelled off by about the 400th step.

Adaptive Annealing Engine Test Screenshot

2008-11-07T14:09:00.003Z

So, finally, some work that I am meant to be doing (click for bigger, as usual):

It's crude, but it is starting to work. The screenshot shows only the first 100 steps of a run, but clearly you can get the right answer (almost by chance, actually) for such a small molecule when you are essentially just permuting the atoms.

Now to test more fully...

Shadows

2008-10-23T23:10:00.004+01:00

Not the most important use of rendering code..

:)

Symbols and Elements

2008-10-20T06:39:00.004+01:00

I have made attempts over the years to draw triangles and circles like, for example, this:

(in fact, this was drawn by David Westhead's code). There is a small design decision to be made between representing drawing elements as basic geometry or as symbols. A symbol may be the same as a geometric element (like a circle), but it might be a combination of them. Consider the 'N' and 'C' in boxes, above. Or this image:

where the mass number and charge are separate text elements on the left, and part of the whole symbol on the right. Neither is a 'better' way of doing things; they each have their advantages - and disadvantages.

It can be a lot clearer to use symbols, as each model object (atom, bond, helix, gene, cell) has a corresponding representation, and then a diagram composes the symbols into a manageable whole. On the other hand, you can re-use elements in different combinations for diagrams. A general vector drawing package would have Line, Text, Rectangle, Triangle, and so on.

So my vote would be to use symbols for the CDK/JCP renderer. The JCP application is not intended to be a full vector drawing package, but a specialised chemical editor.

The Trouble with Tribble Visitors

2008-10-17T16:43:00.005+01:00

So, I'm now partly sold on the power of a Visitor approach to rendering. Consider this snippet:


          if (clicked) diagram.accept(new DropShadowVisitor(g, 5, -5);
          else         diagram.accept(new DrawVisitor(g));

What this is doing is drawing a diagram normally unless the mouse is clicked, when it draws with a drop shadow. I see that the beauty of this is the ability to manipulate functionality as a block (just as in languages where you can pass around functions...).

However, I should point out that the approach has its tricky rapids as well as such smooth sailing. The image below is a spot-the-difference (click for bigger):

On the left is the version drawn by a naive first try at the drop visitor. Its methods look like this, the visitText(Text text) method:

    g.setColor(Color.LIGHT_GRAY);
   g.drawString(text.text, text.x, text.y);
   g.setColor(Color.BLACK);
   g.drawString(text.text, text.x + dx, text.y + dy);

The problem with this code is subtle - the elements are visited once each, and the shadow is rendered at the same time as the element. The change that had to be made was to have a boolean "drawingShadow" and to visit the elements twice

        this.drawingShadow = true;
   for (DiagramElement element : diagram.children) { element.accept(this); }
   this.drawingShadow = false;
   for (DiagramElement element : diagram.children) { element.accept(this); }

The shadow has to be drawn first - BUT first for the whole diagram, not just first for each element. So the new text method is :


        if (this.drawingShadow) {
   g.setColor(Color.LIGHT_GRAY);
   g.drawString(text.text, text.x, text.y);
  } else {
   g.setColor(Color.BLACK);
   g.drawString(text.text, text.x + dx, text.y + dy);
  }

Arvid's Renderer Design

2008-10-16T17:10:00.005+01:00

This is a sketch of my understanding of Arvid's renderer design:

I say 'my understanding' with good reason - a diagram can be biased, even a formal(ish) one like this, so I don't claim that this is definitive!

It is interesting, anyway, as it seems to be a combination of a Composite and a Visitor pattern. The RenderingModel implements IRenderingElement (and Iterable of IRenderingElements), which is Composite.

The Modules are Elements in the Visitor pattern, and the Elements are Visitors.

edit: What nonsense I talk! The IRenderingElements are Elements and the IRenderingModules are Visitors. That's better.

Font Management and GlyphVectors

2008-10-15T23:09:00.004+01:00

Fonts and text are always a pain when drawing vector graphics!

After hitting myself over the head over a stupid bug in my zoom tests (like this:)

JButton inButton = new JButton("IN");
inButton.setActionCommand("IN");
JButton outButton = new JButton("OUT");
inButton.setActionCommand("OUT");

I got center scaling implemented in my branch. Which lets me test a different approach to managing fonts. I had tried to be clever and do something like:

GlyphVector glyphs = font.createGlyphVector("N");
ArrayList shapes = new ArrayList();
for (int i = 0; i < glyphs.getNumberOfGlyphs(); i++) {
shapes.add(affineTransform.getTransformedShape(glyphs.get(i)));
}

but letter shapes fill horribly with lots of missing pixels. Anyway, I eventually settled on just storing the Glyphs themselves. Less pure, but it works, and it allows the size of TextSymbols to be computed and used by the class in between drawing.

So, for fonts there had to be a way in my architecture to change the font size when the whole molecule is scaled. The nicer approach is the obvious one - just to transform the Graphics object. However, it turns out there are advantages to the more cumbersome approach. Here's an image:

which shows rings at various sizes, with the numbers at their centers showing the font size at that scale. An important one is the 'size' (9-5) - which is 4, but there is no readable size of font at that size. However, the FontManager class keeps track of how far below (or above) the minimum and maximum font sizes, and then returns to that size at the appropriate point.

A correction, and a simpler example

2008-10-14T15:47:00.005+01:00

There was of course a deliberate mistake in the previous post:

as the highlighting shows, the matrix for cubane numbered PVR style is neatly divided. A simpler example is cyclobutane:

which also shows the important point that, even if the sum of the row-numbers is the same as for other possible numberings, the PVR numbering is still the only ordered one.

On Canonical Numberings

2008-10-12T18:18:00.007+01:00

So, after reading* this (2005) paper : "On Canonical Numbering of Carbon Atoms in Fullerenes : C60 Buckminsterfullerene" (link) I made some pictures to illustrate the difference between it and the numbering scheme used for SMILES (as described here). Er, which is used in the CDK.

Anyway, the point is that the scheme used by Plavšić, Vukičević, and Randić (or PVR as I will refer to them, I hope they don't mind!) numbers the atoms in a way that produces an adjacency matrix with a particular property. If you consider the rows of the matrix to be binary numbers, then the set of numbers is the smallest possible. So, for example:

The structure on the left is cubane, with its adjacency matrix on the right. The column on the far right shows the rows of the matrix in base 10. They are clearly in order. Now what happens for the SMILES? Well:
Here, the rows are neither in order (I'm not sure from their paper whether the ordering is an expected outcome for all structures, nor have I checked...) nor is their sum less than for PVR scheme - 765 vs 753.

Of course, the PVR labelling would be useless for generating SMILES for cubane since there is no way to get a path from it. Indeed, the labels are designed to be maximally unfriendly by pairing the highest with the lowest.

Furthermore their scheme goes on to label bonds and rings:

Which also look quite random; or, as they say :

"... we admit that the final labels ... do not appear »orderly« but one has to recognise that there is no »simple« labelling in [fullerenes] that will appear simple"

which makes sense. Obviously, not for cubane here, but for C60/C70 and so on it does.

*(probably because the name follows the "On X" paper naming scheme)

KDTree in SymbolTree

2008-10-10T17:59:00.002+01:00

So, the new Symbol tree base class now has a working implementation of a KDTree. (Or as close as the test can tell). What this means is that, should you ever need to edit (not just display) a very large 2D structure, you would do something like:

SymbolTree tree = createTree(); // get a tree somehow
tree.setHitDistance(minDistance); // the minimum distance from mouse to symbol
tree.useKDTree(true); // essential
tree.highlightClosestSymbol(point); // will now do a fast search for the closest

Of course, this is really only useful for a) large trees and b) when doing many "closest symbol" operations. For example, highlighting when moving the mouse.

The best situation would be to dynamically call the "useKDTree" method when the size was above some threshold (100 atoms?). The interface is the same, either way.

Steps to move Bioclipse to the new CDK plugins

2008-10-08T12:06:00.003+01:00

Thanks to egon the monolithic cdk plugin is now about 40 small plugins. Here's how to migrate...

In no particular order:

Step 1: Delete plugin org.openscience.cdk from the workspace
Step 2: Download the third-party library plugins from bioclipse2/trunk/cdk-externals/
Step 3: Download the (many) new cdk plugins from bioclipse2/trunk/cdk-externals/trunk/
Step 4: Update those plugins that complain about missing org.openscience.cdk dependency
(Step 5: Possibly change the Require-bundle for org.apache.log4j : as in this bug report)

It should now all compile...