The Long Tail - vertical table partitioning I

DISCLAIMER: This is just an idea, I don't have tried this in a production environment!

Having said that, any input is welcome. :-)

Storing tons of rows in a table of which only a small percentage of rows are frequently queried is a common scenario for RDBMS, especially with databases that have to keep historical information just in case they might be audited, e.g. in environments regulated by law.

With the advent of the first affordable 8TB harddisk and fast SSDs still being much more expensive per GB, I wondered if such long-tailed tables could be split over a SSD holding the frequently accessed pages and a near-line storage HDD keeping the archive - elegantly - with PostgreSQL.

With elegant, I mean without fiddling around with VIEWs, INSTEAD OF triggers and exposing a clean and familiar interface to the developer.

OK, since PostgreSQL already supports horizontal partitioning, spreading one table transparently over many parallel tables, how about vertical partitioning, spreading one table over a hierarchy of speed?

The speed zones can be mapped to tablespaces:

CREATE TABLESPACE fast LOCATION '/mnt/fastdisk';
CREATE TABLESPACE slow LOCATION '/mnt/slowdisk';

Next comes the table(s):

CREATE TABLE the_table
(
  id integer NOT NULL,
  value real
)
WITH (
  OIDS=FALSE
)
TABLESPACE fast;

CREATE TABLE the_table_archive
(
)
INHERITS (the_table)
WITH (
  OIDS=FALSE
)
TABLESPACE slow;

Table inheritance in PostgreSQL is so cool...

And a function to move data from fast to slow:

CREATE OR REPLACE FUNCTION move_longtail()
  RETURNS boolean AS
$BODY$
DECLARE worked BOOLEAN;
DECLARE rowcount INTEGER;
BEGIN
worked := false;
rowcount := count(*) FROM ONLY the_table WHERE id >= 5000000;
IF (rowcount > 100000) THEN
INSERT INTO the_table_archive SELECT * FROM ONLY the_table WHERE id >= 5000000;
DELETE FROM ONLY the_table WHERE id >= 5000000;
worked := true;
END IF;
RETURN worked;
END;
$BODY$
  LANGUAGE plpgsql VOLATILE STRICT;

This function runs only if a minimum of movable rows qualify. This is recommended since SMR disks like large contiguous writes due to how SMR technically works.

Notice the ONLY keyword. This allows to control precisely which partition of the table is affected by the DML statement. Did I say already that table inheritance in PostgreSQL is so cool?

And basically that's it. All the developer sees is SELECT * FROM the _table; or something, being oblivious to the underlying machinery.

Some kind of automatic maintenance is missing. On to part II...

pgchem::tigress 3.2 released

pgchem::tigress 3.2 is finally out!

• This builds against PostgreSQL 9.4.x and OpenBabel 2.3.2 on Linux.
• It contains all fixes and contributions of the previous versions.
• Windows is not supported anymore - and since it builds and runs way better on Linux, probably never will be again.
• Depiction functions have been removed. Their run time was too unpredictable to be run inside a database server.
• Theoretical isotope pattern generation with MERCURY7 is now available with isotopes for 39 elements.

So: CREATE EXTENSION pgchem_tigress;

Counting rows again - inheritance strangeness solved

Taking a second look at the execution plans, I've noticed that the scan on slow read twice the number of pages from disk than the one on fast:

Buffers: shared read=442478 vs. Buffers: shared read=221239

Since I loaded all rows into slow first and then moved 50% of them into fast, this makes sense, I guess.
If I understand it correctly, those pages in slow are now empty, but PostgreSQL keeps them for future use.

So I tried a VACUUM FULL on slow and ran my queries again. That changed the plans:

Buffers: shared read=221239 vs. Buffers: shared read=221239

And execution times are now about equal. The observed effect had nothing to do with table inheritance at all.

So VACUUM FULL can help, when moving big chunks of data between tables with otherwise low write activity.

Counting rows again - inheritance strangeness

Hm, for something I'm trying at the moment on 9.4, I created two identical tables where the second one inherits all from the first:

CREATE UNLOGGED TABLE everything.slow
(
  id integer NOT NULL,
  value real
)
WITH (
  OIDS=FALSE
)

CREATE UNLOGGED TABLE everything.fast
(
)
INHERITS (everything.slow)
WITH (
  OIDS=FALSE
)

No indexes. Default tablespace.

If I then load 50 million records into each of those tables and query them individually using the ONLY restrictor, a count(*) on the parent table (slow) is slower than on the descendant (fast):

select count(*) from only everything.slow;

