Thursday, March 20, 2014

Eclipse compiler produces faster FP code?

While testing this very simple FP code in Java:

package the.plateisbad;

public class Simple {

  public static void main(String[] args) {
    final long arg_flops;
    double i,x = 0.0d, y = 0.0d;
    final long start, end;

    if (args.length != 2) {
      System.out.println("Usage: the.plateisbad.Simple ");
      return;
    }

    arg_flops = Long.parseLong(args[0]);
    y = Long.parseLong(args[1]);

    System.out.println("Thinking really hard for " + arg_flops + " flops...");

    start = System.currentTimeMillis();

    for (i = 0; i < arg_flops; i++) {
      x = i * y;
    }

    end = System.currentTimeMillis();

    System.out.println("We calculated: " + x + " in " +(end-start)+ " ms");
  }
}



I've stumbled over the fact, that it runs considerably faster when compiled with the Eclipse ECJ compiler compared to a standard javac.

With ECJ, executed with JDK 1.7:

java -server the.plateisbad.Simple 1000000000 3

Thinking really hard for 1000000000 flops...
We calculated: 2.999999997E9 in 1964 ms

With javac, executed with JDK 1.7:

java -server the.plateisbad.Simple 1000000000 3

Thinking really hard for 1000000000 flops...
We calculated: 2.999999997E9 in 3514 ms

With the new JDK 1.8, there is no noticeable difference between javac and ECJ:

java -server the.plateisbad.Simple 1000000000 3

Thinking really hard for 1000000000 flops...
We calculated: 2.999999997E9 in 3727 ms

but it is always the slowest of the three. The Bytecode tells me that ECJ builds a tail controlled loop which loops while i is < arg_flops:

  64: invokestatic  #52                 // Method java/lang/System.currentTimeMillis:()J
      67: lstore        9
      69: dconst_0     
      70: dstore_3     
      71: goto          84
      74: dload_3      
      75: dload         7
      77: dmul         
      78: dstore        5
      80: dload_3      
      81: dconst_1     
      82: dadd         
      83: dstore_3     
      84: dload_3      
      85: lload_1      
      86: l2d          
      87: dcmpg        
      88: iflt          74
      91: invokestatic  #52                 // Method java/lang/System.currentTimeMillis:()J


while javac builds a head controlled loop that exits if i >= arg_flops:

      67: invokestatic  #13                 // Method java/lang/System.currentTimeMillis:()J
      70: lstore        9
      72: dconst_0     
      73: dstore_3     
      74: dload_3      
      75: lload_1      
      76: l2d          
      77: dcmpg        
      78: ifge          94
      81: dload_3      
      82: dload         7
      84: dmul         
      85: dstore        5
      87: dload_3      
      88: dconst_1     
      89: dadd         
      90: dstore_3     
      91: goto          74
      94: invokestatic  #13                 // Method java/lang/System.currentTimeMillis:()J

And ECJ uses StringBuffer while javac uses StringBuilder for the String operations, but since these are not in the loop, that should not make any difference.

Does somebody know what is going on here?

UPDATE: This seems to be an anomaly. SciMark 2.0 shows now significant differences between ECJ and javac and jdk1.7 and jdk1.8 - with 1.8 being slightly faster.

Wednesday, March 5, 2014

A suggestion to all architects of high-security buildings ;->

Please, don't put the restrooms outside the security gates!

Wednesday, February 12, 2014

Arbitrary parallel (well, almost) ad-hoc queries with PostgreSQL

Contemporary PostgreSQL lacks the ability to run single queries on multiple cores, nodes etc., i.e. it lacks automatic horizontal scaling. While this seems to be under development, what can be done today?

PL/Proxy allows database partitioning and RUN ON ALL executes the function on all nodes simultaneously. PL/Proxy is limited to the partitioned execution of functions and has good reasons for this design. But PostgreSQL can execute dynamic SQL within functions, so let's see how far we can get.

Worker function (on all worker nodes):

CREATE OR REPLACE FUNCTION parallel_query(statement text)
  RETURNS SETOF record AS
$BODY$
DECLARE r record;
BEGIN
IF lower($1) LIKE 'select%' THEN
FOR r IN EXECUTE $1 LOOP
RETURN NEXT r;
END LOOP;
ELSE
RAISE EXCEPTION 'Only queries allowed';
END IF;
END
$BODY$
  LANGUAGE plpgsql VOLATILE;


Proxy function (on all head nodes):

CREATE OR REPLACE FUNCTION parallel_query(statement text)
  RETURNS SETOF record AS
$BODY$
 CLUSTER 'head'; RUN ON ALL;
$BODY$
  LANGUAGE plproxy VOLATILE;


Table (on all worker nodes):

CREATE TABLE users
(
  username text NOT NULL,
  CONSTRAINT users_pkey PRIMARY KEY (username)
)
WITH (
  OIDS=FALSE
);


With 10000 rows in two nodes, partitioned by username hash (~5000 on each node)

select * from parallel_query('select * from users') as (username text);

returns all 10000 rows. Since the nodes can be databases within the same server, there is no need for additional hardware, server installations etc. But if more performance is required in the future, adding more boxes is possible.

All it takes is logical partitioning and a bit of PL/pgSQL if you really need to run parallel queries.

There are some differences though. Take the following query:

select * from execute_query('select max(username) from users') as (username text);

"user_name_9995"
"user_name_9999"

It now returns two maximums, one for each partition. To get the expected result a second stage is needed:

select max(username) from execute_query('select max(username) from users') as (username text);

"user_name_9999"

The same applies for other aggregation functions like avg() etc.

The proxy function can finally be hidden in a VIEW:

CREATE OR REPLACE VIEW "users" AS select * from parallel_query('select * from users') as (username text);

