pandas

Creation of data objects

One of the most basic operations is the creation of a data object. Here are examples of how to create three different data objects: datarray.DataArray, pandas.DataMatrix, and la.larry.

First let's set up the data we will use to create the data objects:

>>> import numpy as np
>>> N = 100
>>> arr = np.random.rand(N, N)
>>> idx1 = map(str, range(N))
>>> idx2 = map(str, range(N))

Create a datarray.DataArray:

>>> import datarray
>>> datarray.__version__
    '0.0.4'
>>> timeit datarray.DataArray(arr, [('row', idx1), ('col', idx2)])
    10000 loops, best of 3: 153 us per loop

Create a pandas.DataMatrix:

>>> import pandas
>>> pandas.__version__
    '0.2.1.dev'
>>> timeit pandas.DataMatrix(arr, idx1, idx2)
    10000 loops, best of 3: 50.1 us per loop

Create a la.larry:

>>> import la
>>> la.__version__
    '0.5.0dev'
>>> timeit la.larry(arr, [idx1, idx2])
    100000 loops, best of 3: 14 us per loop

Most of the time used to create a larry is spent doing an integrity check of the labels. If you skip the integrity check, you get a 10x speed up:

>>> timeit la.larry(arr, [idx1, idx2], integrity=False)
    1000000 loops, best of 3: 1.33 us per loop

Tabular data, mutability

Many statistical data sets arrive in tabular (2D) form-- typically a collection of vectors with associated labels. R provides the mutable data.frame class:

> df <- data.frame(a=rnorm(5), b=rnorm(5))
> df
           a           b
1 -0.4454339 -0.87899052
2  0.1119211  0.08378914
3  0.5823613  1.51314787
4 -0.4235883  1.51682599
5 -0.4429238 -1.86833929
> df$a
[1] -0.4454339  0.1119211  0.5823613 -0.4235883 -0.4429238
> df$c <- 1
> df
           a           b c
1 -0.4454339 -0.87899052 1
2  0.1119211  0.08378914 1
3  0.5823613  1.51314787 1
4 -0.4235883  1.51682599 1
5 -0.4429238 -1.86833929 1

For users of other statistical environments, I'd be interested to hear about the tabular data structures and what they can do. Is such mutability useful? I've found it quite useful in cases like this. The R data.frame does not provide any sort of row-oriented data alignment (see below), but it's worth asking what are the features of such a 2D data structure that are the most useful.

The pandas code for the above is nearly identical, with pythonic notation of course:

In [22]: df = DataFrame({'a' : randn(5), 'b' : randn(5)})

In [23]: df
Out[23]:
     a              b
0    -1.894         -1.373
1    -0.8486        -0.5785
2    1.838          -3.761
3    -0.6144        -0.7812
4    0.8473         0.9544


In [24]: df['a']
Out[24]:
0    -1.89423416608
1    -0.848596359606
2    1.83799383312
3    -0.614350513146
4    0.84729188782

In [25]: df['c'] = 1

In [26]: df
Out[26]:
     a              b              c
0    -1.894         -1.373         1
1    -0.8486        -0.5785        1
2    1.838          -3.761         1
3    -0.6144        -0.7812        1
4    0.8473         0.9544         1

In the context of misaligned data, I've seen a lot of pandas users use DataFrame purely as a container for aligning labeled data:

In [74]: df
Out[74]:
                       A              B
2000-01-03 00:00:00    -0.2306        0.3456
2000-01-04 00:00:00    -0.882         -0.1626
2000-01-05 00:00:00    0.2474         -0.1887
2000-01-06 00:00:00    -0.2704        0.1398
2000-01-07 00:00:00    1.498          1.628
2000-01-10 00:00:00    -1.645         0.09256
2000-01-11 00:00:00    0.5666         -1.128
2000-01-12 00:00:00    -0.5082        -0.5381
2000-01-13 00:00:00    0.2205         -0.5728
2000-01-14 00:00:00    -0.4916        -0.3333


In [76]: s
Out[76]:
2000-01-03 00:00:00    -0.115323455775
2000-01-04 00:00:00    -0.441017336819
2000-01-05 00:00:00    0.123677568341
2000-01-06 00:00:00    -0.135191054412
2000-01-07 00:00:00    0.749180511001
2000-01-10 00:00:00    -0.822250171143
2000-01-11 00:00:00    0.283300117597
2000-01-12 00:00:00    -0.254115521874

In [77]: df['C'] = s