"Aggregate  (cost=1067783.10..1067783.11 rows=1 width=0) (actual time=4973.812..4973.813 rows=1 loops=1)"
"  Output: count(*)"
"  Buffers: shared read=442478"
"  ->  Seq Scan on everything.slow  (cost=0.00..942722.08 rows=50024408 width=0) (actual time=1012.708..3416.349 rows=50000000 loops=1)"
"        Output: id, value"
"        Buffers: shared read=442478"
"Planning time: 0.118 ms"
"Execution time: 4973.901 ms"

select count(*) from only everything.fast;

"Aggregate  (cost=846239.00..846239.01 rows=1 width=0) (actual time=3988.235..3988.235 rows=1 loops=1)"
"  Output: count(*)"
"  Buffers: shared read=221239"
"  ->  Seq Scan on everything.fast  (cost=0.00..721239.00 rows=50000000 width=0) (actual time=0.101..2403.813 rows=50000000 loops=1)"
"        Output: id, value"
"        Buffers: shared read=221239"
"Planning time: 0.086 ms"
"Execution time: 3988.302 ms"

This works with other aggregates like avg() too.

I had expected some overhead when querying without ONLY on slow, because of the traversal of the inheritance hierarchy, but not when I restrict the query to a specific table with ONLY...

Can someone explain this? UPDATE: Maybe I found it myself. See next post for my explanation.

Finding mass spectra with PostgreSQL: Indexing

When naively searching for a specific spectral contrast angle, e.g.

select id from spectrum, round(spectral_contrast('[{"m/z":144.1150218,"intensity":0.908209840784744},{"m/z":145.118459841064,"intensity":0.0841924148716925}]'::json, id),2) as spectral_contrast where spectral_contrast = 1.0;

performance is not really stellar on my ~ 270k spectra database:

Execution time: 18695.619 ms for 71 hits.

This is no surprise, since the database has to calculate the spectral contrast angle for all rows and filter each row. How about indexing?

The simplest method I can envision, is to precalculate the number of peaks for each spectrum, index the column and add it to the query. Since the spectral contrast angle requires the number of peaks of the compared spectra to be equal, this is a safe restriction.

select id from spectrum, round(spectral_contrast('[{"m/z":144.1150218,"intensity":0.908209840784744},{"m/z":145.118459841064,"intensity":0.0841924148716925}]'::json, id),2) as spectral_contrast where spectral_contrast = 1.0 and num_peaks = 2;

Execution time: 797.346 ms for 71 hits.

OK, but can we do better? Yes, by using a more selective restriction, like the lower and upper bounds of the m/z values of the spectra, index the column and add it to the query. These can easily be obtained by min(m/z)and max(m/z) and then be stored in a numrange column. Since the spectral contrast angle requires the m/z values of the peaks of the compared spectra to be equal, this again is a safe restriction.

select id from spectrum, round(spectral_contrast('[{"m/z":144.1150218,"intensity":0.908209840784744},{"m/z":145.118459841064,"intensity":0.0841924148716925}]'::json, id),2) as spectral_contrast where spectral_contrast = 1.0 and mzrange = numrange(144.1150218,145.118459841064,'[]');

Execution time: 35.128 ms for 71 hits.

About 534 times faster.

Using a range type instead of storing the lower and upper bound in discrete columns gives more flexibility when querying with other criteria, like the Tanimoto index. Then, PostgreSQL's range operators, like &&  (overlap), might come in handy.

Finding mass spectra with PostgreSQL: Spectra as function arguments

If the spectra are represented as tables/recordets, like in my last posts on the topic, the question arises how to use a table as an argument to a function in PostgreSQL without creating a custom datatype or using temporary tables.

1. As one dimensional array of m/z values followed by the corresponding peak values. The array can then be cut in half and unnested inside the function into a set of records:

select unnest(sdarr[1:1]) as "m/z", unnest(sdarr[2:2]) as intensity from ...;

2. As two dimensional array of m/z values and their corresponding peak values. The partial arrays can then be unnested inside the function into a set of records:

select unnest(tdarray[1:3]) as "m/z", unnest(tdarray[4:6]) as intensity from ...;

3. As JSON of m/z values and their corresponding peak values. The JSON can then be converted inside the function into a set of records:
 
select * from json_to_recordset( '[{"m/z":1.0,"intensity":2.2},{"m/z":3.3,"intensity":4.8}]') as x("m/z" double precision, intensity double precision);

All statements above generate a recordset like:

m/z    intensity
1.0    2.2
3.3    4.8
...    ...

In an application, I'd go with JSON, since it has the most understandable structure for a developer against this API and it does not require fiddling around with array support functions of the driver, e.g. like createArrayOf() in JDBC.

count(*) is faster than count(1) - no more with nine.four

In PostgreSQL 9.4, count(1) speed is equally fast to count(*).

