Monday, October 27, 2014

The case for map/reduce with relational databases

I was at pgconf.eu and held a lightning talk about map/reduce with PostgreSQL. Upfront, I was asked "Why do you want to do that anyway?" and my initial response was like, "Because I can.". :-)

But that got me thinking about the real case behind the idea. What is the heart of map/reduce? Citing from the original paper:

"MapReduce is a programming model and an associated implementation for processing and generating large data sets." Note the word "sets"?

"Users specify a map function that processes a key/value pair to generate a set of intermediate key/value pairs, and a reduce function that merges all intermediate values associated with the same intermediate key.", or, expressed in a more formal way

map (k1,v1)                -> list(k2,v2)
reduce (k2,list(v2))    -> list(v2)


Please note, that this is a bit more precise than the initial definition. The paper explains:

"I.e., the input keys and values are drawn from a different domain than the output keys and values. Furthermore, the intermediate keys and values are from the same domain as the output keys and values."

This is important, because the authors are introducing a domain transformation of the input data here. That is, in my opinion, already the heart of map/reduce.

Going back to the initial definition, this is basically what all RDBMS already do when processing parallel queries, be it by builtin ability or bolted on like with PL/Proxy + PostgreSQL: In the first step the input set is broken down to partitions, then the query runs on that partitions in parallel and produces intermediate result sets and finally that intermediate result sets are aggregated to the final result set. But the formal definition above adds a little twist, the domain transformation.

To clarify this, I'll use the canonical example, counting words in a text. The map function converts semi structured data, a text with lines of arbitrary length, into a well structured set of key (a word) and value (its count) tuples. This is the difference and the key to the power of map/redcue.
The ability to handle semi structured data which the relational model usually does not handle very well. (And I won't say unstructured data. Truly unstructured data is statistical noise.)

But modern RDBMS, especially PostgreSQL, often already have functions to transform semi structured data into relations and/or allow for user defined functions to extend their capabilities and that allows for running a map/reduce type job inside a RDBMS. Still, why would somebody want to do this?

1.) Integration
An awful lot of data is stored in relational models and will stay there. Simultaneously, especially for analytical workloads which become more and more important, the need for integrating relational and semi-structured data grows. Why handle them in different systems when one will do?
This decision of course heavily depends on the real world requirements. But rest assured that the datacenter guys who have to run the show will like to operate one database better than 2..n.

2.) Sets
Remember the word "sets" from the initial definition? And now the definition of a "relation" in a RDBMS:

"R is a relation on these n domains if it is a set of elements of the form (d1, d2, ..., dn) where dj ∈ Dj for each j=1,2,...,n." (E. F. Codd (Oct 1972). "Further normalization of the database relational model". "Data Base Systems". Courant Institute: Prentice-Hall. ISBN 013196741X.)

If a relation is a set of tuples with values from some domain D and map/reduce does domain transformation on key/value pairs (aka. tuples) what does that call for? Right, a very efficient set processor. Since relational DBMS are very efficient set processors by nature, they allow for
writing compact map/reduce functions that are also less error prone due to the declarative nature of SQL.

To clarify what I mean take a look at the following map and reduce functions for wordcount written for MongoDB in JavaScript from here:

var map = function() { 
    var summary = this.summary;
    if (summary) {
        // quick lowercase to normalize per your requirements
        summary = summary.toLowerCase().split(" ");
        for (var i = summary.length - 1; i >= 0; i--) {
            // might want to remove punctuation, etc. here
            if (summary[i])  {      // make sure there's something
               emit(summary[i], 1); // store a 1 for each word
            }
        }
    }
};

var reduce = function( key, values ) {   
    var count = 0;   
    values.forEach(function(v) {           
        count +=v;   
    });
    return count;
}

db.so.mapReduce(map, reduce, {out: "word_count"})


PostgreSQL:

For the code, see my previous post.

Well, while it seems to require more code than MongoDB, there is a subtle difference. The most PostgreSQL code is standard boilerplate to write a set returning function and to make PL/Proxy work. Once you got that right, you usually never have to look back. The actual work is done in two lines of SQL:

SELECT TRIM(both from word),count(1) FROM (SELECT regexp_split_to_table(line, E'\\W+') as word FROM kjb) w GROUP BY word

and 

SELECT word,sum(count) FROM map_kjb() AS (word text, count bigint) WHERE word != '' GROUP BY word

IMHO, the ability to express this with extensively tried and tested functions instead of having to implement them yourself combined with strong typing is worth so much, that one should give it a try before looking somewhere else. Granted, raw performance may become so paramount that an alternative technology might be called for, but if you already use PostgreSQL now there is another reason to like and not to leave it.

And this is it, the long answer I didn't already had ready at the conference.

P.S. It also allows to move computation instead of data around.

Monday, October 20, 2014

pgconf.eu 2014

On the worker nodes:

 CREATE TABLE kjb
(
  line text
);



CREATE OR REPLACE FUNCTION map_kjb()
  RETURNS SETOF record AS
$BODY$
DECLARE r record;
BEGIN
FOR r IN (SELECT TRIM(both from word),count(1) FROM (SELECT regexp_split_to_table(line, E'\\W+') as word FROM kjb) w GROUP BY word) LOOP
RETURN NEXT r;
END LOOP;
RETURN;
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 1000
  ROWS 1000;


On the head node:

CREATE TYPE t_mr_wc AS (word text, count bigint);

CREATE OR REPLACE FUNCTION map_kjb()
  RETURNS SETOF record AS
$BODY$
CLUSTER 'head';
RUN ON ALL;
$BODY$
  LANGUAGE plproxy VOLATILE
  COST 1000
  ROWS 1000;


CREATE OR REPLACE FUNCTION reduce_kjb()
  RETURNS SETOF t_mr_wc AS
$BODY$ DECLARE result public.t_mr_wc; BEGIN
FOR result IN select word,sum(count) from map_kjb() AS (word text,count bigint) where word != '' group by word LOOP RETURN NEXT result; END LOOP; END; $BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100
  ROWS 100;


1. Load data into tables.

2. Get word count:

SELECT * FROM reduce_kjb();

3. Get top 21 word counts: 

SELECT * FROM reduce_kjb() ORDER BY 2 DESC LIMIT 21; 

Wednesday, September 3, 2014

Using Code::Blocks with the Intel C/C++ compiler on Windows in five minutes

Code::Blocks apparently does not support the latest Intel C/C++ compiler out of the box. Trying to set up the necessary include paths for both the Visual Studio compiler and the Intel compiler leads to madness.

But there is an easier way:

Make two Intel compiler entries in Code::Blocks for the 32 and 64 bit version with the correct base directory, e.g. C:\Program Files (x86)\Intel\Composer XE\bin.
The 32 bit toolchain executables are prefixed ia32\, the 64 bit ones intel64\.

When you have installed the Intel compiler, there is an 'Command Prompt with Intel Compiler' entry in the start menu for 32 and 64 bit targets each. This opens a command prompt with all necessary paths and environment variables set correctly.

Just start codeblocks.exe from this command prompt and it inherits the environment. Then you can select the 32 or 64 bit Intel compiler option in the 'Project->Build options' and it works without further ado.

The only thing you cannot do is switching between 32 and 64 bit targets on the fly, you have to start Code::Blocks with the correct environment first.

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.