In [78]: df
Out[78]:
                       A              B              C
2000-01-03 00:00:00    -0.2306        0.3456         -0.1153
2000-01-04 00:00:00    -0.882         -0.1626        -0.441
2000-01-05 00:00:00    0.2474         -0.1887        0.1237
2000-01-06 00:00:00    -0.2704        0.1398         -0.1352
2000-01-07 00:00:00    1.498          1.628          0.7492
2000-01-10 00:00:00    -1.645         0.09256        -0.8223
2000-01-11 00:00:00    0.5666         -1.128         0.2833
2000-01-12 00:00:00    -0.5082        -0.5381        -0.2541
2000-01-13 00:00:00    0.2205         -0.5728        NaN
2000-01-14 00:00:00    -0.4916        -0.3333        NaN

To combine two larrys, use merge:

>> lar1 = larry(randn(5,2))
>> lar2 = larry(randn(2,1), [[0, 4], ['c']])

>> lar1.merge(lar2)
label_0
    0
    1
    2
    3
    4
label_1
    0
    1
    c
x
array([[ 0.19271985,  0.38086157,  0.89416946],
       [ 0.17987281, -0.20143654,         NaN],
       [ 0.56567726,  2.13787551,         NaN],
       [-1.92165255,  0.22968882,         NaN],
       [ 1.14440607,  0.07242789, -0.56150006]])

If the two larrys overlap, you can use the update option:

>> lar3 = larry(np.ones((2,1)), [[0, 4], [0]])
>> lar1.merge(lar3, update=True)
label_0
    0
    1
    2
    3
    4
label_1
    0
    1
x
array([[ 1.        ,  0.38086157],
       [ 0.17987281, -0.20143654],
       [ 0.56567726,  2.13787551],
       [-1.92165255,  0.22968882],
       [ 1.        ,  0.07242789]])

Data alignment

• Union versus intersection of data labels
• Performance for differently-indexed data

Here is a pretty pathological but realistic example:

# Setup
import numpy as np
import pandas
import la
N = 10000
arr1 = np.random.randn(N)
arr2 = np.random.randn(N)
idx1 = range(N)
idx2 = range(N)[::-1]

# pandas
series1 = pandas.Series(arr1, index=idx1)
series2 = pandas.Series(arr2, index=idx2)

# larry
lar1 = la.larry(arr1, [idx1])
lar2 = la.larry(arr2, [idx2])

Here is with latest svn/bzr trunks on 5/30/10:

# aligned
In [5]: timeit series1 + serie1
10000 loops, best of 3: 36.4 us per loop
In [6]: timeit lar1 + lar1
10000 loops, best of 3: 72.5 us per loop

# unaligned
In [9]: timeit series1 + series2
100 loops, best of 3: 2.72 ms per loop
In [8]: timeit lar1 + lar2
100 loops, best of 3: 4.97 ms per loop

A 2d example:

# Setup
import numpy as np
import pandas
import la
N = 1000
arr1 = np.random.randn(N, N)
arr2 = np.random.randn(N, N)
idx1 = range(N)
idx2 = range(N)[::-1]

# pandas
dma1 = pandas.DataMatrix(arr1, idx1, idx1)
dma2 = pandas.DataMatrix(arr2, idx2, idx2)

# larry
lar1 = la.larry(arr1, [idx1, idx1])
lar2 = la.larry(arr2, [idx2, idx2])

Here is with latest svn/bzr trunks on 5/30/10:

# aligned
>> timeit dma1 + dma1
100 loops, best of 3: 16 ms per loop
>> timeit lar1 + lar1
100 loops, best of 3: 15.9 ms per loop

# unaligned
>> timeit dma1 + dma2
10 loops, best of 3: 93.5 ms per loop
>> timeit lar1 + lar2
10 loops, best of 3: 91.7 ms per loop

Reduction operations

Find the maximum value in a 2d data object of shape (1000, 1000):

import pandas
import la

N = 1000
arr = np.random.rand(N, N)
dma = pandas.DataMatrix(arr, range(N), range(N))
lar = la.larry(arr)

# max

dma.max().max()
0.99999943789619117
lar.max()
0.99999943789619117

>> timeit dma.max().max()
10 loops, best of 3: 18 ms per loop
>> timeit lar.max()
100 loops, best of 3: 14.1 ms per loop
>> timeit np.nanmax(arr)
100 loops, best of 3: 14.1 ms per loop

IO