Thursday, January 23, 2014

From palloc() to palloc0()

I've just replaced all palloc()/memset() pairs for getting zeroed-out memory with palloc0(). I doubt that this will show a significant speedup, but it fits better into the PostgreSQL memory model and reduces the lines of code.

E.g.:

retval = (text *) palloc (len + VARHDRSZ);
memset(retval,0x0,len + VARHDRSZ);

becomes:

retval = (text *) palloc0 (len + VARHDRSZ);

The changes have been comitted to Github.

Tuesday, December 17, 2013

It compiles...

>But, since Postgresql 9.3 just came out, I'll see if I can manage at least to compile pgchem::tigress against 9.3, OpenBabel 2.3.2 and Indigo 1.1.11 before 2013 ends...

Short update: I have merged the two pull requests from Steffen Neumann and Björn Grüning into the repository and then made some minor corrections. LANGUAGE C (case insensitive but without the '') solves the one problem, replacing int4 with int32 the other. Actually this was no gcc issue, but the PostgreSQL guys have apparently decided to remove int4 as an PostgreSQL internal datatype from 9.2 to 9.3.

I had always wondered why there were so many duplicate but internally identical (e.g. int4 (number of bytes) and int32 (number of bits)) datatypes anyway. Probably some legacy...

Also, pgchem_tigress can now be installed as an relocatable extension with the CREATE EXTENSION mechanism. No need to run various installation scripts anymore (but still supported).

I still need to put an 'install' target into the Makefile to fully integrate pgchem_tigress into the extension system. Currently, the shared objects still have to be copied by hand before calling CREATE EXTENSION The PGXS Makefile already takes care of that - there is progress towards a new release.

Thursday, October 24, 2013

Natural vs. Chemical

Last week I attended a very interesting lecture about metabolic pathways.

Maybe it would help all the people who like to make a firm distinction between 'natural products' -> good and 'chemical products' -> bad to attend such a lecture from time to time...

Wednesday, September 18, 2013

Long time no see

Hi all,

the development of pgchem::tigress - and this Blog - have been quiet for some time now.

Professionally, I'm currently moving into two new spheres of activity, document management for product registration (i.e. license to sell) and bioinformatics, and this occupies most of my time. That my wife had to go to the hospital for a few months (nothing life threatening, but quite hindering) didn't help much either. ;->

But, since Postgresql 9.3 just came out, I'll see if I can manage at least to compile pgchem::tigress against 9.3, OpenBabel 2.3.2 and Indigo 1.1.11 before 2013 ends...

Wednesday, September 26, 2012

Bekloppt

Drug discovery as a Facebook game? WTF...

Friday, December 23, 2011

Christmas presents

Jérôme Pansanel has completed the new serialization code for mychem and pgchem, so there is no need to handle stereo and non-stereo queries differently anymore. I have moved the index functions to GCC's vector extensions where applicable, and the first result is that index build times have been roughly cut by half while substructure search times have improved, but not that much.

Index build times


SystemIndex build time
pgchem with OpenBabel or Indigo352137 ms
pgchem with OpenBabel or Indigo vectorized192815 ms

OpenBabel with binary storage and FP2 fingerprint vectorized


QueryHitsno IndexHitswith Index
GH2484098416 ms484017044 ms
GH726094053 ms2601564 ms
GH13580113690 ms58034504 ms
GH162691099365 ms2691055154 ms

Merry Christmas and a happy new year!

Saturday, December 17, 2011

Benchmark data published

The base data used in the previous benchmarks can be found here: pubchem_0_100000.

Thursday, December 15, 2011

Selected GH17 results for 10^6 structures

GH17 substructure search speed


OpenBabel with binary+SMILES storage and FP2 fingerprint


QueryHitsno IndexHitswith Index
GH24840108807 ms484021164 ms
GH7260105050 ms2601934 ms
GH13580118978 ms58052416 ms
GH1626910109886 ms2691064742 ms

Indigo with binary storage and ext+sub fingerprint


QueryHitsno IndexHitswith Index
GH24840213075 ms484027887 ms
GH7410178963 ms4104451 ms
GH13580251938 ms58039134 ms
GH1627100172534 ms2710080523 ms

Bingo 1.7beta2 with molfiles as text storage


QueryHitsno IndexHitswith Index
GH24710647889 ms471021733 ms
GH7410538784 ms4106658 ms
GH13580675093 ms58012418 ms
GH1627100528891 ms2710028541 ms

Index build times


SystemIndex build time
pgchem with OpenBabel or Indigo352137 ms
Bingo3458681ms


Again, Bingo without it's index is apparently killed by the overhead of parsing text into the internal molecule format. With index it's a mixed bag, while it shines at GH13 and GH16, pgchem is about equal or faster at GH2 and GH7.

Wednesday, December 14, 2011

GH17 results

Mikhail Rybalkin from GGA Software asked me to do this, so here it is...

The GH17 test queries used


There is a small set of queries used in article Chemical substructure search in SQL by Golovin and Henrick. These queries were lately reused in other articles:

GH1 ONC1CC(C(O)C1O)[n]2cnc3c(NC4CC4)ncnc23
GH2 Nc1ncnc2[n]cnc12
GH3 CNc1ncnc2[n](C)cnc12
GH4 Nc1ncnc2[n](cnc12)C3CCCC3
GH5 CC12CCC3C(CCC4=CC(O)CCC34C)C1CCC2
GH6 OC2=CC(=O)c1c(cccc1)O2
GH7 Nc1nnc(S)s1
GH8 C1C2SCCN2C1
GH9 CP(O)(O)=O
GH10 CCCCCP(O)(O)=O
GH11 N2CCC13CCCCC1C2Cc4c3cccc4
GH12 s1cncc1
GH13 C34CCC1C(CCC2CC(=O)CCC12)C3CCC4
GH14 CCCCCCCCCCCP(O)(O)=O
GH15 CC1CCCC1
GH16 CCC1CCCC1
GH17 CCCC1CCCC1