Another reason to upgrade to 9.4.

count(<column>) is still slower if no index is used...

count(*) is faster than count(1) revisited

Whoa, I didn't expect that much resonance to the last post. :-)

"Can you post the `explain analyze` output of both queries? I'm pretty sure this is normal "distribution" caused by other stuff going on in the system"

OK. Since my 'test' was on a VirtualBox VM with a self-compiled PostgreSQL and this might not be a reliable test setup, so I hijacked one of our real PostgreSQL servers:

uname -a
Linux  2.6.32-358.23.2.el6.x86_64 #1 SMP Sat Sep 14 05:32:37 EDT 2013 x86_64 x86_64 x86_64 GNU/Linux

Intel(R) Xeon(R) CPU E5-2670 v2 @ 2.50GHz ,16 GB RAM, storage on SAN, ext4.

select version();
"PostgreSQL 9.3.5 on x86_64-unknown-linux-gnu, compiled by gcc (GCC) 4.4.7 20120313 (Red Hat 4.4.7-4), 64-bit"

create table hmr(id serial, value real);

insert into hmr (value) select random()*10000000 from generate_series(1,10000000);

CREATE INDEX idx_id
  ON hmr
  USING btree
  (id);

vacuum analyze hmr;

-----------------------------------------------------------------------

explain analyze select count(1) from hmr;

"Aggregate  (cost=169248.00..169248.01 rows=1 width=0) (actual time=1439.521..1439.521 rows=1 loops=1)"
"  ->  Seq Scan on hmr  (cost=0.00..144248.00 rows=10000000 width=0) (actual time=0.006..642.276 rows=10000000 loops=1)"
"Total runtime: 1439.548 ms"

select count(1) from hmr; --811 msec.

-----------------------------------------------------------------------

explain analyze select count(*) from hmr;

"Aggregate  (cost=169248.00..169248.01 rows=1 width=0) (actual time=1299.034..1299.034 rows=1 loops=1)"
"  ->  Seq Scan on hmr  (cost=0.00..144248.00 rows=10000000 width=0) (actual time=0.006..644.912 rows=10000000 loops=1)"
"Total runtime: 1299.061 ms"

select count(*) from hmr; --670 msec 

-----------------------------------------------------------------------

explain analyze select count(id) from hmr;

"Aggregate  (cost=169248.00..169248.01 rows=1 width=4) (actual time=1576.046..1576.046 rows=1 loops=1)"
"  ->  Seq Scan on hmr  (cost=0.00..144248.00 rows=10000000 width=4) (actual time=0.004..636.076 rows=10000000 loops=1)"
"Total runtime: 1576.069 ms"

select count(id) from hmr; --920 msec

----------------------------------------------------------------------- 

And the winner is: still count(*) with count(1) second and count(id) third.

This could be an explanation for count(1) being slower:

"So count(1) is explicitly passing the 1 in and checking it for being NULL before incrementing a counter while count(*) is simply incrementing the counter."

And count(id) might be slower, because it takes additional time to consider using the index?
Or because it's actually reading data: width=4 vs. width=0?

I'm out of my league here.

count(*) is faster than count(1)

Today I've discovered incidentally, that count(*) is slightly faster than count(1) on PostgreSQL 9.3.5:

create table  hmr(id serial, value real);

insert into hmr (value) select random()*10000000 from generate_series(1,10000000);

select count(1) from hmr; --731 msec

select count(*) from hmr; --661 msec

Used repeatedly, e.g. in a function, this can really add up and make a difference.
In the spectral contrast angle CTE from the previous post, it reduces runtime by two seconds or 24%.

I would have expected the constant to be faster than actually reading the row...

Finding mass spectra with PostgreSQL: Spectral contrast angle in SQL

The spectral contrast angle of two spectra is another, according to literature, one of the methods for comparing spectra by similarity.

The spectral contrast angle S is calculated by building a intensity vector in N-dimensional space for each spectrum, where N is the number of m/z peaks in the spectrum, and finding the cosine of the angle between the two vectors. If both are the same, S is 1.0 if the two are orthogonal, S is 0.0.

The formula is S = sum(a.intensity*b.intensity) / sqrt(sum(a.intensity^2)*sum(b.intensity^2)) with a and b being the spectra to compare.

On a table like

CREATE TABLE <spectra_table>
(
  id integer NOT NULL,
  "m/z" numeric NOT NULL,
  intensity numeric NOT NULL
)

this can easily be done in PostgtreSQL with a common table expression:

WITH SQ AS (select "m/z", intensity from <spectra_table> where id = 1 order by "m/z" ASC), ST AS (select  "m/z", intensity from <spectra_table> where id = 2 order by "m/z" ASC),
qc as (select count(1) as qp from sq), tc as (select count(1) as tp from st)
select sum(sq.intensity*st.intensity)/sqrt(sum(sq.intensity^2)*sum(st.intensity^2)) as spectral_contrast from SQ,ST,qc,tc where qc.qp = tc.tp and SQ."m/z" = ST."m/z"

One interesting property of S is, that it evaluates to 1.0 if the intensities of the two spectra are different, but only by a multiple of an integer. If I understood it right, this is correct, because this means that you have measured the same composition at different concentrations.

Put into a function:

CREATE OR REPLACE FUNCTION spectral_contrast(int, int) RETURNS numeric
    AS ' WITH SQ AS (select "m/z", intensity from <spectra_table> where id = $1 order by "m/z" asc), ST AS (select  "m/z", intensity from <spectra_table> where id = $2 order by "m/z" asc),
qc as (select count(1) as qp from sq), tc as (select count(1) as tp from st)
select coalesce(sum(sq.intensity*st.intensity)/sqrt(sum(sq.intensity^2)*sum(st.intensity^2)) ,0.0) as spectral_contrast from SQ,ST,qc,tc where qc.qp = tc.tp and SQ."m/z" = ST."m/z"'
LANGUAGE SQL
    IMMUTABLE
    RETURNS NULL ON NULL INPUT;

Please not the additional coalesce() to return 0.0 when the SQL evaluates to NULL because of the join.

In a real application, the calculation of S can of course be accelerated, e.g. by prefiltering spectra by their number of peaks or their m/z range.

If you want to see the angle in degrees, acosd() is your friend.

That's it for today.

Finding mass spectra with PostgreSQL: Tanimoto similarity in SQL

Now that we can store mass spectra in PostgreSQL, the next question is how to find them?

There is more than one solution to this problem. One is the Tanimoto coefficient or Jaccard index of the peaks.

The Tanimoto coefficient T can be calculated as T = c / (a+b-c) with c the number of peaks both spectra have in common and a and b the number of peaks in each spectrum.

Similarity is then expressed as a number between 0 and 1, with 0 meaning no and 1 maximum similarity. Likewise the Tanimoto distance 1-T expresses the distance between two spectra.

On a table like

CREATE TABLE <spectra_table>
(
  id integer NOT NULL,
  "m/z" numeric NOT NULL,
  intensity numeric NOT NULL
)

this can easily be done in SQL with a common table expression:

WITH SQ AS (select "m/z" from <spectra_table> where id=1),
ST AS (select "m/z" from <spectra_table> where id=2),
SC AS (select count(1)::float as common from SQ, ST where sq."m/z" = st."m/z")
select (select common FROM SC) / ((select count(1) from sq) + (select count(1) from st) - (SELECT common FROM SC))
as spectral_tanimoto;

If desired, further properties like the peak intensity can be used to narrow the match:

WITH SQ AS (select "m/z", intensity from <spectra_table> where id=1),
ST AS (select "m/z", intensity from <spectra_table> where id=2),
SC AS (select count(1)::float as common from SQ, ST where sq."m/z" = st."m/z" and sq.intensity = st.intensity)
select (select common FROM SC) / ((select count(1) from sq) + (select count(1) from st) - (SELECT common FROM SC))
as spectral_tanimoto;

The problem here is, that the join between the peaks relies on exact matching properties. 'Fuzzy' joins can relax this to some extent. The simple variant is to cast the m/z values to integers:


WITH SQ AS (select "m/z"::int from <spectra_table> where id=1),
ST AS (select "m/z"::int from <spectra_table> where id=2),
SC AS (select count(1)::float as common from SQ, ST where sq."m/z" = st."m/z")
select (select common FROM SC) / ((select count(1) from sq) + (select count(1) from st) - (SELECT common FROM SC))
as spectral_tanimoto;

PostgreSQL's range types allow finer grained control:

WITH SQ AS (select "m/z" from <spectra_table> where id=1),
ST AS (select "m/z" from <spectra_table> where id=2),
SC AS (select count(1)::float as common from SQ, ST where numrange(sq."m/z"-0.25, sq."m/z"+0.25,'[]') && numrange(st."m/z"-0.25, st."m/z"+0.25,'[]'))
select (select common FROM SC) / ((select count(1) from sq) + (select count(1) from st) - (SELECT common FROM SC))
as spectral_tanimoto;

Neat. But this comes with a serious caveat, since 'fuzzy' joins may produce excess tuples by their very nature. Given a table with two entries 1.1 and 0.9, joining with 1.1 and a range of +/- 0.2, now both tuples will match, generating two result tuples!