Users of data objects often need the ability to quickly archive their data. Here is the time taken to save and load two data objects of shape (N, N):

Could someone run this benchmark on a Linux machine with NumPy > 1.3? I (Wes) ran these on Windows and I don't trust the OS.

In [8]: roundtrip_archive(100)
Numpy (npz)    0.0107 seconds
larry (HDF5)   0.0064 seconds
pandas (HDF5)  0.0050 seconds

In [9]: roundtrip_archive(500)
Numpy (npz)    0.1026 seconds
larry (HDF5)   0.0908 seconds
pandas (HDF5)  0.0836 seconds

In [10]: roundtrip_archive(1000)
Numpy (npz)    0.5945 seconds
larry (HDF5)   0.3517 seconds
pandas (HDF5)  0.2611 seconds

In [11]: roundtrip_archive(2000, iterations=2)
Numpy (npz)    4.6487 seconds
larry (HDF5)   0.8843 seconds
pandas (HDF5)  0.7653 seconds

In [12]: roundtrip_archive(4000, iterations=1)
Numpy (npz)   19.7412 seconds
larry (HDF5)   2.9064 seconds
pandas (HDF5)  3.3371 seconds

And here is the code used:

This needs pandas r180 and PyTables.

import time, os
import numpy as np

import la
import pandas

def timeit(f, iterations):
    start = time.clock()

    for i in xrange(iterations):
        f()

    return time.clock() - start

def roundtrip_archive(N, iterations=10):

    # Create data
    arr = np.random.randn(N, N)
    lar = la.larry(arr)
    dma = pandas.DataMatrix(arr, range(N), range(N))

    # filenames
    filename_numpy = 'c:/temp/numpy.npz'
    filename_larry = 'c:/temp/archive.hdf5'
    filename_pandas = 'c:/temp/pandas_tmp'

    # Delete old files
    try:
        os.unlink(filename_numpy)
    except:
        pass
    try:
        os.unlink(filename_larry)
    except:
        pass
    try:
        os.unlink(filename_pandas)
    except:
        pass

    # Time a round trip save and load
    numpy_f = lambda: numpy_roundtrip(filename_numpy, arr, arr)
    numpy_time = timeit(numpy_f, iterations) / iterations

    larry_f = lambda: larry_roundtrip(filename_larry, lar, lar)
    larry_time = timeit(larry_f, iterations) / iterations

    pandas_f = lambda: pandas_roundtrip(filename_pandas, dma, dma)
    pandas_time = timeit(pandas_f, iterations) / iterations

    print 'Numpy (npz)   %7.4f seconds' % numpy_time
    print 'larry (HDF5)  %7.4f seconds' % larry_time
    print 'pandas (HDF5) %7.4f seconds' % pandas_time

def numpy_roundtrip(filename, arr1, arr2):
    np.savez(filename, arr1=arr1, arr2=arr2)
    npz = np.load(filename)
    arr1 = npz['arr1']
    arr2 = npz['arr2']

def larry_roundtrip(filename, lar1, lar2):
    io = la.IO(filename)
    io['lar1'] = lar1
    io['lar2'] = lar2
    lar1 = io['lar1']
    lar2 = io['lar2']

def pandas_roundtrip(filename, dma1, dma2):
    from pandas.io.pytables import HDFStore
    store = HDFStore(filename)
    store['dma1'] = dma1
    store['dma2'] = dma2
    dma1 = store['dma1']
    dma2 = store['dma2']

def pandas_roundtrip_pickle(filename, dma1, dma2):
    dma1.save(filename)
    dma1 = pandas.DataMatrix.load(filename)
    dma2.save(filename)
    dma2 = pandas.DataMatrix.load(filename)

Sometimes you only want to load a small part of a large data object from your archive. larry supports loading pieces of a data object. Here are some examples:

io['lar2'][:,-1]
io['lar2'][0,:]

Joining / merging data sets

In pandas, it's very similar to elsewhere

In [65]: df1
Out[65]:
                       A              B
2000-01-03 00:00:00    -0.1174        -0.941
2000-01-04 00:00:00    -0.6034        -0.008094
2000-01-05 00:00:00    -0.3816        -0.9338
2000-01-06 00:00:00    -0.3298        -0.9548
2000-01-07 00:00:00    0.9576         0.4652
2000-01-10 00:00:00    -0.7208        -1.131
2000-01-11 00:00:00    1.568          0.8498
2000-01-12 00:00:00    0.3717         -0.2323
2000-01-13 00:00:00    -1.428         -1.997
2000-01-14 00:00:00    -1.084         -0.271