GH17 substructure search speed


OpenBabel with binary+SMILES storage and FP2 fingerprint


QueryHitsno IndexHitswith Index
GH109517 ms025 ms
GH24848519 ms484111 ms
GH3638632 ms6343 ms
GH458950 ms548 ms
GH53610020 ms3678 ms
GH608696 ms032 ms
GH7268279 ms2631 ms
GH81708454 ms17056 ms
GH93488068 ms34871 ms
GH10368820 ms3621 ms
GH11669113 ms6652 ms
GH128317920 ms831124 ms
GH13589864 ms58448 ms
GH1449998 ms436 ms
GH1530088549 ms3008555 ms
GH1626918665 ms2691501 ms
GH1722908717 ms2290560 ms

Indigo with binary storage and ext+sub fingerprint


QueryHitsno IndexHitswith Index
GH1028161 ms032 ms
GH248419188 ms484187 ms
GH36820498 ms6878 ms
GH4523388 ms533 ms
GH53625887 ms3662 ms
GH6021034 ms031 ms
GH74116302 ms4131 ms
GH817016629 ms17078 ms
GH937314005 ms37398 ms
GH103716881 ms3721 ms
GH116624091 ms6678 ms
GH1282914842 ms829210 ms
GH135824470 ms58133 ms
GH14421403 ms436 ms
GH15304715817 ms3047749 ms
GH16271016767 ms2710732 ms
GH17230417524 ms2304788 ms

Bingo 1.7beta2 with molfiles as text storage


QueryHitsno IndexHitswith Index
GH1070277 ms0125 ms
GH247156821 ms471156 ms
GH36857754 ms68125 ms
GH4560067 ms5125 ms
GH53665586 ms36140 ms
GH67957188 ms79125 ms
GH74152134 ms41125 ms
GH817051685 ms170140 ms
GH937347613 ms373138 ms
GH103749961 ms37110 ms
GH116661176 ms66125 ms
GH1277450281 ms829156 ms
GH135861108 ms58156 ms
GH14453636 ms4140 ms
GH15304750213 ms3047343 ms
GH16271051227 ms2710327 ms
GH17230451495 ms2304362 ms

Fingerprint efficiency (with regard to false positives)


FP2


QueryCandidates screenedHits matchedfalse positivesEfficiency
GH10001.000
GH248548410.998
GH3696360.913
GH4375320.135
GH512036840.300
GH6790790.000
GH74126150.634
GH817717070.960
GH9377348290.923
GH10373610.973
GH111236610.537
GH1283183101.000
GH1313465812880.043
GH14204160.200
GH15376030087520.800
GH16330526916140.814
GH17330522907150.762

ext+sub


QueryCandidates screenedHits matchedfalse positivesEfficiency
GH10001.000
GH248448411.000
GH3686801.000
GH45501.000
GH54736110.766
GH60001.000
GH7414101.000
GH817017001.000
GH937337301.000
GH10373701.000
GH11666601.000
GH1282982901.000
GH13259582010.224
GH14204160.200
GH1530613047140.995
GH1627202710100.996
GH17272023044160.847

Index build times


SystemIndex build time
pgchem with OpenBabel or Indigo25690 ms
Bingo336319 ms


Indigo's ext+sub fingerprint is truly more selective than FP2. Still, OpenBabel with binary storage shows the better prformance because of its faster matcher lower query overhead.

Also, the result for GH3, GH7, GH9, GH12, GH15, GH16, and GH17 are different between OpenBabel and Indigo, and Bingo finds 79 hits for GH6 where pgchem finds zero.

Index building on pgchem is 13 times faster than Bingo, but since pgchem (currently) does not support features like tautomer searching or SMARTS searching with index support this comparison is a bit like apples and oranges.

The slow performance of Bingo without index, comparable to pgchem without binary storage, is quite likely a result of the storage of molecules in textual representation. Parsing text to binary molecules is a first class performance killer. Unfortunately, there is no way to convert molecules into native format directly with Bingo for PostgreSQL, but Bingo does the conversion implicitly when building the index.

Friday, December 9, 2011

OBMol de-/serialization revisited

The current serialization/deserialization mechanism in pgchem and mychem does not preserve stereochemistry of OBMol objects properly. As Tim Vandermeersch wrote:

The OBChiralData isn't used anymore. Also the functions OBAtom::IsClockwise, and OBAtom::IsAntiClockwise are obsolate. Instead, you should serialize the OBCisTransStereo and OBTetrahedralStereo data objects.