In the context of the Tanimoto coefficient, this can lead to the anomaly that the numerator c becomes larger than the denominator. If this happens, T > 1 or T < 0 may occur, which are strictly undefined. So handle with care.

The nicer properties of the calculation are, that you can easily transform spectra measured with different ionization levels and that it can be made run very fast, if you know the target T you like to filter with beforehand.

1.) Since an m/z of 50 measured at an ionization of +1 means a mass of 50 u but a m/z of 50 measured at ionization +2 means a mass of 100 u (m/z is the ratio between mass and charge of the Ion measured), the statements above would give false results if the ionization was different and not taken into account. This can be done at the join level:

...where sq."m/z"*2 = st."m/z"...
  
Usually the MS software does that for you so you have normalized spectra, but anyway...

2.) If you know the number of peaks in each spectrum and the target T you want to filter with, similarity searches can be made very fast. Here is the paper on how to do it.

That's it for today. Spectral contrast in SQL next.

oracle_fdw: 'cannot use LDAP' means 'don't even touch'

The README for the oracle_fdw clearly says: "You cannot use LDAP functionality both in PostgreSQL and in Oracle, period."

Actually this means: "DO NOT EVEN TOUCH!"

I found this out yesterday when tracking down frequent sigsevs in one of our PostgreSQL servers.

The Oracle instant client was set up with NAMES.DIRECTORY_PATH=(TNSNAMES, LDAP) in sqlnet.ora.

Even when the connect string used in the FDW is in the //<server>:<port>/<service> form, i.e. no Oracle name resolution is required, just allowing LDAP in sqlnet.ora will kill your server.

Calculating theoretical isotope patterns for mass spectrometry with PostgreSQL and pgchem::tigress

New year, new challenge: mass spectra.

To calculate the theoretical isotope pattern for a given molecule in PostgreSQL you need two things:

1. a cheminformatics toolkit to determine the element composition of the molecule and it's exact mass
2. a usably fast and accurate algorithm to infer the isotope pattern from the composition

For 1. I take OpenBabel, for 2. the Mercury7 algorithm as implemented in libmercury++.

Then, the two have to be made working together. Like so:

 

typedef struct
{
    size_t num_entries;
    double *mz;
    double *intensity;
    double *intensity_normalized;
    unsigned int *md;
} _ISOTOPE_PATTERN;

 

extern "C" _ISOTOPE_PATTERN *
ob_molfile_to_isotope_pattern(char *molfile, int charge, double normal)
{
    _ISOTOPE_PATTERN *retval;
    vector msa_mz;
    vector msa_abundance;
    vector composition;
    vector m;
    const double limit = 10e-30;
    OBMol mol;
    OBConversion conv;
    string tmpStr (molfile);
    istringstream molstream (tmpStr);
    double mass_diff, monoisotopic_mass, scale = 1.0, delta = numeric_limits::max();
    int i=0;

    //Initialize composition vector to all zeroes
    for(unsigned int i=0; i         composition.push_back(0);

    conv.SetInAndOutFormats("MDL",NULL);

    conv.Read(&mol, &molstream);

    if (mol.Empty ())
        return NULL;

    // Add missing hydrogens, if any
    mol.AddHydrogens(false, false); // All Hs, no PH correction

    //Iterate over atoms and update the composition vector accordingly
    FOR_ATOMS_OF_MOL(atom, mol)
    {
        switch (atom->GetAtomicNum())
        {
        case 1: // H
            composition[0]++;
            break;
        case 6: // C
            composition[1]++;
            break;
        case 7: // N
            composition[2]++;
            break;
        case 8: // O
            composition[3]++;
            break;
        case 16: // S
            composition[4]++;
            break;
        case 35: // Br
            composition[5]++;
            break;
        case 17: // Cl
            composition[6]++;
            break;
        case 53: // I
            composition[7]++;
            break;
        case 15: // P
            composition[8]++;
            break;
        case 9: // F
            composition[9]++;
            break;
        default: // Others ignored
            break;
        }
    }

    // Calculate isotope patterns with MERCURY7
    if(0 != mercury::mercury(msa_mz, msa_abundance, composition, charge, limit))
    {
        return NULL;
    }

    monoisotopic_mass = mol.GetExactMass();

    // Allocate return value and copy values
    retval = (_ISOTOPE_PATTERN*) calloc(1,sizeof(_ISOTOPE_PATTERN));

    retval->num_entries = msa_mz.size();
    retval->mz=(double*) calloc(msa_mz.size(), sizeof(double));
    retval->intensity=(double*) calloc(msa_abundance.size(), sizeof(double));
    retval->intensity_normalized=(double*) calloc(msa_abundance.size(), sizeof(double));
    retval->md=(unsigned int*) calloc(msa_abundance.size(), sizeof(unsigned int));

    copy(msa_mz.begin(), msa_mz.end(), retval->mz);
    copy(msa_abundance.begin(), msa_abundance.end(), retval->intensity);
    copy(msa_abundance.begin(), msa_abundance.end(), retval->intensity_normalized);

    for(std::vector::iterator it = msa_mz.begin(); it != msa_mz.end(); it++)
    {
        mass_diff = *it-monoisotopic_mass;

        m.push_back(mass_diff);

        if(abs(mass_diff) < delta)
        {
            delta = abs(mass_diff);
            scale = msa_abundance[i];
        }
        i++;
    }

    scale = normal / scale;

    for(i=0; retval->num_entries; i++)
    {
        retval->intensity_normalized[i] *= scale;
    }

    copy(m.begin(), m.end(), retval->md);

    return retval;
}
 

This parses a V2000 or V3000 molfile into a molecule object, iterates over the atoms and generates the elemental composition for Mercury7, calculates the exact mass and the theoretical isotope pattern. It then stores the m/z vector, the peak vector and the derived peak vector normalized to some value, e.g. 100%, and the mass distance relative to the base peak m/z value. All of this is stored in the C structure _ISOTOPE_PATTERN in order to transfer it from C++ to the C code of PostgreSQL.

 

Over to the PostgreSQL side of things :-)

 