In [67]: df2
Out[67]:
                       C              D
2000-01-03 00:00:00    0.2833         -0.1937
2000-01-05 00:00:00    1.868          1.207
2000-01-07 00:00:00    -0.8586        -0.7367
2000-01-11 00:00:00    2.121          0.9104
2000-01-13 00:00:00    0.7856         0.9063


In [68]: df1.join(df2)
Out[68]:
                       A              B              C              D
2000-01-03 00:00:00    -0.1174        -0.941         0.2833         -0.1937
2000-01-04 00:00:00    -0.6034        -0.008094      NaN            NaN
2000-01-05 00:00:00    -0.3816        -0.9338        1.868          1.207
2000-01-06 00:00:00    -0.3298        -0.9548        NaN            NaN
2000-01-07 00:00:00    0.9576         0.4652         -0.8586        -0.7367
2000-01-10 00:00:00    -0.7208        -1.131         NaN            NaN
2000-01-11 00:00:00    1.568          0.8498         2.121          0.9104
2000-01-12 00:00:00    0.3717         -0.2323        NaN            NaN
2000-01-13 00:00:00    -1.428         -1.997         0.7856         0.9063
2000-01-14 00:00:00    -1.084         -0.271         NaN            NaN

In [70]: df1.join(df2, how='inner')
Out[70]:
                       A              B              C              D
2000-01-03 00:00:00    -0.1174        -0.941         0.2833         -0.1937
2000-01-05 00:00:00    -0.3816        -0.9338        1.868          1.207
2000-01-07 00:00:00    0.9576         0.4652         -0.8586        -0.7367
2000-01-11 00:00:00    1.568          0.8498         2.121          0.9104
2000-01-13 00:00:00    -1.428         -1.997         0.7856         0.9063

The index (row labels) are the default key for joining, but a column can also be used for a similar SQL-like join:

In [73]: df2
Out[73]:
                       C              D              key
2000-01-03 00:00:00    0.2833         -0.1937        0
2000-01-05 00:00:00    1.868          1.207          1
2000-01-07 00:00:00    -0.8586        -0.7367        0
2000-01-11 00:00:00    2.121          0.9104         1
2000-01-13 00:00:00    0.7856         0.9063         0


In [74]: df3 = DataFrame({'code' : {0 : 'foo', 1 : 'bar'}})

In [75]: df3
Out[75]:
     code
0    foo
1    bar


In [76]: df2.join(df3, on='key')
Out[76]:
                       C              D              code           key
2000-01-03 00:00:00    0.2833         -0.1937        foo            0
2000-01-05 00:00:00    1.868          1.207          bar            1
2000-01-07 00:00:00    -0.8586        -0.7367        foo            0
2000-01-11 00:00:00    2.121          0.9104         bar            1
2000-01-13 00:00:00    0.7856         0.9063         foo            0

Groupby-type operations

Time series / date functionality

Handling missing data

Statistics in computational biology and bioinformatics

I am writing here in response to Vincent Davis' message on the Biopython mailing list about feedback on statistical data structure design issues in biopython. While I have been a Biopython developer for many years, different users have different needs and perspectives, so what I write here is just my personal opinion rather than an official position of the Biopython consortium.

With respect to Python/NumPy/SciPy, computational biology and bioinformatics is different from other numerical sciences in three important ways:

• Bioinformatics is dominated by Perl, and therefore people new to bioinformatics tend to choose Perl as a default rather than Python. Many people do their statistical analysis by a combination of Perl and R. For us in Biopython this means that we should make using Python as straightforward as possible.
• People in bioinformatics come from a wide range of educational backgrounds. Some come from hard-core computational sciences such as physics, mathematics, or computer science, whereas others have a more traditional background in molecular biology and may not have a solid background in numerics. Some people starting in bioinformatics do not know what a compiler is (but otherwise they can be great scientists). Unfortunately, this means that tools available in SciPy are often not accessible simply because it is too difficult to install.
• Some numerical libraries may be more important in bioinformatics than in other fields and vice versa. For example, to be able to calculate p-values it is important to have a library of special functions, whereas linear algebra may be less important.

What I really miss in NumPy are special functions and some basic statistical functions such as kernel density estimators. I know that these are exist in SciPy, but because of the second point this won't do me much good. I hope that some day at least the special functions can be made available in NumPy instead of SciPy.

See below for a discussion of MaskedArrays (numpy.ma)