Fortunately (for me, since I don't have the time at the moment), Jérôme Pansanel has:

...started the serialization of the OBCisTransStereo and OBTetrahedralStereo objects.

For the time being, I have tried the following workaround and it seems to be working well. First, I have removed now unneccessary code from the unserialization:

bool unserializeOBMol(OBBase* pOb, const char *serializedInput)
{
  OBMol* pmol = pOb->CastAndClear<OBMol>();
  OBMol &mol = *pmol;
  unsigned int i,natoms,nbonds;

  unsigned int *intptr = (unsigned int*) serializedInput;

  ++intptr;

  natoms = *intptr;

  ++intptr;

  nbonds = *intptr;

  ++intptr;

  _ATOM *atomptr = (_ATOM*) intptr;

  mol.ReserveAtoms(natoms);

  OBAtom atom;
  int stereo;

  for (i = 1; i <= natoms; i++) {
    atom.SetIdx(atomptr->idx);
    atom.SetHyb(atomptr->hybridization);
    atom.SetAtomicNum((int) atomptr->atomicnum);
    atom.SetIsotope((unsigned int) atomptr->isotope);
    atom.SetFormalCharge((int) atomptr->formalcharge);
    stereo = atomptr->stereo;

    if(stereo == 3) {
      atom.SetChiral();
    }

    atom.SetSpinMultiplicity((short) atomptr->spinmultiplicity);

    if(atomptr->aromatic != 0) {
      atom.SetAromatic();
    }

    if (!mol.AddAtom(atom)) {
      return false;
    }

    atom.Clear();

    ++atomptr;
  }

  _BOND *bondptr = (_BOND*) atomptr;

  unsigned int start,end,order,flags;

  for (i = 0;i < nbonds;i++) {
    flags = 0;

    start = bondptr->beginidx;
    end = bondptr->endidx;
    order = (int) bondptr->order;

    if (start == 0 || end == 0 || order == 0 || start > natoms || end > natoms) {
      return false;
    }

    order = (unsigned int) (order == 4) ? 5 : order;

    stereo = bondptr->stereo;

    if (stereo) {
      if (stereo == 1) {
        flags |= OB_WEDGE_BOND;
      }
      if (stereo == 6) {
        flags |= OB_HASH_BOND;
      }
    }

    if (bondptr->aromatic != 0) {
      flags |= OB_AROMATIC_BOND;
    }

    if (!mol.AddBond(start,end,order,flags)) {
      return false;
    }

    ++bondptr;
  }

  intptr = (unsigned int*) bondptr;

  mol.SetAromaticPerceived();
  mol.SetKekulePerceived();

  return true;
}
Then, when matching, I check if the query might contain tetrahedral stereo. If yes, I build the target OBMol fresh from SMILES, if no directly from the serialized object:
 if (strchr (querysmi, '@') != NULL)
    {
        //Match against an OBMol generated from SMILES
    }
    else
    {
         //Match against an OBMol deserialized from binary
    }

Tuesday, December 6, 2011

First Light

As promised, here the first results for substructure searching c1ccccc1Cl on the first 100000 compounds from Pubchem: select * from pubchem.compound where compound >= 'c1ccccc1Cl'::molecule

Substructure search speed


RankBuildStorageFingerprintHitsno IndexHitswith Index
1OpenBabelbinary+SMILESFP280709236 ms8067936 ms
2Indigobinaryext+sub804924821 ms80492418 ms
3OpenBabelSMILESFP2807057971 ms80675432 ms

I've checked why the OpenBabel FP2 fingerprint eliminates three structures that otherwise would pass: Using the VF2 OBIsomorphismMapper instead of the OBSmartsPattern, it's also 8067 without index. But it's about four times slower, 90526 ms without index, 8798 with index.
The structures in question are: 80944, 83450 and 99925 and I'm pretty sure it's caused by differences in aromaticity detection.

Stereochemistry in substructure searches


CheckQueryExpectedIndigo binaryOpenBabel SMILESOpenBabel binary
R/S differentselect 'C[C@H]([C@@H](C(=O)O)N)O'::molecule <= 'O=C(O)[C@H](N)[C@@H](O)C'::moleculefalsepassfailfail
R/S sameselect 'C[C@H]([C@@H](C(=O)O)N)O'::molecule <= 'C[C@H]([C@@H](C(=O)O)N)O'::moleculetruepasspasspass
E/Z different select 'C(=C\Cl)\Cl'::molecule <= 'C(=C/Cl)\Cl'::moleculefalsepassfailfail

OpenBabel fails the 'E/Z different' and 'R/S different' checks, but these are known issues up to version 2.3.1. More disturbing are the issues with 'R/S different'. For SMILES I can reproduce it with obgrep, so it's not my code causing it.

Indigo has no obvious issues with matching R/S and E/Z stereochemistry. Speedwise, OpenBabel with binary storage would be the leader of the pack, but has inconsistent behaviour compared with it's SMILES storage sibling. So the serialization/deserialization code clearly is missing something. But I've found a simple workaround: If the query contains chirality information it uses SMILES, otherwise binary. Speedwise, OpenBabel with binary storage now is the leader of the pack.

Fingerprint efficiency


Fingerprint typeCandidates screenedHits matchedfalse positivesEfficiency
ext+sub80908049410.995
FP281458067780.99

Pretty close.

Wednesday, November 30, 2011

What's up, doc?


OK, I have now a working (i.e. not crashing & giving meaningful results) hybrid version of pgchem::tigress with OpenBabel and Indigo, where the OpenBabel functions in the search engine have been replaced by Indigo functionality.

Why?

  • To see if it can be done :-)
  • To learn more about programming with Indigo
  • To see which is faster
  • Indigo has better reaction support as it seems

Interestingly I could keep the external interface almost stable. As only a dozen or so functions outside the search engine are affected by the change, I think with some #IFDEFs I can make a common codebase that compiles either to 100% OpenBabel or the hybrid version.

Monday, November 14, 2011

R-Matrix implementation with Indigo

Long time no see.

I have been trying programming with Indigo recently and here is my first, straightforward (to avoid the term 'naive') implementation of an Ugi like R-Matrix (Gasteiger/Engel, Chemoinformatics, 2003, Wiley-VCH, p. 185ff.) for adding make/break/change bond detection to its reaction mapper:
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include "indigo.h"

static inline int mp(const int column, const int row, const int maxRows)
{
    return column-1+((row-1)*maxRows);
}