PG_FUNCTION_INFO_V1 (pgchem_isotope_pattern);

Datum
pgchem_isotope_pattern (PG_FUNCTION_ARGS)
{
    _ISOTOPE_PATTERN *isotope_pattern = NULL;
    _ISOTOPE_PATTERN *pattern_copy = NULL;
    FuncCallContext *ctx=NULL;
    MemoryContext original_ctx=NULL;
    DECOMPRESSED_DATA *original_data=NULL;
    char *tmpMolfile=NULL;

    if(SRF_IS_FIRSTCALL())
    {
        MOLECULE *arg_molecule = PG_GETARG_MOLECULE_P (0);
        int32 arg_charge = PG_GETARG_INT32 (1);
        float8 arg_normal = PG_GETARG_FLOAT8 (2);
        ctx = SRF_FIRSTCALL_INIT();

        original_ctx = MemoryContextSwitchTo(ctx->multi_call_memory_ctx);

        if(arg_molecule->original_format==FORMAT_V2000 || arg_molecule->original_format==FORMAT_V3000)
        {
            original_data = decompress_data(CIPTR(arg_molecule),arg_molecule->compressed_sizeo, arg_molecule->sizeo);
            tmpMolfile = strdup(original_data->decompressed_data);
        }
        else
        {
            tmpMolfile = ob_smiles_to_V2000 (SMIPTR(arg_molecule));
        }

        isotope_pattern = ob_molfile_to_isotope_pattern(tmpMolfile, arg_charge, arg_normal);

        if(NULL!=tmpMolfile)
        {
            free(tmpMolfile);
        }

        if(NULL == isotope_pattern)
        {
            elog (ERROR, "Could not generate isotope pattern!");
            PG_RETURN_NULL();
        }

        pattern_copy = (_ISOTOPE_PATTERN*) palloc0(sizeof(_ISOTOPE_PATTERN));
        pattern_copy->mz = (double*) palloc0(isotope_pattern->num_entries*sizeof(double));
        pattern_copy->intensity = (double*) palloc0(isotope_pattern->num_entries*sizeof(double));
        pattern_copy->intensity_normalized = (double*) palloc0(isotope_pattern->num_entries*sizeof(double));
        pattern_copy->md = (unsigned int*) palloc0(isotope_pattern->num_entries*sizeof(unsigned int));

        pattern_copy->num_entries = isotope_pattern->num_entries;

        memcpy(pattern_copy->mz, isotope_pattern->mz, isotope_pattern->num_entries*sizeof(double));
        memcpy(pattern_copy->intensity, isotope_pattern->intensity, isotope_pattern->num_entries*sizeof(double));
        memcpy(pattern_copy->intensity_normalized, isotope_pattern->intensity_normalized, isotope_pattern->num_entries*sizeof(double));
        memcpy(pattern_copy->md, isotope_pattern->md, isotope_pattern->num_entries*sizeof(unsigned int));

        if(NULL != isotope_pattern->md)
            free(isotope_pattern->md);

        if(NULL != isotope_pattern->mz)
            free(isotope_pattern->mz);

        if(NULL != isotope_pattern->intensity)
            free(isotope_pattern->intensity);

        if(NULL != isotope_pattern->intensity_normalized)
            free(isotope_pattern->intensity_normalized);

        if(NULL != isotope_pattern)
            free(isotope_pattern);

        isotope_pattern = NULL;

        ctx->max_calls = pattern_copy->num_entries;
        ctx->user_fctx = pattern_copy;

        if(get_call_result_type(fcinfo, NULL, &ctx->tuple_desc) != TYPEFUNC_COMPOSITE)
        {
            elog (ERROR, "Calling the isotope pattern SRF in a context that does not expect tuples in return is not supported!");
        }

        BlessTupleDesc(ctx->tuple_desc);

        MemoryContextSwitchTo(original_ctx);
    }

    ctx = SRF_PERCALL_SETUP();

    isotope_pattern = ctx->user_fctx;

    if(ctx->call_cntr < ctx->max_calls)
    {
        HeapTuple rettuple = NULL;
        Datum *retvals;
        bool *retnulls;

        retvals = (Datum*) palloc0(4*sizeof(Datum));
        retnulls = (bool*) palloc0(4*sizeof(bool));

        retvals[0] = Float8GetDatum(isotope_pattern->mz[ctx->call_cntr]);
        retvals[1] = Float8GetDatum(isotope_pattern->intensity[ctx->call_cntr]);
        retvals[2] = Float8GetDatum(isotope_pattern->intensity_normalized[ctx->call_cntr]);
        retvals[3] = UInt8GetDatum(isotope_pattern->md[ctx->call_cntr]);

        retnulls[0] = false;
        retnulls[1] = false;
        retnulls[2] = false;
        retnulls[3] = false;

        rettuple = heap_form_tuple(ctx->tuple_desc, retvals, retnulls);

        pfree(retvals);
        pfree(retnulls);

        SRF_RETURN_NEXT(ctx, HeapTupleGetDatum(rettuple));
    }
    else
    {
        SRF_RETURN_DONE(ctx);
    }
}
 