int main(int argc, char *argv[])
{

    int index, reaction, reactant, product, atom, bond, atomIter, reactantIter, productIter, neighborIter, totalAtoms, mappingNumber, targetMappingNumber, bondMakeBreak;
    int *bMatrix, *eMatrix;//, *rMatrix;

    totalAtoms=0;

    qword session = indigoAllocSessionId();

    indigoSetSessionId(session);

    printf("%s\n",indigoVersion());

    reaction = indigoLoadReactionFromFile("c:/temp/testrxn.rxn");

    //indigoFoldHydrogens(handle);

    indigoDearomatize(reaction);

    //printf("1:%s\n",indigoGetLastError());

    indigoAutomap(reaction,"discard");

    reactantIter = indigoIterateReactants(reaction);

    while((reactant = indigoNext(reactantIter)))
    {
        totalAtoms+=indigoCountAtoms(reactant);
    }

    indigoFree(reactantIter);

    bMatrix = (int*) calloc(totalAtoms*totalAtoms,sizeof(int));
    eMatrix = (int*) calloc(totalAtoms*totalAtoms,sizeof(int));
    //rMatrix = (int*) calloc(totalAtoms*totalAtoms,sizeof(int));

    reactantIter = indigoIterateReactants(reaction);

    while((reactant = indigoNext(reactantIter)))
    {
        atomIter = indigoIterateAtoms(reactant);
        while((atom = indigoNext(atomIter)))
        {
            mappingNumber = indigoGetAtomMappingNumber(reaction,atom);
            //printf("%i\n",mappingNumber);
            //indigoGetExplicitValence(atom, &valence);
            //printf("%i\n",valence);
            //bMatrix[mp(mappingNumber, mappingNumber, totalAtoms)]=valence;

            neighborIter = indigoIterateNeighbors(atom);

            while((atom = indigoNext(neighborIter)))
            {
                targetMappingNumber = indigoGetAtomMappingNumber(reaction,atom);

                index = mp(mappingNumber, targetMappingNumber, totalAtoms);

                if(bMatrix[index]==0)
                {
                    bond=indigoBond(atom);
                    bMatrix[index]=indigoBondOrder(bond);
                }
            }

            //printf("\n");

            indigoFree(neighborIter);
        }
        indigoFree(atomIter);
    }

    indigoFree(reactantIter);

    productIter = indigoIterateProducts(reaction);

    while((product = indigoNext(productIter)))
    {
        atomIter = indigoIterateAtoms(product);
        while((atom = indigoNext(atomIter)))
        {
            mappingNumber = indigoGetAtomMappingNumber(reaction,atom);
            //printf("%i\n",mappingNumber);
            //indigoGetExplicitValence(atom, &valence);
            //printf("%i\n",valence);
            //eMatrix[mp(mappingNumber, mappingNumber, totalAtoms)]=valence;

            neighborIter = indigoIterateNeighbors(atom);

            while((atom = indigoNext(neighborIter)))
            {
                targetMappingNumber = indigoGetAtomMappingNumber(reaction,atom);

                index = mp(mappingNumber, targetMappingNumber, totalAtoms);

                //printf("%i",indigoBondOrder(bond));
                if(eMatrix[index]==0)
                {
                    bond=indigoBond(atom);
                    eMatrix[index]=indigoBondOrder(bond);
                }
            }

            //printf("\n");

            indigoFree(neighborIter);
        }
        indigoFree(atomIter);
    }

    indigoFree(productIter);

    for(int i=0; i<(totalAtoms*totalAtoms); i++)
    {
        bondMakeBreak=eMatrix[i]-bMatrix[i];

        //rMatrix[i]=bondMakeBreak;

        if(bondMakeBreak>0)
        {
            if(bMatrix[i]>0)
                printf("change: %i to %i\n",(i / totalAtoms)+1,(i % totalAtoms)+1);
            else
                printf("make: %i to %i\n",(i / totalAtoms)+1,(i % totalAtoms)+1);
        }
        else if(bondMakeBreak<0)
        {
            if(bMatrix[i]>1)
                printf("change: %i to %i\n",(i / totalAtoms)+1,(i % totalAtoms)+1);
            else
                printf("break: %i to %i\n",(i / totalAtoms)+1,(i % totalAtoms)+1);
        }
    }

    //printf("2:%s\n",indigoGetLastError());

    indigoSaveRxnfileToFile(reaction,"c:/temp/testrxn_mapped.rxn");

    //printf("3:%s\n",indigoGetLastError());

    indigoFree(reaction);
    indigoReleaseSessionId(session);
    free(bMatrix);
    free(eMatrix);
    //free(rMatrix);

    return 0;
}
At the moment it has the following shortcomings:
  1. 50% too much memory consumption, since only half of the R-Matrix would suffice
  2. Each bond is found twice (Atom1 to Atom2 and Atom2 to Atom1)
  3. The make/break/change info is not added to the mapped reaction file

Have fun!

Monday, May 23, 2011

Putting things into perspective

How to do it the wrong way:

Radiation in Japan

I am no physicist, but this is so straightforward ignorant of factors like long term exposure effects, ingestion of radionuclides and non-uniform contamination that it leaves me speechless.

Dear silicon valley executives with geiger counters, leave them at home. Please.

Monday, May 9, 2011

JSDraw available again

Apprently, the HTML/JavaScript chemical editor JSDraw is available again:

JSDraw

As being haunted by 'issues' with plugin type editors based on Java or ActiveX for years, it's nice to see that browser-native editing/depiction of chemical structures is gaining more and more momentum.

The other players I know of are:

ChemWriter
Ketcher
Sketcher 2D

Friday, April 29, 2011

Java in-memory compression

Sometimes it would be nice to have a means to compress seldom used objects in memory to time for space (since compression/decompression takes a bit of overhead time). Especially in chemoinformatics applications where objects often carry substantial compactible information, e.g. textual representations of structures, the compression ratio is quite noticable.

Sure, one could do it for every variable individually, but how about compressing whole Objects using the standard serializing mechanism?

If you use WEKA, you have the SerializedObject.

If not - enter the CompressedReference:

import java.io.ByteArrayInputStream;
import java.io.ByteArrayOutputStream;
import java.io.IOException;
import java.io.ObjectInputStream;
import java.io.ObjectOutputStream;
import java.io.Serializable;
import java.util.zip.GZIPInputStream;
import java.util.zip.GZIPOutputStream;

public class CompressedReference<T extends Serializable> implements Serializable {

private static final long serialVersionUID = 7967994340450625830L;

private byte[] theCompressedReferent = null;

public CompressedReference(T referent) {
try {
compress(referent);
} catch (IOException e) {
// TODO Auto-generated catch block
e.printStackTrace();
}
}

public int size() {
return theCompressedReferent.length;
}

public T get() {
try {
return decompress();
} catch (IOException e) {
// TODO Auto-generated catch block
e.printStackTrace();
} catch (ClassNotFoundException e) {
// TODO Auto-generated catch block
e.printStackTrace();
}
return null;
}

private void compress(T referent) throws IOException {

ByteArrayOutputStream bos = new ByteArrayOutputStream();
GZIPOutputStream zos = new GZIPOutputStream(bos);
ObjectOutputStream ous = new ObjectOutputStream(zos);

ous.writeObject(referent);

zos.finish();

bos.flush();

theCompressedReferent = bos.toByteArray();

bos.close();
}

@SuppressWarnings("unchecked")
private T decompress() throws IOException, ClassNotFoundException {
T tmpObject = null;
ByteArrayInputStream bis = new ByteArrayInputStream(theCompressedReferent);
GZIPInputStream zis = new GZIPInputStream(bis);
ObjectInputStream ois = new ObjectInputStream(zis);
tmpObject = (T) ois.readObject();

ois.close();

return tmpObject;
}
}

A quick test shows 528 byte size for a String of 250 characters (since Unicode needs two bytes per char) and 64 bytes after compression, a ratio of about 8:1. The only requirement is that the Object stored in the reference has to implement Serializable.

And yes, I have to clear those TODO reminders... ;-)

Update
The compression ratio for a CDK Molecule representation of 2-Fluoronaphthalene is about 3:1.
The compression ratio for the V2000 mofile of 2-Fluoronaphthalene is about 7:1.

Friday, November 26, 2010

pgchem::tigress sets new world record!

The IPB now has successfully loaded ~28 million structures from pubchem into a pgchem::tigress / PostgreSQL database.

This is the largest installation I know of and therefore sets a new world record - for pgchem. :-)

Friday, November 12, 2010

Development of pgchem::tigress 2.0 has started

I finally managed to get OpenBabel 2.3.0 and PostgreSQL 9.x to talk to each other on Windows. This only works when building with the MinGW Makefiles, not with MSYS, as PostgreSQL does not like mixing with the MSYS runtime libraries. And still a static build of libopenbabel and libinchi is required.

The overall good news is that the Makefile could be simplified a lot with help from Steffen Neumann from the Leibniz Institute of Plant Biochemistry and all current functions are still working with OB 2.3. The build process on Windows is still a bitch compared to building on Posix systems, though.

The roadmap for 2.0:

  • A lot of code needs cleanup, especially the make_molecule() function, which is unreadable and error prone. Mr. Neumann has also identified various possible memory leaks that need to be fixed
  • Allow 2D coordinate generation for Molfile output
  • Using OBQuery for non-SMARTS substructure searching
  • Add Spectrophore(TM) output
  • Add support for exporting Andrew Dalke's chemfp formats  (at least FPS)
  • Move SVG output from Dingo to OB, having one library less to care for
  • Add the fingerprint  linear optimizer in a generally usable implementation
  • See if the new stereochemistry inside OB fixes various issues that were reported
  • Basic reaction support
  • Win64 builds ?

This will take some time, so don't hold your breath. Since there are already some bugs reported for OB 2.3.0 I guess I'll have to wait for those being fixed before a production release of pgchem::tigress 2.0 becomes available.

Tuesday, November 9, 2010

GCC6 Day Three (09.11)

So, I rediscovered my notes from day three. Unfortunately there are only two lectures I can really write something about:

  1. Noel O'Boyles lecture about the in-silico design of polymers with optimal properties for organic solar cells. As far as I understood this, the predicted efficency for organic solar cells is about 13 %. Experimentally 6 % have been reached yet. By using OpenBabel together with cclib and Gaussian, synthetically accessible polymers with 11 % predicted efficiency could be designed.
  2. Roger Sayle from NextMove showed his recent work on a chemistry aware spell checker, based on the observation, that the current problems with chemical text mining do not come from poor OCR or poor name2structure conversion, but mainly from bad input because of typos etc. in the source texts and common spell checking software cannot cope with this.

Monday, November 8, 2010

GCC6 Day Two (08.11)

So, I rediscovered my notes from day two. Unfortunately there are only two lectures I can really write something about:
  1. The WizePairZ  algorithm was quite interesting. It aims to solve the problem that often unwanted properties of drug candidates scale in correlation with wanted properties. The idea behind WizePairZ is to automatically find transformations from one moelcule to another that reduce unwanted properties (ideally while improving wanted properties). Then such transformations are applied to molecules which are close to the boundary between 'unwanted' and 'wanted' space in order to push them over the boundary.
  2. The in-silico prediction of Phototocicity aims to find molecular descriptors and associated models that allow to predict the potential phototoxicity of substances, especially their UV absorption between 290 - 450 nm wavelength where the human skin doesn't filter anymore and the ozone layer doesn't filter yet.


GCC6 Day One (07.11)


Arrived 14:00 in Goslar

What's new in Knime presentation:
Mainly usability improvements, like conditional paths through workflows and annotations. This is really getting somewhere. Release planned for December 6th 2010