This is a set returning function, meaning that it can return multiple rows in one call, i.e. it returns a complete table. Basically it retrieves the molecule argument from PostgreSQL, calls ob_molfile_to_isotope_pattern(), constructs a set of tuples and returns them to PostgreSQL. The basics of set returning functions can be found here and here.

 

Now we can declare the function to PostgreSQL:

 

CREATE OR REPLACE FUNCTION isotope_pattern(IN molecule, IN integer DEFAULT 0, IN double precision DEFAULT 100.0)
  RETURNS TABLE("m/z" double precision, intensity double precision, intensity_normalized double precision, mdtbp integer) AS
'libpgchem', 'pgchem_isotope_pattern'
  LANGUAGE c IMMUTABLE STRICT;
 

and use it:

select * from isotope_pattern('c1ccccc1I'::molecule) where intensity >= 0.01;
 

m/z	intensity	intensity_normalized	mdtbp
203.9436	0.934689933304048	100	0
204.946987674881	0.0634781699018538	6.79136124612619	1

 

  

Visualized it looks like this:

Nice. And everything (except the GUI) in the database.

Remember the Therac!

Personal opinion warning!!!

Every time I read something like this, I immediately think of the famous Therac-25.

But let's start from the bottom...

"If you are deploying technology in your healthcare organization, ask yourself first if your people are ready. If they are not, you are setting the stage for failure."

This is so generic, Mr. Crounse should win a a prize. I'm sure, there is a prize somewhere for phrases like that.

How about:

"If you are deploying technology in your organization, ask yourself first if your people are ready. If they are not, you are setting the stage for failure."

Almost universally true but the gain of knowledge is essentially zilch.

But what really irks me is stuff like this (on the training of physicians):

"We spend our young adulthood immersed in the scholarly pursuit of a medical degree. We take four or more years in specialty training. [...] It’s all about following a certain process, a definitive kind of workflow"

Call me oldschool, but I expect a physician to be well trained and following certain, proven processes, when I lay my fate into her hands!

IT has already a major role in healthcare. Without advanced signal processing and medical imaging, MRI would be useless. Genome data simply cannot be analyzed without the help of computers. The human body can be simulated to some extent to reduce the need of testing new drugs on animals and human beings alike. Particle therapy needs absolute precision controls, provided by computers. And so on and on...

And what do I read?

"They are embracing cloud technologies to streamline IT resources and focus more of those limited resources on that which healthcare systems are all about—providing care to patients and increasingly, focusing on ways to improve population health and disease prevention."

OK, people are expensive, care is expensive so let's use IT to squeeze more out of the staff we have.

There is no single sentence about medical progress in this article. It's from top to bottom about cost-cutting.

I have to give credit to the overall assumption, that something good could be gained  from the integration of health-care data, though. A lot of knowledge about how diseases work and what weapon works best to fight them with could be gained from a holistic view on all the data available.