My presentation:
I was not boohed from the stage

MOSGrid presentation:
It's a molecular simulation grid driven by a consortium of academia and industry, lead by the University of Cologne

The planned usage is, that you can specify a task through a web interface in MSML (molecular simulation markup laguage), which describes a task for the grid, then the MSML is translated into input specific to a program (e.g. Gaussian).
The program is then run on the grid and the program specific output is translated again into MSML and returned to the user.

Currently

  • Gaussian
  • Turbomole

are supported.

Tasks can be chained into workflows. Lots of expertise and manpower, and if buerocracy doesn't step in the way, this could work and is worth a look. Unfortunately, the licensing of the used programs in a grid seems to be an unsolved legal problem.

Demonstration of MOSGrid:
A basic web interface (to Gaussian I guess) was shown and how to run quantum calculations on the grid from there. The interface is designed for humans. Currently there are no plans to offer this as a service to machines. Licensinsg issues of using commercial software in a grid must be ironed out to make it widely available, but beta testing will start in December to a limited number of users.

Friday, November 5, 2010

Thursday, November 4, 2010

The SMARTSViewer

This is a nifty service for visualizing SMARTS patterns from the University of Hamburg:


Probably useful when your SMARTS pattern does not match what you think it should...

Wednesday, October 27, 2010

Nasty thoughts

Typically reading the FDA Warning Letters is pretty dry stuff, but this looks like a fun product.

"The effects of ... will promote a thrilling energy, more stamina, constant readiness, nasty thoughts, prolonged arousal, feelings of well being, romantic and sensual experience."

"Nasty thoughts" - those guys from Natural Wellness, LLC really should try Mefloquine...

Wednesday, August 11, 2010

Searching Structures in ChemSpider from the Android Browser...

...is now implemented as empirical formula search.

Not what was ultimately desired, but at least it works.

Friday, August 6, 2010

Searching Structures in ChemSpider from the Android Browser...

...does not work at the moment.

Having, with the friendly help of ChemSpider support, found out how to use their web API, I tried to call it from the android browser.

The HTML form works at least in Firefox, IE 6.x and Chrome, but from the Android browser it either opens the ChemSpider start page - or sometimes only a blank page.

It seems that I'll have to monitor what the Android browser really sends to ChemSpider with Wireshark to find out what the problem is.


Friday, July 23, 2010

Thesis abstract

The storage and retrieval of chemical graphical datatypes such as structures and reactions in relational database systems is a common technique used in academia and industry alike. Due to the computationally intensive algorithms used for (sub)graphisomorphism detection, such systems commonly use faster screening mechanisms in order to reduce the set of potentional match positives before applying aforementioned algorithms.

Widely used screening mechanisms are based on numerical and binary vectors, called fingerprints, with a clear dominance of binary fingerprints due to the raw speed advantage of bitwise operations and compactness in storage. The two most commonly used types of binary fingerprints are path-generated and substructure-generated, both of which have specific shortcomings, especially blind spots.

To overcome this shortcomings, the Pgchem::Tigress chemistry extension to the PostgreSQL object-relational database management system uses a hybrid binary fingerprint, consisting of an invariant path-generated part and an substructure-generated part which is externally configurable through a dictionary of substructure patterns.

This thesis presents a novel approach of using dynamic discrete optimization to find an optimized dictionary configuration for the substructure-generated part of the fingerprint for arbitrary sets of structural data.

By means of applying the method developed in this thesis, the computational power neccessary to run a chemical information system can be reduced by 42 percent on average. By improving the query throughput, upgrading the server hardware to the next level of computational power can be avoided and thus opportunity revenues of the operating costs are realized.

Update: It is now readable online.

Monday, July 19, 2010

Thesis defence

The date of my defence is 21.07.2010, 10:00-11:15 UTC+2.

If somebody's willing to cross their fingers for me, it'll be appreciated. :-)

Friday, July 9, 2010

ChemSpider Web API, anyone?

I'm trying to do a structure search using the published ChemSpider Web API:  http://www.chemspider.com/WebAPI.aspx

I thought it would be sufficient to get the form they provide, put a molfile into the designated area, add an submit button and it works.

But all I get back are blank pages with no content. No error messages, nothing useful.

How is the darn thing meant to be used?

UPDATE: Adding a method="post" to the form did the trick, but apparently only substructure searches with unlimited result set are supported. This results in very poor performance with structures of low selectivity. Still not good

UPDATE 2:

  1. I had an extra line in my molfile.
  2. There is really a bug in the API. Substructure search triggers exact search and vice versa. Since the API is heavily in use with workarounds to this bug, fixing this would break many applications.


Monday, July 5, 2010

The curious case of the infinite canvas

After some fruitless attempts to extend the usable canvas of DCE ChemPad with some kind of self-made view port, in order to allow drawing of larger-than-screen structures, I think I now found a more painless way.

Simply wrapping the DCE Control in a HorizontalScrollView was not enough, because the ScrollView intercepts too much MotionEvents, but after subclassing HorizontalScrollView and overriding onInterceptTouchEvent(), it seems to work as expected. The code needs some more polishing though to call it working though.

































As you can see on the screenshot, the structure now can be scrolled left and right as needed to draw more than would fit on the screen. The thin line on the bottom separates the drawing from the scrolling sensitive area, because otherwise the whole screen would be scroll sensitive and drawing while scrolling is not very precise. The screenshot is from a NexusOne with Android 2.2.

Wednesday, June 2, 2010

DCE ChemPad Update 1.1

A few bugs have been fixed (sorry, no autolayout yet :-)), all controls have been moved into the menu and no titlebar anymore to save screen estate. The vibration feedback for fusing atoms works now as expected and you can send the molfile via the phone's messaging systems. Since it seems that attachment handling is broken not fully implemented on Android, the molfile is embedded as text in the message itself. A localized help function was added to the application.