Yet, in many countries, there is a strict regulation on medical data. And here I have to give credit again - it's all about people. People building systems, people using systems and people having trust in systems.

In 2010 and 2011, physicians manipulated medical data of their patients waiting for a liver transplant in order to push them higher up the waiting list. This so called Organspendeskandal led to the lowest trust and henceforth lowest number of potential organ donors in Germany since 2002.

I am a registered bone marrow donor at DKMS. Yet, we are only a scant 25 million worldwide.  Already way to few. Guess what would happen, if there was a major security breach in this database? Or in any other large health data database?

"Information is no candy bar" they say in Bombardiers, (hilarious and highly recommended) If you eat it the candy bar its gone, but information can be stolen without notice.

Let's assume somebody would steal all the paper patient files from a hospital, or the insurance company. This could be noticed, because a) they are gone and b) somebody might ask what the four trucks full of paper are all about. Electronic files can be stolen without notice on an USB stick.

We, the people, still cannot build reliable and secure software systems at a large scale. SQL Injection is still #1 at OWASP. As it was 2010. 2007 it was on #2. Unbelievable.

Even in highly specialized and controlled areas, we fail, fail, fail.

So, where does this strict regulation on medical data come from again?

Whenever you work on hard or software for healthcare applications - please remember the Therac.

So much useful knowledge could be gained, but so much trust can be lost.

Personal opinion warning!!!

A working RAMDisk for PostgreSQL on Windows

Aside from the debate if this is useful or dangerous, I tried several free RAMDisks for Windows, Gavotte RamDisk, AMD Radeon RAMDisk and SoftPerfect RAM Disk - and none of them works well with PostgreSQL.

Either PostgreSQL refuses to use them at all for datafiles or it seems to work, but at some point in time the server tells you, that there is something wrong with the datafiles on the ramdrive and data corruption bites.

Apparently most ramdrives on Windows provide drive emulations only good enough to support simple file storage.

ImDisk is the first I tried, that seems to work well with PostgreSQL.

It's free, has a GUI and a command line, and the driver is properly signed for x64 Windows.

There are people := {a,b} for a ≠ b IN (SELECT people FROM world): Those who understand SQL, and those who have to wait.

Couldn't resist. ;->

I just pushed a view on an Oracle Database from 464 seconds down to 3 seconds (~ 157x speedup or 0.0065% of the original execution time) by removing a few pointless ORDER BYs and replacing one expensive JOIN with subqueries.

ToroDB Adventures: Adding unimplemented commands

From

> db.createCollection("test")
{ "ok" : 0, "errmsg" : "Unimplemented command: create", "code" : 1000002 }

to

> db.createCollection("test")
{ "ok" : 1 }

in five simple steps:

1.) Find out how the command is called internally. toroDB tells you that in the message:

"Unimplemented command: create"

so it is "create".

2.) Look up the "create" command. Start in QueryCommandProcessor.java with the QueryCommandGroup enum which contains all enums of all known commands. Since it is an administration command, you'll find it in AdministrationQueryCommand.java.

3.)  "create" does nothing at the moment, so add some code.
I just copied it from "createIndexes" so

create,

becomes now

create {
        @Override
        public void doCall(RequestBaseMessage queryMessage, BSONDocument query, ProcessorCaller caller) throws Exception {
            caller.create(query);
        }
    },

4.) Now, ProcessorCaller needs to know about the new command. It's  enclosed in QueryCommandProcessor.java, so by adding

public void create(@Nonnull BSONDocument document) throws Exception {
queryCommandProcessor.create(document, messageReplier);
}

to ProcessorCaller AND the prototype

public void create(@Nonnull BSONDocument document, @Nonnull MessageReplier messageReplier) throws Exception;

to the enclosing QueryCommandProcessor class, make it known.

5.) Implement the command in ToroQueryCommandProcessor.java, which is the actual implementation of QueryCommandProcessor:

@Override
    public void create(BSONDocument document, MessageReplier messageReplier) {
      Map keyValues = new HashMap();
      keyValues.put("ok", MongoWP.OK);
       
      String collectionName = (String) document.getValue("create");
      Boolean capped = (Boolean) document.getValue("capped");
      Boolean autoIndexId = (Boolean) document.getValue("autoIndexId");
      Boolean usePowerOf2Sizes = (Boolean) document.getValue("usePowerOf2Sizes");
      Double size = (Double) document.getValue("size");
      Double max = (Double) document.getValue("max");
       
      BSONDocument reply = new MongoBSONDocument(keyValues);
      messageReplier.replyMessageNoCursor(reply);
    }

And that's pretty much it. As of now it just reads all allowed values from the command and acknowledges OK. But now everything is set to make it a "real" command if needed.

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.

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; 

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.