Wednesday, May 26, 2010

DCE ChemPad 1.0 is out in the wild

DCE ChemPad is released for free on the Android Market!

It shows the capabilities of the DCE Chemistry Editor Control to add a chemical editor to arbitrary Android applications.

It was tested on the HTC Nexus One, ACER Liquid S 100, Motorola Milestone/Droid and the Emulator and should work with all Android versions >= 1.5. It does reportedly not work on the Motorola Cliq.

If you have acceess to the Market, please try it and tell me what you think...

Ah, screenshots:

Thursday, April 22, 2010

Solving the 'big finger vs. tiny screen' problem

  • Learn how to make custom View controls - check
  • Learn to make compound controls - check
  • Learn how to paint on the canvas - check
  • Design a effective 2D rendering pipeline for undirected graphs- check
  • Remember basic planar trigonometry - check
  • Design a fuzzy lock-on selection method for the touchscreen - check
  • Design a robust backing model - check
  • Design an effective UI for the touchscreen - check
  • Implement all the nasty details - check

Wednesday, April 14, 2010

NexusOne: Big finger vs. tiny screen

I now have a NexusOne at hand.

And while it is much faster than the Emulator, a very profane problem has come up. It is impossible to precisely draw a chemical structure with an editor designed for the mouse! It is just to sensitive to use it with a finger on a tiny screen. And a pen won't work with a capacitive touchscreen...

Maybe i have to take an intensive look at the android.graphics package.

Thursday, April 1, 2010

JavaScript molecule editor roundup

Next to JsDraw which had it's 1.0 release recently, I've found two more pure JavaScript molecule editors:

WebCME which has a lot of features, notably a large library of molecules but is quite painful to use. This is because it's developers have chosen a system of 'select two atoms, add bond, deselect them, select another two atoms, add bond, oh wrong one, delete bond, add correct bond...' for drawing.

The ChemDoodle web components on the other hand have a totally minimalistic, yet powerful UI. Atoms are drawn by mouseclick, bonds drawn by mouse drag. Atom types can be changed via keyboard, bond types by mouseclick and delete is done by the Backspace or Delete key.

Unfortunately none of the three works on mobile browsers. Either they don't work at all or only by half.

In contrast, the jsMolEditor does work even in Android WebViews, but seems to be not under development anymore. I suspect that it was abandoned in favour of JsDraw.

What a pity. Having a (even simple) but working JavaScript molecule editor, that works in Smartphone WebViews would open a whole new world of applications for those devices. They now have the computational power to handle chemical data, but who wants to enter SMILES strings by hand...

pgchem::tigress 1.2 is out

http://pgfoundry.org/projects/pgchem/

Built and tested against PostgreSQL 8.4.2 with OpenBabel 2.2.3 on XP 32 bit, Windows 7 64 bit and Ubuntu 8.04 LTS 32 bit.

MACCS166 binary fingerprints.

Dice and Tversky similarity.

Small bug fixes.

Tuesday, March 16, 2010

How to run mx on Android

mx runs on Android!

To make it work, you need to build from sources and remove all references to javax.swing (which doesn't seem to break the rest of the code btw.), since Android does not contain AWT or Swing.

Then repackage the jar and it can be used in Android applications.

Chemistry on the smartphone, yay. :-)

Wednesday, March 10, 2010

Changing the public API between minor versions of PostgreSQL

I hate when they do this.

Between 8.3 and 8.4, the API for CREATE OPERATOR CLASS has changed. Now the RECHECK flag is obsolete, letting the index dynamically decide if it is lossy or not.

While this in itself is an improvement, it generates an incompatibility between GiST C code and scripts written for 8.3 and 8.4. Fortunately, the fix seems to be an easy one...

Wednesday, March 3, 2010

Wrapping native libraries with JNA: Dingo

Even with the advent of pure Java chemoinformatics toolkits like MX or the CDK, there is a lot of interesting native code floating around on the net. Unfortunately, wrapping native code with JNI is no real fun.

JNA comes to the rescue. It does all the neccessary loading and marshalling stuff dynamically in the background for you. All you need is a declaration of the native interface, the rest is magic.

Here's an incomplete but working example for Dingo 1.0:

package your_package_here;
import com.sun.jna.Native;
public class NativeDingoWrapper {

  static {
    Native.register("dingo");
  }

  public static native int dingoSetOutputFormat(String anOutputFormat);
  public static native int dingoSetColoring(int aColoringFlag);
  public static native int dingoSetHighlightColorEnabled(int aHighlightFlag);
  public static native int dingoSetHighlightThicknessEnabled(int aHighlightThicknessFlag);
  public static native int dingoSetStereoOldStyle(int aStereoFlag);
  public static native int dingoSetImageSize(int aWidth, int aHeight);
  public static native int dingoLoadMolFromString(String aMol);
  public static native int dingoLoadMolFromFile(String aFile);
  public static native int dingoSetOutputFile(String anOutputFile);
  public static native int dingoRender();
  public static native int dingoMoleculeIsEmpty();
}


And that's it.

The only drawback of JNA is that it needs a glue DLL specific to the operating system, so theoretically it is more platform limited than JNI.

But since "JNA has been built and tested on OSX (ppc, x86, x86_64), linux (x86, amd64), FreeBSD/OpenBSD (x86, amd64), Solaris (x86, amd64, sparc, sparcv9) and Windows (x86, amd64). It has also been built for windows/mobile and Linux/ppc64, although those platforms are not included in the distribution." this is a quite limited limitation for most cases.

I have successfully wrapped Dingo and Barsoi with JNA and up to now it just works as advertised.