Python for Bioinformatics

MathJax

2012-01-22T13:51:00.001-05:00

Let's follow up on a few previous posts (a calculation of d/dt of the total energy in a gravitational field, a post about typesetting math, and an introductory exploration of LaTeX). I'm grateful for a comment from a reader, with a link to this page, which makes it pretty clear that, yes, you can use their server to get the script in your pages.

So that's what we'll do. I had to modify the LaTeX commands a bit for this but it's mostly the same as before.

[ UPDATE: I do see a problem, now. The script applies to all LaTeX on the page, which if you go to the main page for the blog, includes the previous post... Just click on the post title, to see the original formatting code. ]

$\mathbf{F} = m\mathbf{a}= m\ddot{\mathbf{r}}
= m\frac{d^2}{dt^2}\mathbf{r}$

$\mathbf{F} = -$$\nabla$$V(\mathbf{r})$

$E = \frac{1}{2}m|\dot{\mathbf{r}}|^2+ V$

$\frac{d}{dt} E = \frac{d}{dt} ( \frac{1}{2}m
|\dot{\mathbf{r}}|^2 + V )$ $= ?$

$|\dot{\mathbf{r}}|^2 = |\dot{\mathbf{r}}|
|\dot{\mathbf{r}}| = \dot{\mathbf{r}}
\cdot\dot{\mathbf{r}}$

$\frac{d}{dt}\frac{1}{2}m|\dot{\mathbf{r}}|^2
= \frac{1}{2}m\frac{d}{dt}(\dot{\mathbf{r}}
\cdot\dot{\mathbf{r}})$$ = m\dot{\mathbf{r}}\cdot\ddot{\mathbf{r}}$$ = \dot{\mathbf{r}}\cdot-($$\nabla$$V$$)$

$\nabla$$V$ $= < \frac{\partial{V}} {\partial{x}},\frac{\partial{V}} {\partial{y}},\frac{\partial{V}} {\partial{z}} >$

$\dot{\mathbf{r}}$ = $<$ $\frac{dx}{dt}$,$\frac{dy}{dt}$,$\frac{dz}{dt}$$ >$

$\nabla$$V$ $\cdot$$\dot{\mathbf{r}}$ $= <\frac{\partial{V}} {\partial{x}} \frac{dx}{dt},\frac{\partial{V}} {\partial{y}} \frac{dy}{dt},\frac{\partial{V}} {\partial{z}}\frac{dz}{dt}>$$=\frac{d}{dt}V$

$\frac{d}{dt} E = \frac{d}{dt} ( \frac{1}{2}m
|\dot{\mathbf{r}}|^2 + V )$ $= ?$

$\frac{d}{dt}$E $= \dot{\mathbf{r}}$$\cdot$(-$\nabla$V) + $\nabla$V $\cdot$ $\dot{\mathbf{r}}$ = 0

Typesetting, again

2012-01-21T19:37:00.003-05:00

Here are the equations from the last post, in LaTeX. It took about two hours, fairly painstaking work, but I think it looks pretty good..

UPDATE:

Here is a nice simple reference for LaTeX. I started with "LaTeX Bootcamp" (pdf).

$\mathbf{F} = m\mathbf{a}= m\ddot{\mathbf{r}} = m\frac{d^2}{dt^2}\mathbf{r}$ $\mathbf{F} = -$$\nabla$$V(\mathbf{r})$ $E = \frac{1}{2}m\norm{\dot{\mathbf{r}}}^2+ V$ $\frac{d}{dt} E = \frac{d}{dt} (\frac{1}{2}m \norm{\dot{\mathbf{r}}}^2 + V)$ \hspace*{2em}$ = ?$\\ $\norm{\dot{\mathbf{r}}}^2 = \norm{\dot{\mathbf{r}}} \norm{\dot{\mathbf{r}}} = \dot{\mathbf{r}} \cdot\dot{\mathbf{r}}$ $\frac{d}{dt}\frac{1}{2}m\norm{\dot{\mathbf{r}}}^2 = \frac{1}{2}m\frac{d}{dt}(\dot{\mathbf{r}} \cdot\dot{\mathbf{r}})$ \hspace*{2em}$ = m\dot{\mathbf{r}}\cdot \ddot{\mathbf{r}}$ \hspace*{2em}$ = \dot{\mathbf{r}}\cdot-($ $\nabla$$V$$)$\\ $\nabla$$V$$ = <\frac{\partial{V}} {\partial{x}},\frac{\partial{V}} {\partial{y}},\frac{\partial{V}} {\partial{z}}>$ $\dot{\mathbf{r}}$ = $<$$\frac{dx}{dt}$,$\frac{dy}{dt}$,$\frac{dz}{dt}$$>$ $\nabla$$V$$\cdot$$\dot{\mathbf{r}}$$=<\frac{\partial{V}} {\partial{x}} \frac{dx}{dt},\frac{\partial{V}} {\partial{y}} \frac{dy}{dt},\frac{\partial{V}} {\partial{z}}\frac{dz}{dt}>$ \hspace*{2em}$=\frac{d}{dt}V$\\ $\frac{d}{dt}$E = $\dot{\mathbf{r}}$$\cdot$(- $\nabla$V ) + $\nabla$V $\cdot$ $\dot{\mathbf{r}}$ = 0

Vector fun

2012-01-21T10:37:00.001-05:00

In Marsden & Tromba's Vector Calculus, I found the following problem. It involves time-derivatives of the position vector r(t), for which we're using the dot notation of the physicists (and Newton).

In the first panel, we have Newton's second law, then the statement that the force is minus the gradient of the gravitational potential V, followed by a calculation of the total energy in the system. The problem is to compute d/dt of the energy. (hint)

The answer is given in the last panel, below. According to the book, it's a "simple calculation."

The first step is to compute d/dt of the kinetic energy. We use the formula from above, plus a trick to convert the squared term back to the dot product of dr/dt with itself. Then we use the chain rule, and finally, the definition of the force in terms of the gradient of the potential.

Calculation of d/dt of the potential energy puzzled me for quite a while, though it really shouldn't have. My solution was to work backward from the answer, as shown.

We put the two results together, and notice that they cancel. Surprise!

Typesetting math

2012-01-21T08:19:00.000-05:00

This post is really just a note to myself about something I need to investigate more. The problem I'd like to solve is how to format mathematical equations for the web. You can see the caveman approach on most of my posts here, like this one.

I just make a table in html

<div style="overflow-x: scroll;">
<table bgcolor="#ffffb0" border="0" padding="4">
<tbody>
<tr><td><pre style="font-family: monaco;">

I like the background color, and vary it depending on whether the content is code or output from a program.

It might be nice to have something prettier. So, looking around, I happened across MathJax, and also this post which explains how to use it on Blogger.

Here's a screenshot of the example:

The method used in the post is to load the script from the mathjax server, but I think what I'm probably supposed to do is direct people to resource on my (nonexistent) server. OTOH, they link to the post on the mathjax site.

What I'll probably do is just look into how to use LateX and then post screenshots.

Any thoughts?

Cycloid

2012-01-21T07:15:00.000-05:00

While working through the MIT ocw lectures on multi-variable calculus (Prof. Denis Auroux, here), I particularly enjoyed his discussion about the cycloid. Above is a graphic from the wikipedia article (actually the graphic is an animated gif, but I grabbed one of the frames). The red curve is generated by the motion of a point on the edge of a rolling circle.

In addition to the beauty of the curve, it turns out that the length and area under the curve have simple values that are relatively easy to calculate. See wikipedia for the details.

One thing the article doesn't explain is how to get the "parametrization" for the curve. This looks hard, but is made easy by using vectors. It's explained in the second half of Auroux's fifth lecture.

Another thing the article doesn't explain is how to integrate

√(2 - 2 cos t)

Start from the double angle formula:

cos 2s = cos²s - sin²s
cos 2s = 1 - 2 sin²s
2 sin²s = (1 - cos 2s)
2 sin²(t/2) = (1 - cos t)

It's straightforward from there.

The Mathworld article is also quite nice, and references a famous challenge in history, the one which led to this quote (in reference to Newton):

"Ah, I know the lion by his paw!"

Pilgrim's progress

2012-01-10T18:05:00.000-05:00

Lately I've been working (again) on understanding multi-variable calculus. I always wished I had time to cover this material in college, but I took lots of courses in biochemistry and molecular biology instead. I wanted to put up the best resources I've found for this so far.

First, and no surprise, the ocw videos at MIT by Denis Auroux (here). I am currently at #25, about to jump off into integrals in 3D space. (It's my second time to get this far). This time I think I really have everything under control. The proof of Green's theorem was perfect, a kind of mathematical satori.

I also work through everything in Gilbert Strang's book as well. I need to read carefully there, because he's so concise, but his insight is just incredible. Check out the proof of the flux version of Green's theorem.

Two more are recently discovered resources that are especially helpful because they develop everything slowly but completely:

Paul's lecture notes
A beautiful set of pages from someone at U Minnesota

All highly recommended.

Getting a little help from Newton

2011-12-14T15:43:00.001-05:00

I've posted here on a number of topics and derivations in fundamental calculus. This is about finding the derivative for powers of x with negative or fractional exponents. It seems a trivial proof, until I got stopped in my tracks by the cofactors of the binomial distribution.

The idea is to find the slope of the curve at x by adding a little bit to x, (x + h), calculating f(x + h), then subtract f(x), and finally divide by h and finding the limit as h gets very small. It works great for y = x². But consider

y = x^1/2

We want to do:

1/h [(x + h)^1/2 - x^1/2]

So thinking about the binomial expansion of the x + h term, I'm trying to figure what is the expansion for the 1/2 power? I mean, what can "n choose k" mean when n is 1/2?

The other thing that is a good approach here is to bring x out of the sum as follows

(x + h)^1/2
= x^1/2(1 + h/x)^1/2

To expand this, we just remember Newton's super-duper binomial expansion from here. If the expression is

(1 + Q)^m/n

the expansion is

1 + (m/n)Q + (m/n)((m-n)/2n)Q² + (m/n)((m-n)/2n)((m-2n)/3n)Q³ + ..

We only need the first two terms because Q² is (h/x)² which we will ignore. So we have for the expansion

(1 + h/x)^1/2
= 1 + 1/2(h/x)

and the derivative is the limit as h => 0 of

1/h [x^1/2(1 + 1/2(h/x)) - x^1/2]
= 1/h (1/2 h) 1/x^1/2
= 1/2 x^-1/2

Similarly for m = -2 and n = 1 the first two terms are 1 and -2 and the expansion we need will be

(1 + (h/x))^-2 = 
= 1 - 2(h/x)

The derivative is the limit as h => 0 of

1/h [x^-2 (1 - 2(h/x)) - x^-2] =
= 1/h x^-2 (-2)(h/x)
= -2 x^-3

Easy when you know how.

Basic Python for Bioinformatics

2011-12-14T14:39:00.000-05:00

I've been silent online for the last few months, due to other commitments associated with my day job. But I haven't forgotten all of you in the ether.

Just before I got busy, I wrote a book about basic programming with Python using lots of biological examples. It's been proofed once or twice, but is not quite a beta version, so let's call it the gamma release. I want to push it out now, even it could certainly use some more polish, because it's not clear when I'll get the time to do that.

The target audience is the novice programmer who is interested in biology. However, it differs from most similar material because there are lots of biological and sequence analysis examples, collected from various sources including posts on this site.

Here is a screenshot of the toc for the middle section of the book.

The book is in html format (link to Dropbox zip). It was built using the Sphinx software. One can generate a pdf version of the book as well, but it suffers from the major limitation that cut-and-paste of the code loses the indentation. The version linked to here has code examples that can be pasted directly into a text file and executed, and often will work in the interpreter as well.

I hope you like it. I would be very grateful for comments, criticisms, corrections, and even complaints. Post 'em as comments here, or find my work email through the aboutme. Enjoy!

Binary multiplication, again

2011-08-31T18:23:00.000-04:00

I had a lot of fun this afternoon, but all along I knew I was "reinventing the wheel." It's easy when you know where you're going.

I implemented binary multiplication with a class that keeps its data as Python strings like '00100011'. It's probably because I have been re-reading Charles Petzold's book Code, which I think is a classic. In any event, I wrote the class "word" which allows objects to add, multiply and exponentiate themselves. The second file is test.py which exercises the class.

It's a great "simple Python" coding project. We do bit-shifting by hand. If you want to be slick, you could go further by trying to implement fast exponentiation. Here is another example by a guy who obviously knows what he's doing.

A bit of sample output:

> python test.py

x     467

y     327

x + y 00000000 00000000 00000011 00011010 = 794

check 467 + 327 = 794

x * y 00000000 00000010 01010100 10000101 = 152709

check 467 * 327 = 152709

z      3

x**z  00000110 00010010 00010010 00001011 = 101847563

check 467 ** 3 = 101847563



x     995

y     402

x + y 00000000 00000000 00000101 01110101 = 1397

check 995 + 402 = 1397

x * y 00000000 00000110 00011010 01110110 = 399990

check 995 * 402 = 399990

z      3

x**z  00111010 10110111 00001100 10111011 = 985074875

check 995 ** 3 = 985074875



x     897

y     490

x + y 00000000 00000000 00000101 01101011 = 1387

check 897 + 490 = 1387

x * y 00000000 00000110 10110100 11101010 = 439530

check 897 * 490 = 439530

z      4

x**z  10111011 11001001 10001110 00000001 = 3150548481

ran into a slight problem

word.py

class word:

    SZ = 32

    def __init__(self,n):

        if type(n) == type('a'):

            b = n.rjust(word.SZ,'0')

        elif type(n) == type(1):

            b = bin(n)[2:]

        else:

            raise ValueError, "can't do that"

        self.b = b.rjust(word.SZ,'0')

        

    def __repr__(self):

        b = self.b

        s = ' '.join([b[:8], b[8:16], b[16:24], b[24:]])

        s += ' = ' + str(eval('0b' + self.b))

        return s

        

    def __add__(self, a):

        carry = '0'

        rL = list()

        for i in range(len(self.b)-1,-1,-1):

            x = self.b[i]

            y = a.b[i]

            t = x+y

            if (t == '01' or t == '10'):

                if carry == '0':

                    rL.append('1')

                if carry == '1':

                    rL.append('0')

            elif t == '00':

                rL.append(carry)

                carry = '0'

            else:

                assert t == '11'

                if carry == '0':

                    rL.append('0')

                    carry = '1'

                else:

                    rL.append('1')

                    carry = '1'

        if carry == '1':

            raise ValueError, 'overflow'

        rL.reverse()

        return word(''.join(rL))

    

    def __mul__(self, a):

        L = list()

        for i in range(len(self.b)-1,-1,-1):

            x = self.b[i]

            if not x == '1':

                continue

            n = word.SZ - i - 1

            r = word(a.b[n:] + '0'*(n))

            L.append(r)

        res = L.pop(0)

        #print 'start:       ', res

        while L:

            next = L.pop(0)

            #print 'next:        ', next

            res = res + next

            #if L:

                #print 'intermediate:', 

            #else:

                #print 'result:      ',

            #print res

        return res



    def __pow__(self, n):

        res = self

        for i in range(n-1):

            res = res * self

        return res

test.py

import random

from word import word

        

R = range(1000)

N = 10

for i in range(N):

    x = random.choice(R)

    y = random.choice(R)    

    xw = word(x)

    print 'x    ', x

    yw = word(y)

    print 'y    ', y

    

    print 'x + y',

    r = xw + yw

    print r

    S = x+y

    assert S == int(str(r).split()[-1])

    print 'check', x, '+', y, '=', S

    

    print 'x * y',

    r = xw * yw

    print r

    P = x*y

    assert P == int(str(r).split()[-1])

    print 'check', x, '*', y, '=', P

    

    for i in range(2,N*2):

        try:

            r = xw**i

        except ValueError:

            i = i-1

            break

    

    z = i

    print 'z     ', z

    print 'x**z ', r

    E = x**z

    try:

        assert E == int(str(r).split()[-1])

    except AssertionError:

        print 'ran into a slight problem\n'

        continue

    print 'check', x,'**', z, '=', E

    print

Dissecting RSA keys in Python (4)

2011-08-31T15:04:00.000-04:00

Still working with RSA keys. I accomplished the most important of the items left undone so far, to use the keys to encrypt a message using my own code.

In the process, I learned a new simple Python fact: the pow function (which is now a built-in), can take a modulus as a third argument.

And, it works where this doesn't seem to: (x**y) % z

(I discovered this by digging into the rsa module source).

python -m timeit -s "import rsa;  \

f = open('id_rsa');  data = f.read(); f.close(); k = rsa.PrivateKey.load_pkcs1(data)"\

 "41330915578951772302369**k.e % k.n"



10000 loops, best of 3: 79.4 usec per loop



python -m timeit -s "import rsa;  \

f = open('id_rsa');  data = f.read(); f.close(); k = rsa.PrivateKey.load_pkcs1(data)"\

 "pow(41330915578951772302369, k.e, k.n)"



10000 loops, best of 3: 70.2 usec per loop

In the example here, both work and have about the same timing. But in the code below, the first version hangs when doing decryption (with a large base). Here's the output, followed by the script.

> python encode.py 

m:   Hello, secret world!

p:   .xyz.Hello, secret world!.xyz.

a:   413309155789517723023698766343791993289928631329

c:   11434905702482726455415220687715293368262190253795 ..

i:   413309155789517723023698766343791993289928631329

r:   Hello, secret world!

encode.py

import rsa



with open('id_rsa') as f:

    data = f.read()

k = rsa.PrivateKey.load_pkcs1(data)

n = k.n

e = k.e

d = k.d



def my_atoi(s):

    L = [ord(c) for c in s]

    k = 256

    iL = L[:]

    iL.reverse()

    x = iL[0]

    for i in iL[1:]:

        x += i*k

        k *= 256

    return x



def my_itoa(i):

    rL = list()

    while i:

        rL.append(i%256)

        i = i/256

    rL.reverse()

    return ''.join([chr(n) for n in rL])



def encrypt(m):

    return m**e % n

    

def decrypt(c):

    # note:  c**d % n fails

    return pow(c,d,n)

    

if __name__ == '__main__':

    m = 'Hello, secret world!'

    pad = '.xyz.'

    p = pad + m + pad

    a = my_atoi(m)

    c = encrypt(a)

    i = decrypt(c)

    r = my_itoa(i)

    r = r.replace(pad,'')



    L = zip('mpacir',[m,p,a,c,i,r])

    N = 50

    for varname, var in L:

        s = str(var)

        print varname + ':  ', s[:N],

        if len(s) > N:  print '..'

        else:  print

Dissecting RSA keys in Python (3)

2011-08-31T08:00:00.000-04:00

I'm working on understanding the structure of RSA keys as output by the ssh-keygen utility (previous posts here and here). At the risk of becoming redundant, I want to try to simplify what we did previously. Let's begin by rehashing old material on bytes in Python:

It's easy to get confused between a byte and its string representation, or for that matter, between an unsigned int and its string representation. Although hexadecimal is natural too, it's perhaps easiest to think of an individual byte in terms of the decimal equivalent (0..255). That's because a Python int converts easily to chr, bin, and hex.

>>> i = 35

>>> c = chr(i)

>>> c

'#'

>>> i == ord(c)

True

>>> b = bin(i)

>>> b

'0b100011'

>>> h = hex(i)

>>> h

'0x23'

Notice that hex and bin both give unpadded output:

>>> hex(1)

'0x1'

>>> bin(1)

'0b1'

Don't be fooled. h and b here are 'str' datatypes. But these string representations of bin and hex values go back to int easily (though we must specify the old base explicitly):

>>> int(b,2)

35

>>> int(h,16)

35

To interconvert hex and bin, it's easiest to go through int:

>>> bin(35)

'0b100011'

>>> hex(35)

'0x23'

>>> hex(int(bin(35),2))

'0x23'

>>> bin(int(hex(35),16))

'0b100011'

The struct module

This module performs conversions between Python values and C structs represented as Python strings. This can be used in handling binary data stored in files or from network connections, among other sources. It uses Format Strings as compact descriptions of the layout of the C structs and the intended conversion to/from Python values.

unpack works on the particular string representation of a byte in which it is a single character:

>>> i = 35

>>> c = chr(i)

>>> c

'#'

>>> unpack('B',c)

(35,)

>>> b = bin(i)

>>> unpack('B',b)

Traceback (most recent call last):

  File "<stdin>", line 1, in <module>

struct.error: unpack requires a string argument of length 1

>>>

To construct a multi-byte value like a 4-byte unsigned int we can do this:

>>> h = '\x23'

>>> type(h)

<type 'str'>

>>> b_little = '\x23\x00\x00\x00'

>>> type(b_little)

<type 'str'>

With multi-byte objects (like ints) we get into endian-ness---'I' is the format for an unsigned int, and '>' means big-endian:

>>> unpack('I', b_little)

(35,)

>>> b_big = '\x00\x00\x00\x23'

>>> unpack('I', b_big)

(587202560,)

>>> unpack('>I', b_big)

(35,)

One can also use the binascii module, but this is more complicated. Let's write two bytes to a file and then read the binary data:

>>> fn = 'temp'

>>> FH = open(fn,'w')

>>> FH.write('##')

>>> FH.close()

>>> FH = open(fn,'rb')

>>> data = FH.read()

>>> FH.close()

This "data" is still a string:

>>> type(data)

<type 'str'>

>>> data

'##'

>>> len(data)

2

>>> b2a_hex(data[0])

'23'

>>> b2a_hex(data)

'2323'

>>> b2a_hex('#Az')

'23417a'

Or you can get fancy, but it seems unnecessary..

>>> import array

>>> L = array.array('B','#Az')

>>> type(L)

<type 'array.array'>

>>> L[0]

35

>>> L[1]

65

>>> L

array('B', [35, 65, 122])

>>> b2a_hex(L)

'23417a'

With all this in mind, we can redo the script from yesterday to read the base64-encoded key data in a way that is simply understandable. The output first:

> python read.py 

dlen 7

dlen 1

dlen 129



209,8,39,27 .. 169,97,145,127

0xd1,0x8,0x27,0x1b .. 0xa9,0x61,0x91,0x7f 



1467871546 .. 1722964351

We decode the data and break it into 3 parts exactly as before. Looking at the third part, we unpack it byte by byte. The result of unpacking is a list of ints. We convert that list directly into the large number n (by doing b * 256**i for each byte, moving along the list in reverse order), or we can convert the ints into hex values, assemble that list into one long hex string, and then call eval as before. It seems perfectly transparent now.

It even seems clear (in retrospect) why there is an additional null byte after the int that specifies the size of the third segment. It is probably to align the base64 encoding to begin with the first data byte.

The only remaining difficulty is that this approach fails for the private key. So that's still a mystery.

read.py

import base64

from struct import unpack



FH = open('data.txt','r')

data = FH.read()

FH.close()

b64_data = ''.join(data.strip().split())

data = base64.b64decode(b64_data)



L = list()

while data:

    dlen = unpack('>I',data[:4])[0]

    L.append(data[4:dlen+4])

    data = data[dlen+4:]

    print 'dlen', dlen

print



# let's look at n

bL = L[2]

bL = bL[1:]    # extra null first, why?

iL = [unpack('B',b)[0] for b in bL]

hL = [hex(i) for i in iL]



# look at the int and hex values

pL = [str(n) for n in iL]

print ','.join(pL[:4]), '..', ','.join(pL[-4:])

print ','.join(hL[:4]), '..',

print ','.join(hL[-4:]), '\n'



# first just do the computation ourselves

m = 256

iL.reverse()

n = iL[0]

for i in iL[1:]:

    n += i*m

    m *= 256

    

s = str(n)

print s[:10], '..', s[-10:]

    

# hex doesn't pad the output

# remove the '0x' and add the extra '0' if needed

for i in range(len(hL)):

    b = hL[i]

    if len(b) < 4:

        hL[i] = '0' + b[-1]

    else:

        hL[i] = b[-2:]  

h = '0x' + ''.join(hL)

nh = eval(h)

assert n == nh

Dissecting RSA keys in Python (2)

2011-08-30T18:09:00.003-04:00

The public key from last time was given (in id_rsa.pub) as:

AAAAB3NzaC1yc2EAAAABIwAAAIEA0QgnG+LBQosx

pbkJIU0eWV9Zf7L63fy3BXp6T9RQH84yISs8SMvF

V4J1NYJQtNpGopNFlvzsuFnnxaRaM8ureFIQVa8k

/dRzv3r7QyIuoGyNQLDk3K8Tse/55fjuIxJFEf5j

voJs3Z/jYdPHreg8wzNCB/uMkbCccYsVt6lhkX8=

(It has been broken into 40 character lines). If we look at the first 32 characters of this base64 encoding:

AAAAB3NzaC1yc2EA

AAABIwAAAIEA0Qgn

the decoded version of this in hex is:

0x0  0x0  0x0 0x7  0x73 0x73 0x68 0x2d 0x72 0x73 0x61 0x0

0x0  0x0  0x1 0x23 0x0  0x0  0x0  0x81 0x0  0xd1 0x8  0x27

The first four bytes are

0x0 0x0 0x0 0x7

interpreted as a 32-bit unsigned int that is equal to 7, and indicates that the next 7 bytes are grouped together to form this element of the data (which is 'rsa-ssh', as I mentioned before).

0x73 0x73 0x68 0x2d 0x72 0x73 0x61

   s    s    h    -    r    s    a

The next 4 bytes after that are:

0x0 0x0 0x0 0x1

indicating the following block is a single byte of data: 0x23. In ASCII-encoding this would be the character '#', but it's actually an unsigned int:

>>> int('0x00000023', 16)

35

This is the value of the exponent for the key pair. More on this in a moment. Finally we get to

0x0 0x0 0x0 0x81

>>> int('0x00000081', 16)

129

It turns out there are exactly 129 bytes left in the data. In the code from the first post (from here), we used the struct module to unpack these three segments of data (we labeled them "parts").

According to the docs, the struct module can be used to interpret strings as packed binary data.

The first argument to unpack is a format string. The first character of the format string can be used to indicate the byte order, size and alignment of the packed data

The '>' indicates big-endian data, data laid out in registers in the same order as it is in memory. While Intel x86 (and current Macs) are little-endian, network order is big-endian and (I suppose) the SSH standard respects this. The next format character 'I' indicates the data following are unsigned ints.

import struct



hL = ['0x0','0x0','0x0','0x81']

iL = [int(h,16) for h in hL]

cdata = ''.join([chr(n) for n in iL])

res = struct.unpack('>I', cdata)

>>> print res

(129,)

It returns a tuple and we want the first value.

In the code from yesterday, we simply copied the correct number of bytes in each segment (as "characters"---not ASCII).

In the last part of this example, we need to transform those bytes into numerical data. We'll use struct.unpack again, except this time specifying binary data. For the exponent we had a single byte '\x23'', corresponding to character '#' which is 35 as an unsigned int:

>>> struct.unpack('B','#')

(35,)

Our function f passed the first element of this tuple back to a big list comprehension.

e = eval('0x' + ''.join(['%02X' % f(x) for x in parts[1]]))

In the case of the exponent (e) it operated on only a single value, so I'll rewrite it to clarify what this does. Remembering that f(x) will return 35, we'll get rid of the list stuff and just do:

>>> s = '%02X' % 35

>>> s

'23'

>>> e = eval('0x' + s)

>>> 

>>> e

35

(I'm not too swift with format strings, but you can read about them here).

And it doesn't do much in this case! But suppose we had two of these bytes in a row. Then:

>>> s = ('%02X' % 35) * 2

>>> s

'2323'

>>> eval('0x' + s)

8995

>>> 35*256 + 35

8995

eval will turn our string of bytes into a multi-byte unsigned int. This is news to me. So, if we actually had the string of bytes in hex divided into two sets of bytes, those up to and including '0x81'

['0x0', '0x0', '0x0', '0x7', '0x73', '0x73', '0x68', '0x2d', '0x72', '0x73', '0x61', '0x0', '0x0', '0x0', '0x1', '0x23', '0x0', '0x0', '0x0', '0x81']

and those after:

L = ['0x0', '0xd1', '0x8', '0x27', '0x1b', '0xe2', 

      '0xc1', '0x42', '0x8b', '0x31', '0xa5', '0xb9', '0x9', '0x21', '0x4d', '0x1e', '0x59', 

      '0x5f', '0x59', '0x7f', '0xb2', '0xfa', '0xdd', '0xfc', '0xb7', '0x5', '0x7a', '0x7a', 

      '0x4f', '0xd4', '0x50', '0x1f', '0xce', '0x32', '0x21', '0x2b', '0x3c', '0x48', '0xcb', 

      '0xc5', '0x57', '0x82', '0x75', '0x35', '0x82', '0x50', '0xb4', '0xda', '0x46', '0xa2', 

      '0x93', '0x45', '0x96', '0xfc', '0xec', '0xb8', '0x59', '0xe7', '0xc5', '0xa4', '0x5a', 

      '0x33', '0xcb', '0xab', '0x78', '0x52', '0x10', '0x55', '0xaf', '0x24', '0xfd', '0xd4', 

      '0x73', '0xbf', '0x7a', '0xfb', '0x43', '0x22', '0x2e', '0xa0', '0x6c', '0x8d', '0x40', 

      '0xb0', '0xe4', '0xdc', '0xaf', '0x13', '0xb1', '0xef', '0xf9', '0xe5', '0xf8', '0xee', 

      '0x23', '0x12', '0x45', '0x11', '0xfe', '0x63', '0xbe', '0x82', '0x6c', '0xdd', '0x9f', 

      '0xe3', '0x61', '0xd3', '0xc7', '0xad', '0xe8', '0x3c', '0xc3', '0x33', '0x42', '0x7', 

      '0xfb', '0x8c', '0x91', '0xb0', '0x9c', '0x71', '0x8b', '0x15', '0xb7', '0xa9', '0x61', 

      '0x91', '0x7f', '0x0']

Note: The extra byte at the end here is caused by our having decoded by hand.

Paste it into the interpreter. Then, they have to be padded out...

for i in range(len(L)):

    s = L[i]

    c = s[-1]

    if len(s) < 4:  

         L[i] = '0' + c

    else:  

        L[i] = s[-2:]

Remove the first and last bytes and turn it into one long hex value:

>>> s = '0x' + ''.join(L[1:-1])

>>> eval(s)

146787154640181728129549360865206451974125178122992752499327844555082827713683060312609720092642289502088305950292885968241446393671268243539585900329449529436170858131743078225065287541399805133121852823856839448386268013284889428740476278604926908637093742329548018673369795976229837442944484597101722964351L

That matches what we got yesterday, and what the rsa module gave us for n from the private key.

This use of eval seems to be widely known, but I would never have appreciated it from just reading the docs.

Here is a quick script that shows all the values for each set of four characters in the base64 encoding, first 9 bytes of output first:

> python script2.py 

['A', 'A', 'A', 'A']

[0, 0, 0, 0]

['0b0', '0b0', '0b0', '0b0']

['000000', '000000', '000000', '000000']

['00000000', '00000000', '00000000']

0x0 0x0 0x0

0 0 0 



['B', '3', 'N', 'z']

[1, 55, 13, 51]

['0b1', '0b110111', '0b1101', '0b110011']

['000001', '110111', '001101', '110011']

['00000111', '01110011', '01110011']

0x7 0x73 0x73

7 115 115 



['a', 'C', '1', 'y']

[26, 2, 53, 50]

['0b11010', '0b10', '0b110101', '0b110010']

['011010', '000010', '110101', '110010']

['01101000', '00101101', '01110010']

0x68 0x2d 0x72

104 45 114

from string import *

v = True



cL = list(uppercase + lowercase + digits + '+/')

D = dict(zip(cL,range(len(cL))))

D['='] = 0



key = '''AAAAB3NzaC1yc2EAAAABIwAAAIEA0QgnG+LBQosxpbkJIU0eWV9Zf7L63fy3BXp6T9RQH84yISs8SMvFV4J1NYJQtNpGopNFlvzsuFnnxaRaM8ureFIQVa8k/dRzv3r7QyIuoGyNQLDk3K8Tse/55fjuIxJFEf5jvoJs3Z/jYdPHreg8wzNCB/uMkbCccYsVt6lhkX8='''



L = list(key)



while L:

    s = L[:4]

    L = L[4:]

    if v:  print s

    dL = [D[c] for c in s]

    if v:  print dL

    dL = [bin(c) for c in dL]

    if v:  print dL

    bL = list()

    for b in dL:

        n = 8 - len(b)

        b = n*'0' + b[2:]

        bL.append(b)

    if v:  print bL

    b = ''.join(bL)

    bL = [b[:8], b[8:16], b[16:]]

    if v:  print bL

    iL = [int('0b' + b,2) for b in bL]

    for i in iL:

        print hex(i),

    print

    for i in iL:

        print i,

    print '\n'

Dissecting RSA keys in Python

2011-08-29T20:25:00.001-04:00

I want to look at RSA keys a little bit more in this post. Let's generate a public/private pair (1024 bits) using the SSH utility keygen (no passphrase):

> ssh-keygen -b 1024 -t rsa 

Generating public/private rsa key pair.

Enter file in which to save the key (/Users/telliott_admin/.ssh/id_rsa): id_rsa

Enter passphrase (empty for no passphrase): 

Enter same passphrase again: 

Your identification has been saved in id_rsa.

Your public key has been saved in id_rsa.pub.

The key fingerprint is:

96:49:5d:40:b1:90:de:99:c5:05:43:f2:92:45:69:67 telliott_admin@c-98-236-78-154.hsd1.wv.comcast.net

The key's randomart image is:

+--[ RSA 1024]----+

|        .o=*B+.  |

|        .o B=.E  |

|       ...==.o   |

|       ..o+.     |

|        S        |

|       .         |

|                 |

|                 |

|                 |

+-----------------+

[ To be honest, at this point, I'm not actually sure this is the randomart that goes with the example in this post.. but I think it could be :) ]

Start with the "public" key. The data in the file looks like this:

id_rsa.pub

ssh-rsa AAAAB3NzaC1yc2EAAAABIwAAAIEA0QgnG+LBQosxpbkJIU0eWV9Zf7L63fy3BXp6T9RQH84yISs8SMvFV4J1NYJQtNpGopNFlvzsuFnnxaRaM8ureFIQVa8k/dRzv3r7QyIuoGyNQLDk3K8Tse/55fjuIxJFEf5jvoJs3Z/jYdPHreg8wzNCB/uMkbCccYsVt6lhkX8= my_identifier

Where my_identifier is

telliott_admin@c-myipaddress_with_dashes.hsd1.wv.comcast.net

The common representation of RSA keys (like we see in the file above) uses base64 format. I was confused about this at first, but found a great explanation here, and of course in wikipedia:

The base64, base32, and base16 encodings convert 8 bit bytes to values with 6, 5, or 4 bits of useful data per byte, allowing non-ASCII bytes to be encoded as ASCII characters for transmission over protocols that require plain ASCII, such as SMTP. The base values correspond to the length of the alphabet used in each encoding. There are also URL-safe variations of the original encodings that use slightly different results.

Consider the byte 1001-1110

We can encode this byte using hexadecimal representation as '9e'

>>> b = 0b10011110

>>> b

158

>>> hex(b)

'0x9e'

I didn't actually know you could do this, i.e. using 0b with no quotes (see SO).

The encoding of those 8 bits (one byte) as ASCII characters takes 2 bytes. The '9' and the 'e' take up that much space. An advantage of this representation of the byte is that it's easy to read, also the information can be transmitted to devices that are expecting an ASCII string, where the highest-value bit should be unset.

A disadvantage is that it doubles the size of the string that we're transmitting. As the quote indicates, the other encodings cram more bits of information into the byte to be transmitted. In base32 we could use (RFC4648) the characters 'A'..'Z' plus '2'..'7'.

In base64 we use uppercase, lowercase, digits plus '+' and '/', except that:

Using standard Base64 in URL requires encoding of '+', '/' and '=' characters into special percent-encoded hexadecimal sequences ('+' = '%2B', '/' = '%2F' and '=' = '%3D'), which makes the string unnecessarily longer.
For this reason, a modified Base64 for URL variant exists, where no padding '=' will be used, and the '+' and '/' characters of standard Base64 are respectively replaced by '-' and '_',

The wikipedia article has a great graphic showing an encoding in detail. We do conversions easily by using the Python base64 module.

>>> import base64

>>> e = base64.b64encode

>>> d = base64.b64decode

>>> e('Man')

'TWFu'

>>> d('TWFu')

'Man'

>>> s = 'Fourscore, and seven years ago'

>>> len(s)

30

>>> c = e(s)

>>> c

'Rm91cnNjb3JlLCBhbmQgc2V2ZW4geWVhcnMgYWdv'

>>> d(c)

'Fourscore, and seven years ago'

The encoding process looks at each sequence of 24 bits in the input (3 bytes) and encodes those same 24 bits spread over 4 bytes in the output. The last two characters [in the next example], the ==, are padding because the number of bits in the original string was not evenly divisible by 24 in this example.

>>> e('Ma')

'TWE='

We break this two byte (16 bit) value up into 6 + 6 + 6 = 18 bits---that makes three base64 values, and the third one is now E rather than F because the bit pattern is 000100, where the last two 00 are padding. Additional adding is added to bring us back to the byte boundary.

In Lincoln's famous phrase, there are 30 bytes, which is evenly divisible by 6, hence no padding is needed.

The second part of the problem was solved by a question I found on SO

Here is my version of the code:

import base64

import struct



data = open('id_rsa.pub').read().split(None)[1]

print data[:24]

keydata = base64.b64decode(data)

print keydata[:18]



parts = []

while keydata:

    dlen = struct.unpack('>I',keydata[:4])[0]

    data, keydata = keydata[4:dlen+4], keydata[4+dlen:]

    parts.append(data)



def f(x):

    return struct.unpack('B', x)[0]

    

e = eval('0x' + ''.join(['%02X' % f(x) for x in parts[1]]))

print 'e', e



n = eval('0x' + ''.join(['%02X' % f(x) for x in parts[2]]))

print 'n', n

I have to confess I haven't sorted this out completely yet. Nevertheless, it allows us to convert the RSA public key from ssh-keygen to an integer:

Output:

> python script.py

AAAAB3NzaC1yc2EAAAABIwAA

ssh-rsa#

e 35

n 146787154640181728129549360865206451974125178122992752499327844555082827713683060312609720092642289502088305950292885968241446393671268243539585900329449529436170858131743078225065287541399805133121852823856839448386268013284889428740476278604926908637093742329548018673369795976229837442944484597101722964351

I want to finish by showing that these values are actually correct, and tie the whole thing back to an earlier discussion of the mathematics of public key cryptography.

I found a module for Python that handles RSA cryptography.

One very nice thing about it is that it's pure Python. It might not be fast, but there's nothing to do to build it. (I used easy_install). It turns out that the rsa module can read the private RSA key that we got earlier (top of the post), but not the public key (which is why I learned everything that's come before this point):

> python

Python 2.7.1 (r271:86832, Jun 16 2011, 16:59:05) 

[GCC 4.2.1 (Based on Apple Inc. build 5658) (LLVM build 2335.15.00)] on darwin

Type "help", "copyright", "credits" or "license" for more information.

>>> import rsa

>>> with open('id_rsa') as f:  data = f.read()

... 

>>> data.split()[4][:50]

'MIICWgIBAAKBgQDRCCcb4sFCizGluQkhTR5ZX1l/svrd/LcFen'

>>>

>>> k = rsa.PrivateKey.load_pkcs1(data)

>>> k.n

146787154640181728129549360865206451974125178122992752499327844555082827713683060312609720092642289502088305950292885968241446393671268243539585900329449529436170858131743078225065287541399805133121852823856839448386268013284889428740476278604926908637093742329548018673369795976229837442944484597101722964351L

>>> k.e

35

>>> k.d

71296617968088267948638260991671705244575086516882194071102095926754516318074629294696149759283397758157177175856544613145845391211758861147798865874304045384903321200809536161346504709634159813307741235395231601606901597492591281205810328850837168587984349087962868939575631849929827946396907884195208774355L

>>> k.p

12214944519618901215070692893221534580187337690484280451921546594832793966376942655508370083135703324858458340368238021755353370324360760471666368429071147L

>>> k.q

12017013618393528112233953119917677248859562342149056793052943111068006176296374093622600767249370271694572460567324457033873003984957249202580919628769533L

>>> k.n == k.p * k.q

True

>>>

You can see that n equals p * q and that it matches what we eventually obtained from the public key file (id_rsa.pub) above. That's pretty cool.

Still to do: use the keys to do encryption and decryption and check that the math is right. Figure out how the hash/digest is done. And figure out how the passphrase is used to modify the private key file. We can even add a passphrase at this point to the old version:

ssh-keygen -f id_rsa -p

I used 'abcde' as a passphrase (weak!):

> cat id_rsa

-----BEGIN RSA PRIVATE KEY-----

Proc-Type: 4,ENCRYPTED

DEK-Info: AES-128-CBC,FB4D3BD371510C4B37552EF6E949CC48



E+KtbKVAwE4Dxj3vKYgfHpbgqF9Nk6O4vJ+C6wogXEHXr8vkM9SF5pz6lKNfGRbC

S93/lurvq+1+rCuldDNmEjdy36cbgrWJ1CIwB6fq5VyT3+v+AxDMqV9nQkvHjnxR

/ON2n+CjOSklnvmaBk3l8+IvCIKar2hGS/1qHqJRKJ3B1fgd5I8i9yH6hx8yl8Qm

AfVPgBJ4+91dCW2Y8odtXENgjstGtWZkYd0O+EqPuSuPHLu+J7M2KSQNe8aKOhGf

DXCEwR4aGW8EXfRNQpULbi5JruYebYlg3JILm2MaFwy+OW39q1VElvBc5wD01awW

BF7Y9x++bi9krX/G49NzkLdSVshkDMSP+eRZ4sncWvDCrkN2xlWqnNfQ06I+XCVW

FMUmyMFI87eNmjcZOgIi1S+fQtWcbl4o9Hx1z6xcMd7WxOnfIxe5zHIu8C24WONW

spcVxSsHAA57jYa8j292prtK6HDDa6UOAa32hBz1lvC8fMWPxwluTJ8eyeeK/7O6

hNI5G9gM4MsMEMau8O49fEm0RuiT/aB8IXZeEVQNbV/FdixyBBhY53DMnFQ5wau8

OJCpL51S5lj2o8eAJc61SHwX/Fc9uF+rmAFy08tthLuCF2wGdZWsE7IJ8NQzFJA+

Us2dTCjY57V/QoCC41rWh+YP+srcEK1Wp7DfvXcq98N1/+tECyCe/aDdPCIM/8QD

glKrx8Ah7Tsx6of3vCHiYG0q8Yu68NNJOiWrnGMEmP0Gq9yhLYqKrc8a3oIP0Lst

Ss+0p989df73GdM3nyqtLUwAzaCJQnDivLFT6gltEnM=

-----END RSA PRIVATE KEY-----

Duly quoted

2011-08-28T16:57:00.000-04:00

Recycled from an earlier post:

Programs must be written for people to read, and only incidentally for machines to execute.

- Abelson & Sussman, SICP

How do we convince people that in programming simplicity and clarity—in short: what mathematicians call "elegance"— are not a dispensable luxury, but a crucial matter that decides between success and failure?

- E. W. Dijkstra

Linux Server SSH

2011-08-28T16:07:00.004-04:00

To continue with my Ubuntu Linux server odyssey, I got SSH working this afternoon for access from my LAN. As always, heed the warnings, these posts are without a warranty of any kind, not least because I am a total amateur. They exist because it's easy to have Google to "index my brain" or rather, what was in my brain at a certain time.

The objective is to enable remote login for my Ubuntu machine running with a local ip address of 10.0.1.2. We begin by setting a stronger password for my account on Ubuntu under: System > Preferences > About Me (weirdly, it's not under Passwords and Encryption Keys ..).

On the client (OS X Snow Leopard Server), I already have some key pairs from a previous test

> ls ~/.ssh

authorized_keys	id_dsa		id_dsa.pub	known_hosts

but they are DSA keys, which the notes I'm going to follow deprecate. We'll make a short 1024 ~~byte~~ bit RSA pair for now:

ssh-keygen -b 1024 -t rsa



> cat id_rsa.pub 

ssh-rsa AAAAB..5Eoec= telliott_@c-98-___-__-154.hsd1.wv.comcast.net

Not sure what the pseudo ip address or "comcast" thing is about.

I also take this opportunity to set a passphrase for the public key (and of course I make a note of it). This will be useful (it says) because we can set up the server to require both a password to gain access to the public key, as well as the corresponding private key before login is allowed.

[UPDATE: I'm still working my way through this, but I suspect that the previous statement isn't correct. The passphrase is used to protect the value of the private key. So, in this setup OS X will decrypt the private key using the passphrase, and the stored value of the private key in id_rsa is encrypted. I'll try to figure all this out soon. ]

The easiest way to get my new public key onto the server is to use SSH (with password, before we disable it). The next question is, do we need to install an SSH server? It appears yes.

sudo apt-get install openssh-server

There is a pre-existing config file /etc/ssh/ssh_config but after the previous command there are more including sshd_config. I make sure to save a copy of this file before I modify it. Now try:

ssh te@10.0.1.2



ssh: connect to host 10.0.1.2 port 22: Connection refused

OK, so we need to modify /etc/ssh/sshd_config. Port22 is already uncommented. Now uncomment:

PermitRootLogin no

ChallengeResponseAuthentication yes

PasswordAuthentication yes   # we'll set it to no eventually

On the server:

sudo /etc/init.d/ssh restart

On the client:

> ssh te@10.0.1.2

The authenticity of host '10.0.1.2 (10.0.1.2)' can't be established.

RSA key fingerprint is d1: .. :83.

Are you sure you want to continue connecting (yes/no)?

Check the fingerprint on the the server:

ssh-keygen -l -f /etc/ssh/ssh_host_rsa_key.pub



2048 d1: .. :83 /etc/ssh/ssh_host_rsa_key.pub (RSA)



Are you sure you want to continue connecting (yes/no)? yes

Warning: Permanently added '10.0.1.2' (RSA) to the list of known hosts.

Connection closed by 10.0.1.2



> ssh te@10.0.1.2

te@10.0.1.2's password: 

Welcome to Ubuntu 11.04 (GNU/Linux 2.6.38-11-generic x86_64)

Now to change to not using our password. We need to copy our RSA public key to the server in some secure way. According to the sshd_config file, the "authorized keys file" is
%h/.ssh/authorized_keys. That's in my home directory. I do:

> scp ~/.ssh/id_rsa.pub te@10.0.1.2:~/.ssh/authorized_keys

te@10.0.1.2's password: 

id_rsa.pub                                    100%  260     0.3KB/s   00:00    



te@VB:~$ logout

Connection to 10.0.1.2 closed.



> ssh te@10.0.1.2

Identity added: /Users/telliott_admin/.ssh/id_rsa (/Users/telliott_admin/.ssh/id_rsa)

Welcome to Ubuntu 11.04 (GNU/Linux 2.6.38-11-generic x86_64)



* Documentation:  https://help.ubuntu.com/



Last login: Sun Aug 28 14:35:49 2011 from osxserver.local

te@VB:~$

Now, finally, be sure to turn off password authentication for the ssh server: PasswordAuthentication no. At first, I'm not being prompted for a password to retrieve my public key on the server.. just getting automatically logged in by:

ssh te@10.0.1.2

I forgot to do:

ChallengeResponsePasswords yes



sudo /etc/init.d/ssh restart

Restart and I get the challenge.. Although it is apparently cached sometimes. Looks like it works.

Linux Server

2011-08-28T11:05:00.006-04:00

A couple of weeks ago I described installation of Ubuntu Linux (11.04), running under Virtual Box on Mac OS X Lion. I got matplotlib (and scipy) installed and running without any difficulty using apt-get alone.

Next, I installed Apache Server 2.2 and configured everything to run Python scripts by using Firefox from Ubuntu on the VM like this:

localhost/cgi-bin/test.py

After struggling a bit, I figured out how to use Virtual Box in NAT mode to do "port forwarding"---by modifying

/etc/apache2/sites-available/default

so then I could launch the script from the host OS by doing

localhost:8080/cgi-bin/test.py

from Safari or by using curl in Terminal.

The next step was to try to access the server from a second machine on my home network. In order to do that, (as I understand it) I can either remain in NAT mode and try to connect via the host machine and port-forwarding, or use "bridging" mode to effectively put the Ubuntu server on the network.

In order to accomplish this, I need to know the ip address of the server as assigned by my Apple Time Capsule router. Now the problem is that this ip address may change with time (since we're getting it from the Time Capsule as a DHCP server). According to the folks I asked on serverfault, I need to switch the Ubuntu network IPv4 settings to manual and ask for a specific ip address and this is correct so far as it goes.

However, the address I ask for in Ubuntu

has to match the range of addresses that the router is set up to hand out. Furthermore, we probably want to have most ip addresses handled dynamically, and only use a static ip for the one machine. If we wanted them all static, we'd go here:

I found the answer here but the setting is actually visible in this screenshot:

Under Internet > DHCP there is a TableView called DHCP reservations. I just reserved the first available ip address. After re-booting Ubuntu and restarting the server, I can run the test script or do:

Now I have to figure out how to access the server from the outside world. I should probably remove the test script from cgi-bin (since it reveals information about paths on the server) and read up on firewall configuration, etc. before I do that. The first problem for this part (not yet solved) is that my ip address from my isp, obtained with

checkip.dyndns.org/

does not seem to be accessible. I get a "Safari can't connect to the server .." message.

Misunderstanding Euclid

2011-08-27T09:25:00.012-04:00

Category: things you thought you knew that are wrong
Today's entry is from mathoverflow:

Many students believe that 1 plus the product of the first n primes is always a prime number. They have misunderstood the contradiction in Euclid's proof that there are infinitely many primes. (By the way, 2 * 3 * 5 * 7 * 11 * 13 + 1 is not prime.)

Let's ask Python:

 
from math import sqrt

def f(N):

    N += 1

    print 'N', N

    n = int(sqrt(N))

    for i in range(2, n+1):

        if not N % i:

            print i, N/i



>>> N = 2*3*5*7*11*13

>>> f(N)

N 30031

59 509

>>> N *= 17

>>> f(N)

N 510511

19 26869

97 5263

277 1843

>>> N *= 19

>>> f(N)

N 9699691

347 27953

What the theorem actually says (wikipedia):

Take any finite list of prime numbers p₁, p₂, ..., p_n. It will be shown that at least one additional prime number not in this list exists. Let P be the product of all the prime numbers in the list: P = p₁p₂...p_n. Let q = P + 1.

Then, q is either prime or not:

If q is prime then there is at least one more prime than is listed.

If q is not prime then some prime factor p divides q. This factor p is not on our list: if it were, then it would divide P (since P is the product of every number on the list); but as we know, p divides P + 1 = q. Then p would have to divide the difference of the two numbers, which is (P + 1) − P or just 1. But no prime number divides 1 so there would be a contradiction, and therefore p cannot be on the list. This means at least one more prime number exists beyond those in the list.

Late nights on the web

2011-08-27T04:48:00.001-04:00

This is just a short post to document some interesting links. I recently posted about getting matplotlib installed with minimum fuss on OS X Lion, and this seems to have helped quite a few folks.

The person who showed me how to do this is named Rui Carmo, and I finally found the time to begin browsing through his extensive site; he has lots of good stuff about Macs and other technology (like here).

Perhaps the most interesting thing I found is what appears to be a very high quality company called Linode, offering Linux Virtual Servers in the cloud, their docs include an extensive library related to the topic. I'm thinking seriously about trying it (it only costs just a bit more than the NY Times).

And the final result (of browsing the web on a sleepless night) is to note the profusion of different sites connected to Stack Overflow, including serverfault and superuser. Very rewarding to browse through Q & A, and of course, when I have a question of my own.

One more thing: a pdf form-filler link that is not Adobe. One of the sites I use for work has Adobe forms on it that I can only read through Safari if I have Reader installed as the default pdf viewer. Maybe it will help with this, since I loathe Adobe.

Pretty code (6)

2011-08-24T18:38:00.023-04:00

This is absolutely the last post about html and code stuff. The whole adventure started here. It's a bit amusing to realize that if I had just understood what <br> tags were about (and <span> and <div>), this would never have happened.

I decided to write my own parser for Python code, so that I could post the results here. On the way, I learned something about "styles" and stylesheets. To organize things, I put this code for a "style" in the template for the blog at the end of <head>:

<style type='text/css'>

div.content

{

font-size:12pt;

font-family:arial;

color: #404040;

}



div.code

{ 

font-size:12pt;

font-family:courier;

margin: 20px;

border-style: solid;

padding: 5px;

border-width: 1px;

border-color:#ff00ff;

}



span.kw   { color:blue; }

span.str  { color:red; }

span.str3 { color:purple; }

span.py   { color:#00b0b0; }

span.cm   { color:green; }

</style>

The parser output looks great when run on itself. The first little bit is here:

But, alas, it was not to be. The Blogger editor is now (still) messing with my code, removing all the <br /> that I carefully put in, so that the informative whitespace is gone in the actual post (though not in the intial preview before posting). It looks hopeless at this point.

Files for the parser on Dropbox (here). The only known bugs are those visible in the comment lines above, where the strings should have been masked by being inside a comment.

Unfortunately, I do not know how to make this all work for Blogger.

What a waste! I am not happy with these guys. The only way I can see is to have my own server.

Pretty code (5)

2011-08-22T22:05:00.000-04:00

This is the first section of the same example from the last post (here, more explanation here), run through Pygments. It is self-contained and doesn't require loading any other resources.

The code to generate it was:

> pygmentize -f html -O full -o test.html pyparser.py

This is the same parser that is used for the official Python docs, which are generated in combination with Sphinx. I have only a hazy idea of how that's done.

It's easily customizable if you want.

import sys

from keyword import kwlist

from string import punctuation as punct

from utils import load_data

import html_tags as H

pywords = ['True','False','list','dict',

           'int','float','append','extend',

           'sys','argv','pop','open','write',

           'close','string','time']



try:  

    fn = sys.argv[1]

except IndexError:  

    fn = 'example.py'

data = list(load_data(fn))

Pretty code (4)

2011-08-22T15:54:00.000-04:00

Explanation of what this is about from last time here. I got triple-quoted strings working, and I'm really pleased with the results. More another time on other options for code-highlighting.


import sys
from keyword import kwlist
from string import punctuation as punct
from utils import load_data
import html_tags as H
pywords = ['True','False','list','dict',
           'int','float','append','extend',
           'sys','argv','pop','open','write',
           'close','string','time']

try:  
    fn = sys.argv[1]
except IndexError:  
    fn = 'example.py'
data = list(load_data(fn))

D = {'is_cm':False,
     'in_str_1':False,
     'in_str_2':False,
     'in_str_3':False }

L = list()

while data:
    c = data.pop(0)
    # comments first
    if c == '#':
        if not (D['in_str_1'] or D['in_str_2']):
            L.extend(list(H.cm_start))
            D['is_cm'] = True
    if c == "\n" and D['is_cm']:
        L.extend(list(H.cm_stop))
        D['is_cm'] = False
    L.append(c)
    
    a = '''triple-quoted-string
    with a continuation and two keywords'''
    
    # triple-quoted strings
    if c == "'" and len(data) > 1:
        if data[0] == "'" and data[1] == "'":
            if (not D['in_str_1'] and not D['in_str_2']):
                if not D['in_str_3']:
                    L.pop()
                    L.extend(list(H.str3_start))
                    L.extend(["'"] * 3)
                    data.pop(0)
                    data.pop(0)
                    D['in_str_3'] = True
                else:
                    data.pop(0)
                    data.pop(0)
                    L.extend(["'"] * 2)
                    L.extend(list(H.str3_stop))
                    D['in_str_3'] = False
                    continue
    
    # single-quoted strings
    if c == "'":
        if D['in_str_3']:
            continue        
        if (not D['in_str_1']):
            if not D['in_str_2']:
                # start a str_1
                L.pop()
                L.extend(list(H.str_start))
                L.append(c)
                D['in_str_1'] = True
            else:
                # already in str_2 or str_3
                pass
        else:
            if D['in_str_1']:
                # terminate str_1    
                L.extend(list(H.str_stop))
                D['in_str_1'] = False
            
    # double-quoted strings
    if c == '"':
        if not D['in_str_2'] and not D['in_str_3']:
            if not D['in_str_1']:
                # start a str_2
                L.pop()
                L.extend(list(H.str_start))
                L.append(c)
                D['in_str_2'] = True
            else:
                # already in str_1 or str_3
                pass
        else:
            if D['in_str_2']:
                # terminate str_2    
                L.extend(list(H.str_stop))
                D['in_str_2'] = False
s = ''.join(L)

# keywords last
pL = list()
D['is_str_3'] = False
for line in s.split('\n'):
    D['is_cm'] = False
    words = line.split()
    for w in words:
        # no kw highlighting in comments
        if w.startswith(H.cm_start):
            D['is_cm'] = True
        if w.startswith(H.str3_start):
            D['is_str_3'] = True
        if H.str3_stop in w:
            D['is_str_3'] = False
        if not D['is_cm'] and not D['is_str_3']:
            if w in kwlist:
                r = H.kw_start + w + H.kw_stop
                line = line.replace(w, r)
            for p in pywords:
                if p in w:
                    L = w.split(p)
                    if L[0] and not L[0][-1] in punct:
                        continue
                    if L[1] and not L[1][0] in punct:
                        continue
                    r = H.py_start + p + H.py_stop
                    line = line.replace(p, r)
    pL.append(line)
    
s = H.br.join(pL)
pL = [H.head, H.hr, s, H.hr, H.tail]
s = '\n'.join(pL)

fn = fn.split('.')[0] + '.html'
FH = open(fn,'w')
FH.write(s + '\n')
FH.close()

Here's a screenshot of html_tags.py:

Pretty code (3)

2011-08-22T06:32:00.013-04:00

I worked a bit more on the Python parser that I described last time (here). Actually, I started over from scratch. The first approach was based on breaking the code into words for processing, but I realized that a "stream" approach of one character at a time is better. The only limitation I'm aware of is that I still don't handle triple-quoted strings, but I think it's doable in the future.

As I mentioned before, I'm using <br /> tags in the code now rather than newlines.

The instructions for doing this said to modify the blog template, but I'm afraid of breaking the formatting for old posts. I put the following just below the <body> tag instead, as shown. It seems to work.

<body>



<style>

  cd { font-size:120%; }

  cm { color: green }

  kw { color: blue; }

  str { color: red; }

</style>

Then, to be totally meta about it, I ran the parser on itself. It is set up to grab HTML tags from a separate file (otherwise the meta application chokes on the tags). They are the ones shown above. Plus, there is a head and a tail to make an independent .html document. But for pasting into the blog, you don't need those guys. I invoked it like this:

python simple_parser2.py simple_parser2.py > example.html

And here is the code:



import sys
from keyword import iskeyword
from utils import load_data
import html_tags as H

try:  
    fn = sys.argv[1]
except IndexError:  
    fn = 'example.py'
data = list(load_data(fn))

D = {'is_cm':False,
     'is_str_1':False,'is_str_2':False }

L = list()
for c in data:
    # comments first
    if c == '#':
        if not (D['is_str_1'] or D['is_str_2']):
            L.extend(list(H.cm_start))
            D['is_cm'] = True
    if c == "\n" and D['is_cm']:
        L.extend(list(H.cm_stop))
        D['is_cm'] = False
    L.append(c)
    
    # single-quoted strings
    if c == "'":
        if not D['is_str_1']:
            if not D['is_str_2']:
                # start a str_1
                L.pop()
                L.extend(list(H.str_start))
                L.append(c)
                D['is_str_1'] = True
            else:
                # already in str_2
                pass
        else:
            # terminate str_1    
            L.extend(list(H.str_stop))
            D['is_str_1'] = False
    # double-quoted strings
    if c == '"':
        if not D['is_str_2']:
            if not D['is_str_1']:
                # start a str_2
                L.pop()
                L.extend(list(H.str_start))
                L.append(c)
                D['is_str_2'] = True
            else:
                # already in str_1
                pass
        else:
            # terminate str_2    
            L.extend(list(H.str_stop))
            D['is_str_2'] = False
s = ''.join(L)

# keywords last
pL = list()
for line in s.split('\n'):
    D['is_cm'] = False
    words = line.split()
    for w in words:
        # no kw highlighting in comments
        if w.startswith(H.cm_start):
            D['is_cm'] = True
        if not D['is_cm'] and iskeyword(w):
            r = H.kw_start + w + H.kw_stop
            line = line.replace(w, r)
    pL.append(line)
s = H.br.join(pL)

pL = [H.head, H.hr, s, H.hr, H.tail]
s = '\n'.join(pL)

try:  
    fn = sys.argv[2]
except IndexError:  
    fn = 'example.html'
FH = open(fn,'w')
FH.write(s + '\n')
FH.close()

Pretty code (2)

2011-08-21T19:56:00.005-04:00

So the question is, what are the good ways to format code for blog posts, when you don't control the machine that serves the pages (like on Blogger). The first---and motivating---post in this series is here.

I always prefer simple to complex. I even prefer simple html pages, although I'll ignore the obvious problem that pages are not simple any more when Blogger gets through with them.

So I liked the idea of controlling the color of text in a page with self-defined tags. For example, we can put this in the head of our html document:

It's a seductively simple idea. But.. it requires a way of parsing Python code to generate the appropriate tags. I thought if I enforced some rules that could be done relatively easily.

After a few hours, I have something that sort of works, but it has too many bugs to be useful yet.

The main simplifying rule is that multi-line strings are not allowed. Also, at this point I assume that indentation is a four-space tab.

You can see the results in the graphic (above). I put the script on Dropbox (simple_parser.py). Without a lot more debugging, this seems like it's probably not the best option, but it was an entertaining exercise for the afternoon.

Formatting problem solved

2011-08-21T08:56:00.001-04:00

[ UPDATE: So I discovered that the line breaks setting below seems to be retroactive, and changing it messes up the old posts. Another solution is to substitute newlines in the code with <br /> tags. That's what I've done here, and seems OK. This makes it difficult to actually edit code samples in the editor, but them's the breaks. :) ]

On the topic of how to format code samples for the blog, I've been doing some testing and found that the combination of


Settings > Formatting
section:  Convert line breaks
No

and


Settings > Basic
section:  Global > Select post editor
Updated editor

appears to have solved the problems that cropped up (my post here) with my standard approach using <table>. A sample:

def f(n):
print "Hello world!"

And I like it that the new editor's preview mode gives an accurate view of the product.

However, there are other solutions to consider, including how to get syntax highlighting without needing to load resources that are out of my control (my SO question is here), or whether to use a code hosting service and "embed" the code here. I'll have to get back to you on that.

Bioonductor: summary

2011-08-20T11:04:00.015-04:00

[ Update: Sorry for the dup. I tried using the new formatting on this one, but found a problem, which I hope is fixed now.]

I promise this is the last one. Here we load the sample data from ExpressionSet. I show (once again) how to obtain all the relevant stuff: expression data, sample names, feature names, and gene names.

One of the nice things about ExpressionSet objects is the very fancy indexing one can do with them. A great example is given on p. 166 of the Bioconductor book.

I made an stab at implementing range selectors here but it's almost trivial. Also, in this one, the Python code does not contain any informative whitespace, so I've used the <code> tag.
The script is first and then the output. If you can't figure something out, just ask.

script.py

import rpy2.robjects as robjects 

from rpy2.robjects.packages import importr  

import rpy2.rinterface as ri  

import bioc.biobase as biobase 



r = robjects.r 

g = robjects.globalenv 

#gd = importr('grDevices') 

dim = r['dim'] 



r('data("sample.ExpressionSet")') 

sample_ExpressionSet = r['sample.ExpressionSet'] 

eset = sample_ExpressionSet 

#print tuple(ri.globalenv) 

print eset 



print dim(eset) 

ed = eset.exprs 

print ed.rx(1,1) 



# no way to do es[x,y] from Python? 

# just wrap it 

r(''' 

do_sel <-function(obj,sel,bycol=FALSE) {

if (bycol) { return(obj[,sel]) } 

else { return(obj[sel,]) } 

} 

''') 



do_sel = r['do_sel'] 

# both work: 

v = robjects.IntVector([2,3,6]) 

v = robjects.BoolVector([True,False]*13) 

eset_sub = do_sel(eset,v,bycol=True) 

eset_sub = do_sel(eset_sub,v) 

print dim(eset_sub) 



# no explicit featureData, just from data rows 

rownames = r['rownames'] 

fd = rownames(ed) 

print fd[:2] 



# get the gene names 

importr('annaffy') 

r(''' 

get.genenames <- function(probeids) {

symbols <- aafSymbol(probeids, 'hgu95av2') 

getText(symbols) } 

''') 

get_genenames = r['get.genenames'] 

gn = get_genenames(fd) 



# slightly fancy, no dups, original order 

from collections import OrderedDict 

z = zip(gn,[None]*len(gn)) 

D = OrderedDict(z) 

L = D.keys() 

for i in range(4):  print L[i*5:i*5+5] 

print 



pd = eset.do_slot('phenoData') 

data = pd.do_slot('data') 

v = robjects.IntVector(range(1,4)) 

print do_sel(data,v)

> python script.py 

ExpressionSet (storageMode: lockedEnvironment) 

assayData: 500 features, 26 samples  

element names: exprs, se.exprs  

protocolData: none 

phenoData 

sampleNames: A B ... Z (26 total) 

varLabels: sex type score 

varMetadata: labelDescription 

featureData: none 

experimentData: use 'experimentData(object)' 

Annotation: hgu95av2  



Features  Samples  

500       26  



[1] 192.742 



Features  Samples  

3        3  



[1] "AFFX-MurIL2_at"  "AFFX-MurIL10_at" 





******* Deprecation warning *******: 



The package 'hgu95av2' is deprecated and will not be supported in the future.  



Instead we strongly reccomend that you should start using the 'hgu95av2.db' package.  



We hope you will enjoy these new packages.  If you still have questions, you can search  

(http://dir.gmane.org/gmane.science.biology.informatics.conductor) or ask the mailing list  

(http://bioconductor.org/docs/mailList.html).   





Loading required package: hgu95av2 

['', 'STAT1', 'GAPDH', 'ACTB', 'TFRC'] 

['C15orf31', 'ACTL6B', 'GLRA1', 'KCNB2', 'MGAT5'] 

['BMP3', 'RPL14', 'ZCWPW2', 'KITLG', 'GPR12'] 

['TRA@', 'ATP8A2', 'C1orf68', 'FAM20B', 'SLC34A1'] 



sex    type score 

A Female Control  0.75 

B   Male    Case  0.40 

C   Male Control  0.73 



>

Bioconductor: multiple plots

2011-08-19T19:26:00.009-04:00

Actually, I have a bit more to do with Bioconductor ExpressionSets, as long as I'm thinking about it. When I first played with the ALL example, I tried other pairs besides ('ALL1/AF4','BCR/ABL')---I'd have to look now, but they are probably in the original paper too. Anyway, in the first talk I gave about this stuff, I made plots of all the pairs. Some are really nice. So that's what I'd like to show now.

The difference is that then, I wrote the code in R, and I did not like the experience.

In Python, we do the standard imports. We factor out the R code into several functions, then load it from a text file in Python, and execute the code in R-space. In Python, we pick pairs to compare and call a plotting function in R.

The R code is a little hackish (I'm thinking of color.map), but this is a good start.

R code:


color.map <- function(g) {

  colors = c('red','green','blue',

    'cyan','magenta','maroon')

  if (g == 'ALL1/AF4') i=1

  if (g == 'BCR/ABL') i=2

  if (g == 'E2A/PBX1') i=3

  if (g == 'NEG') i=4

  if (g == 'NUP-98') i=5

  if (g == 'p15/p16') i=6

  colors[i]

}



library('limma')

library('annaffy')



get.genenames <- function(eset) {

  probeids <- featureNames(eset) 

  symbols <- aafSymbol(probeids, 'hgu95av2') 

  getText(symbols)

}



compare.factors <- function(eset,L) {

  eset.sub = eset[,eset$mol.biol %in% L]

  f <- factor(as.character(eset.sub$mol.biol))

  design <- model.matrix(~f)

  fit <- eBayes(lmFit(eset.sub,design))

  pv = fit$p.value[,2]

  sel <- p.adjust(pv) < 0.001

  eset.sel <- eset.sub[sel,]

}



plot.eset <-function(eset,...) {

  patient.colors <- unlist(lapply(

    eset$mol.biol,color.map))

  genenames = get.genenames(eset)

  heatmap(exprs(eset),

    col=topo.colors(100),

    ColSideColors=patient.colors,

    labRow=genenames)

}

The R code is in a file Rcode.txt.

Python code:


import rpy2.robjects as robjects

from rpy2.robjects.packages import importr 

import bioc.biobase as biobase



r = robjects.r

g = robjects.globalenv

gd = importr('grDevices')



importr('ALL')

r('data("ALL")')

ALL = g['ALL']

eset = biobase.ExpressionSet(ALL)

#--------------------------------------------



FH = open('Rcode.txt','r')

s = FH.read()

FH.close()

r(s)

compare_factors = r['compare.factors']

plot_eset = r['plot.eset']

#--------------------------------------------

targets = [ ['ALL1/AF4','NEG'],

            ['ALL1/AF4','BCR/ABL'],

            ['E2A/PBX1','BCR/ABL'] ]



def get_ofn(t):

    #print t, len(t)

    u,v = t

    s = u + '_' + v

    ofn = s.replace('/','-') + '.pdf'

    return ofn

    

def plot(eset_sub,t):

    ofn = get_ofn(t)

    gd.pdf(ofn)

    plot_eset(eset_sub,ofn)

    gd.dev_off()

               

for t in targets:

    eset_sub = compare_factors(eset,t)

    plot(eset_sub,t)

The first plot is at the top, and the second and third are below. It's pretty clear that there are significant commonalities according to mol.biol. BTW, I adjusted the target p-value to keep the number of hits down.

I also added code to substitute the gene names in the plots. That's pretty cool.

And I have to figure out what to do about code listings. <pre> looks good in preview, but we get the extra newlines; <code> doesn't give the right spacings, but it doesn't have the double newlines. I decided to go with <pre> today.

There are problems with the alternatives. That's for another time.

Learning to use RPy2 (5)

2011-08-19T15:13:00.003-04:00

Here is the result of my experiment with the rpy2-bioconductor-extensions-0.2.



import rpy2.robjects as robjects

from rpy2.robjects.packages import importr 

import bioc.biobase as biobase

r = robjects.r

g = robjects.globalenv



importr('ALL')

r('data("ALL")')

data = g['ALL']

eset = biobase.ExpressionSet(data)

print eset

Run as a script.



> python script.py

First (above) section output:



ExpressionSet (storageMode: lockedEnvironment)

assayData: 12625 features, 128 samples 

  element names: exprs 

protocolData: none

phenoData

  sampleNames: 01005 01010 ... LAL4 (128 total)

  varLabels: cod diagnosis ... date last seen (21 total)

  varMetadata: labelDescription

featureData: none

experimentData: use 'experimentData(object)'

  pubMedIds: 14684422 16243790 

Annotation: hgu95av2

Part 2:



expr_data = eset.exprs

print len(expr_data)

L = tuple(expr_data[:5])

print [round(t,3) for t in L]

output:



1616000

[7.597, 5.046, 3.9, 5.904, 5.925]

Part 3:



print list(eset.slotnames())

phenoData = eset.do_slot('phenoData')

print 'phenoData', id(phenoData)

print phenoData

output:



['experimentData', 'assayData', 'phenoData', 'featureData', 'annotation', 'protocolData', '.__classVersion__']

phenoData  0x10bbc35a8

An object of class "AnnotatedDataFrame"

  sampleNames: 01005 01010 ... LAL4 (128 total)

  varLabels: cod diagnosis ... date last seen (21 total)

  varMetadata: labelDescription

Part 4:



pd = eset.get_phenodata()

print 'pd', id(pd)

print pd

print list(pd.slotnames())

print

output:



pd 0x10bbc3488

An object of class "AnnotatedDataFrame"

  sampleNames: 01005 01010 ... LAL4 (128 total)

  varLabels: cod diagnosis ... date last seen (21 total)

  varMetadata: labelDescription



['varMetadata', 'data', 'dimLabels', '.__classVersion__']

Part 5:



data = pd.do_slot('data')

print 'data'

print type(data)

mb = pd.pdata.rx2('mol.biol')

print list(mb)[:10]

levels = r['levels']

print levels(mb)

print list(data.rx2('mol.biol'))[:10]

output:



data

<class 'rpy2.robjects.vectors.DataFrame'>

[2, 4, 2, 1, 4, 4, 4, 4, 4, 2]

[1] "ALL1/AF4" "BCR/ABL"  "E2A/PBX1" "NEG"      "NUP-98"   "p15/p16" 



[2, 4, 2, 1, 4, 4, 4, 4, 4, 2]

Part 6:



fd = eset.get_featuredata

print fd

# fd() gives AttributeError, no _featureData

# 

featureNames = robjects.r['featureNames']

fn = featureNames(eset)

print len(fn), fn[0]

output:



<bound method ExpressionSet.get_featuredata of >

12625 1000_at

In the last part, the "new" way doesn't seem to work yet. So I did as before.
Well, that's it for Bioconductor and Rpy2 for now. The presentation could be a lot better organized, but at least we figured out how to do most things.

Learning to use RPy2 (4)

2011-08-19T08:00:00.002-04:00

This is a brief post about how to access the components (attributes) of S3 and S4 classes from R, using RPy2. Just a few notes on what I've found, but YMMV, of course. There's a pretty detailed post about this kind of thing here. And of course, the R documentation is extensive (but more complex).

Note: the Bioconductor classes can be handled by using extensions to RPy2. That's for a later post. First some standard imports and a bit of setup:



>>> import rpy2.robjects as robjects 

>>> from rpy2.robjects.packages import importr 

>>> r = robjects.r 

>>> g = robjects.globalenv 

>>>

We define a toy S3 class starting from a vector, simply by giving it the "class" attribute, and then define a version of print that can work on it:



>>> r('''

... x <- rep(0:1, c(3,6))

... myS3_class <- x

... class(myS3_class) <- 'myS3_class'

... print.myS3_class <- function(x, ...) {

...   cat("This is my vector:\n")

...   cat(paste(x[1:5]), "...\n")

... }

... ''')

<SignatureTranslatedFunction - Python:0x1446648 / R:0xc5bf34>



>>> myS3_class = r['myS3_class']

>>> myS3_class.rclass

<rpy2.rinterface.SexpVector - Python:0x194550 / R:0x2c51048>

>>> myS3_class.r_repr()

'structure(c(0L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 1L), class = "myS3_class")'

>>> attributes = r['attributes']

>>> aL = attributes(myS3_class)

>>> type(aL)

<class 'rpy2.robjects.vectors.ListVector'>

>>> aL.rx2('class')

<StrVector - Python:0x23b00a8 / R:0x2c51048>

['myS3_class']

>>> len(aL)

1

>>> aL.rx2('structure')

rpy2.rinterface.NULL

To get at the actual numeric data:



>>> len(myS3_class)

9

>>> v = r('as.integer(myS3_class)')

>>> v

<IntVector - Python:0x23b0fa8 / R:0xb77c80>

[       0,        0,        0, ...,        1,        1,        1]

>>> structure = r['structure']

>>> structure(myS3_class)

<IntVector - Python:0x23b71e8 / R:0x2864a00>

[       0,        0,        0, ...,        1,        1,        1]

>>> tuple(myS3_class)

(0, 0, 0, 1, 1, 1, 1, 1, 1)

Now, look at symbols in the global "environment":



>>> import rpy2.rinterface as ri

>>> tuple(ri.globalenv)

('myS3_class', 'print.myS3_class', 'x')

Let's try a simple S4 class:



>>> r('''

... setClass(

...   Class="myS4_class",

...   representation=representation(

...     times = "numeric",

...     names = "character"

...   )

... )

... ''')

<StrVector - Python:0x23b0af8 / R:0xdc3068>

['myS4_class']

>>>

and instantiate one of them:



>>> r('''

... myS4_class = new(

...   Class="myS4_class", 

...   times=c(1,2,4),

...   names=c('a','b')

... )

... ''')

<RS4 - Python:0x23b0ad0 / R:0xa4d550>

>>>

Now, to access its slots:



>>> myS4_class = r['myS4_class']

>>> slotNames = r['slotNames']

>>> slotNames(myS4_class)

<StrVector - Python:0x23b0a58 / R:0xae3de8>

['times', 'names']

>>> myS4_class.do_slot('times')

<FloatVector - Python:0x23b0ad0 / R:0x2c755b8>

[1.000000, 2.000000, 4.000000]

>>>

Assign to one:



>>> v = robjects.StrVector('abcde')

>>> myS4_class.do_slot_assign('names', v)

>>> print(myS4_class)

An object of class "myS4_class"

Slot "times":

[1] 1 2 4



Slot "names":

[1] "a" "b" "c" "d" "e"



>>> v = myS4_class.do_slot('names')

>>> v



['a', 'b', 'c', 'd', 'e']

That's it for this example.

Learning to use RPy2 (3)

2011-08-18T22:31:00.006-04:00

So, I'd been thinking that I'd do the Bioconductor example in Python. The first group of commands in R is:



library('ALL') 

data('ALL') 

library('limma')

In Python we do:



>>> import rpy2.robjects as robjects

>>> from rpy2.robjects.packages import importr

>>> ALL = importr('ALL')

>>> limma = importr('limma')

and to get the data, we can (cheat a bit and) do



>>> robjects.r('data("ALL")')

<StrVector - Python:0x10eb15a28 / R:0x7f8371183f38>

['ALL']

>>> data = robjects.globalenv['ALL']

The second group of R commands is:



eset <- ALL[,ALL$mol.biol %in% c('BCR/ABL','ALL1/AF4')] 

f <- factor(as.character(eset$mol.biol)) 

design <- model.matrix(~f) 

fit <- eBayes(lmFit(eset,design)) 

sel <- p.adjust(fit$p.value[,2]) < 0.05 

esetSel <- eset[sel,]

I figured out how to take apart the data object above, and implement the first line. That was enough to convince me that it's not worth it.

I realized now that the use case is to get the underlying data into and back out of Python. That is, eventually we'll want to learn about how to construct an ExpressionSet object with data starting in Python data structures. But first:



>>> featureNames = robjects.r['featureNames']

>>> fn = featureNames(data)

>>> len(fn)

12625

>>> fn[:5]

<StrVector - Python:0x10eb19ef0 / R:0x7f8371eba940>

['1000_at', '1001_at', '1002_f..., '1003_s..., '1004_at']

>>> fn.rx(1)

<StrVector - Python:0x10eb19440 / R:0x7f8371f85208>

['1000_at']

>>> tuple(fn[:4])

('1000_at', '1001_at', '1002_f_at', '1003_s_at')

That's what we'll want.

The expression data:



>>> exprs = robjects.r['exprs'] 

>>> ed = exprs(data)

>>> ed

<Matrix - Python:0x10eb10710 / R:0x1137e6000>

[7.597323, 5.046194, 3.900466, ..., 3.095670, 3.342961, 3.842535]

>>> dim = robjects.r['dim']

>>> dim(ed)

<IntVector - Python:0x10eb19cb0 / R:0x7f8372023548>

[   12625,      128]

>>> t = tuple(ed)

>>> t[:2]

(7.597322981163869, 5.046194285620063)

>>> len(t)

1616000

>>> import numpy as np

>>> A = np.array(list(t))

>>> A.shape = tuple(dim(ed))

>>> A.shape

(12625, 128)

>>> A[:5,:3]

array([[ 7.59732298,  5.04619429,  3.90046642],

       [ 5.20211564,  5.81351983,  5.81160593],

       [ 6.14364259,  7.1470674 ,  6.77438247],

       [ 6.76609498,  4.61191719,  4.91830898],

       [ 7.12201406,  7.61231255,  6.82803053]])

The phenotypic data including mol.biol:



>>> phenoData = robjects.r['phenoData']

>>> pd = phenoData(data)

>>> mb = pd.rx2('mol.biol') 

Traceback (most recent call last):

  File "", line 1, in 

AttributeError: 'RS4' object has no attribute 'rx2'

UPDATE: this way



>>> slotNames = r['slotNames']

>>> phenoData = r['phenoData']

>>> pd = phenoData(ALL)

>>> slotNames(pd)

<StrVector - Python:0x1023f77a0 / R:0x103140528>

['varMetada..., 'data', 'dimLabels', '.__classV...]

>>> 

>>> df = pd.do_slot("data")

>>> df.rx2('sex')

<FactorVector - Python:0x1023e2170 / R:0x10268da40>

[       2,        2,        1, ...,        2,        2, NA_integer_]

>>> mb = df.rx2('mol.biol')

>>> levels(mb)

<StrVector - Python:0x1023e2ea8 / R:0x103742d88>

['ALL1..., 'BCR/..., 'E2A/..., 'NEG', 'NUP-..., 'p15/...]

>>>

pData works:



>>> pData = robjects.r['pData']

>>> pd = pData(data)

>>> mb = pd.rx2('mol.biol') 

>>> mb

<FactorVector - Python:0x1154c5cf8 / R:0x7f8374e72010>

[       2,        4,        2, ...,        4,        4,        4]

>>> len(mb)

128

>>> mb.levels

<StrVector - Python:0x1154c5f38 / R:0x7f8371a5c660>

['ALL1..., 'BCR/..., 'E2A/..., 'NEG', 'NUP-..., 'p15/...]

>>> L = list(mb)

>>> L[:5]

[2, 4, 2, 1, 4]

>

Notice how the factor changed into ints (see the Rpy2 docs).

Now, to select the entries we want:



>>> i = mb.levels.index('ALL1/AF4') + 1

>>> j = mb.levels.index('BCR/ABL') + 1

>>> L = [True if e in [i,j] else False for e in L]

>>> sum(L)

47





>>> from rpy2.robjects import BoolVector

>>> bv = BoolVector(L)

>>> robjects.r('f <- function(e,bv) { e[,bv] }')

<SignatureTranslatedFunction - Python:0x1154cf2d8 / R:0x7f8371fa4400>

>>> f = robjects.r['f']

>>> eset = f(data,bv)

>>> dim(eset)

<IntVector - Python:0x1154cf638 / R:0x7f8375052c68>

[   12625,       47]

>>> eset.rclass

<rpy2.rinterface.SexpVector - Python:0x10eb17ee8 / R:0x7f83750f5698>

And that's what we set out to do.

Constructing an R object from data in Python will have to wait.

Another look at the selected data:



>>> t = tuple(exprs(eset))

>>> A = np.array(list(t))

>>> A.shape = tuple(dim(eset))

>>> A.shape

(12625, 47)

>>> B = A[:100,]

Might as well plot something while we're here:



>>> fig = plt.figure()

>>> plt.hot()

>>> plt.pcolormesh(B)

<matplotlib.collections.QuadMesh object at 0x10d727610>

>>> plt.colorbar() 

<matplotlib.colorbar.Colorbar instance at 0x10d707560>

>>> plt.savefig('example.png')

That is what is at the top of the post.

Learning to use RPy2 (2)

2011-08-18T16:06:00.005-04:00

We're working on RPy2. Before I continue with Bioconductor, we should work through more of the examples in the docs.

I rashly asked a question on Stack Overflow this morning (after being stumped for several hours), then stumbled across the answer. The question is, how to load the data from ALL. In R we would just do:



data('ALL')

and then the objects would be available to us through their names. It turns out that in Python it is almost that simple. I had actually seen the answer, which is a kind of "shortcut" in which we just execute R code directly, and then grab a variable through the name assigned in R. (I used this device several times in the previous post).

(And I should note that there is apparently a more elegant solution in the bioconductor extensions to rpy2.



>>> import rpy2.robjects as robjects

>>> from rpy2.robjects.packages import importr

>>> base =  importr('base')

>>> ALL = importr('ALL')

>>> data = robjects.r('data(ALL)')

>>> data.rclass

<rpy2.rinterface.SexpVector - Python:0x258b190 / R:0xdf86c8>

>>> data = robjects.globalenv['ALL']

>>> data

<RS4 - Python:0x2591490 / R:0x29a2134>

>>> data.rclass

<rpy2.rinterface.SexpVector - Python:0x258b3b0 / R:0xdf85c8>

>>> featureNames = robjects.r('featureNames')

>>> featureNames(data)

<StrVector - Python:0x23e4f08 / R:0x304fc00>

['1000..., '1001..., '1002..., ..., 'AFFX..., 'AFFX..., 'AFFX...]

>>> exprs = robjects.r['exprs']

>>> e = exprs(data)

>>> e

<Matrix - Python:0x23b2da0 / R:0x84d8000>

[7.597323, 5.046194, 3.900466, ..., 3.095670, 3.342961, 3.842535]

>>>

Anyway, we want not to make any more silly mistakes. Another elementary thing I happened upon in the docs (in Vectors) is

Access to R-style extracting/subsetting is granted though the two delegators rx and rx2, representing the R functions [ and [[ respectively.

It seems to return a strange little beast:



>>> from rpy2.robjects import r

>>> m = r.matrix(range(4),nrow=2)

>>> value = m.rx(3)

>>> value

<ListVector - Python:0x1023db5a8 / R:0x102cf9c28>

[IntVector]

  <no name>: 

  <IntVector - Python:0x1023db368 / R:0x102cf99e8>

[       2]

>>> type(value )

<class 'rpy2.robjects.vectors.ListVector'>

>>> value[0]

<IntVector - Python:0x1023db518 / R:0x102cf99e8>

[       2]

>>> value[0][0]

2

We copy the code to make an R list, directly out of the docs (see ?list):



>>> r('''e1 <- new.env(); e1$a <- 10; e1$b <- 20; L = as.list(e1)''')

<ListVector - Python:0x1023da4d0 / R:0x103306140>

[FloatVector, FloatVector]

  b: 

  <FloatVector - Python:0x1023da680 / R:0x102d1bd58>

[20.000000]

  a: <class 'rpy2.robjects.vectors.FloatVector'>

  <FloatVector - Python:0x1023da3b0 / R:0x102d1bd88>

[10.000000]

>>> L = r['L']

>>> L.rx2('a')

<FloatVector - Python:0x1023da998 / R:0x102d1bd88>

[10.000000]

Here is a somewhat fancier example involving a linear model (from the Intro). Look in the docs for the R code. In Python:



>>> from rpy2.robjects import FloatVector

>>> from rpy2.robjects.packages import importr

>>> from rpy2.robjects import globalenv as g

>>> 

>>> stats = importr('stats')

>>> base = importr('base')

>>> 

>>> ctl = FloatVector([4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14])

>>> trt = FloatVector([4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69])

>>> weight = ctl + trt

>>> 

>>> group = base.gl(2, 10, 20, labels = ["Ctl","Trt"])

>>> g["weight"] = weight

>>> g["group"] = group

>>> lm_D9 = stats.lm("weight ~ group")

>>> print(stats.anova(lm_D9))

Analysis of Variance Table



Response: weight

          Df Sum Sq Mean Sq F value Pr(>F)

group      1 0.6882 0.68820  1.4191  0.249

Residuals 18 8.7292 0.48496

>>> 

>>> print(lm_D9.names)

 [1] "coefficients"  "residuals"     "effects"       "rank"         

 [5] "fitted.values" "assign"        "qr"            "df.residual"  

 [9] "contrasts"     "xlevels"       "call"          "terms"        

[13] "model"        



>>> print names(lm_D9)

Traceback (most recent call last):

  File "<stdin>", line 1, in 

NameError: name 'names' is not defined

>>>

The next example in this section fails with two errors:

ImportError: cannot import name NA_real.

If we just use zeroes, it fails with



in __setitem__

    globalenv_ri)

TypeError: All keywords must be strings (or None).

I'm just going to leave it for now.

PCA example



>>> from rpy2.robjects.packages import importr

>>> from rpy2.robjects import globalenv as g

>>> from rpy2.robjects import r

>>> 

>>> base     = importr('base')

>>> stats    = importr('stats')

>>> graphics = importr('graphics')

>>> 

>>> # each of these x.y could be replaced by `r`

... m = base.matrix(r.rnorm(100), ncol=5)

>>> pca = stats.princomp(m)

>>> graphics.plot(pca, main="Eigen values")

rpy2.rinterface.NULL

The graphic is at the top of the post. I didn't both to print this properly. See the first post on the docs here. But the idea is that although these functions come from the imported modules they are present in the global R namespace (if that's what they call it). In general, it's better to document your code by using the qualified name. (I know, I know, I used r here extensively! Do as I say, not as I do).

kwargs example



>>> from rpy2.robjects import r

>>> r('''f <- function(x='foo',y='blue') { 

...   if (x=='foo')  print(x) 

...   else print (y) }''')

<SignatureTranslatedFunction - Python:0x1023d76c8 / R:0x100ef6b58>

>>> f = r['f']

>>> kwargs = {'x':"bar", 'y':"red"}

>>> f(**kwargs)

[1] "red"

<StrVector - Python:0x1023db170 / R:0x102d0fc28>

['red']

I guess that's enough for one post.

Random Quotes, get your Random Quotes

2011-08-18T13:31:00.003-04:00

A nice collection of quotes here. Perhaps my favorite:

"A month in the laboratory can often save an hour in the library."

-- F. H. Westheimer

Advice to a neuroscientist.

2011-08-18T13:14:00.002-04:00

Money quote for folks at my place:

Third: learn how to do your own data analysis. Know statistics well. Know at least some basic programming/scripting in Python, R, Matlab, etc. This will be of immense value in helping you get your research done efficiently and correctly, without needing to rely on other people's code (and time and commitment). This will become more important as our field becomes more data driven.

--Bradley Voytek (here)
h/t Tyler Cowen

My prior on the probability that they'll notice? 0.01 :)

Learning to use the RPy2 module (1)

2011-08-18T08:00:00.003-04:00

Before we tackle turning the Bioconductor example into Python code, I need to review some basic usage of RPy2 as given in the docs. I posted a bit about this before, but my exploration wasn't systematic, and also the posts contain issues which I solved and described in later posts.

This exercises the first half of what's on the linked page. Python code followed output in bold.



import rpy2.robjects as robjects

from rpy2.robjects.packages import importr



r = robjects.r

g = robjects.globalenv

print type(r)



<class 'rpy2.robjects.R'>



# from the R base package

# any other name comes from .globalEnv

pi = r['pi']

print pi[0]

print type(pi)



3.14159265359

<class 'rpy2.robjects.vectors.FloatVector'>



# r is callable with code to be evaluated

piplus = r('piplus = pi + 1')

print type(piplus)

print piplus

print piplus[0]



<class 'rpy2.robjects.vectors.FloatVector'>

[1] 4.141593



4.14159265359

[1] 4.141593



# explicitly get from .globalEnv

piplus_from_g = g['piplus']

print piplus_from_g

print piplus == piplus_from_g

print piplus[0] == piplus_from_g[0]



False

True



# define an R function and call it

r('''

    f <- function(start=1,stop = 5) {

    n = 0

    for (j in start:stop) { n = n + j }

    print (n)

    }

''')



# it's available in .globalEnv

f = g['f']

f()





[1] 15



# but also from the running R process

f = r['f']

f(4)



[1] 9



# interpolating an R object into code

letters = robjects.r['letters']

s = letters[1:6].r_repr()

rcode = 'paste(%s, collapse="-")' %(s) 

res = robjects.r(rcode)

print(res)



[1] "b-c-d-e-f"



# more on calling R functions

rsum = r['sum']

print rsum(robjects.IntVector([1,2,3]))[0]



6



# with keyword

rsort = r['sort']

iv = robjects.IntVector([3,1,2])

res = rsort(iv, decreasing=True)

print res.r_repr()



c(3L, 2L, 1L)



# plotting

gd = importr('grDevices')

ofn = '/Users/telliott_admin/Desktop/plot.pdf'

gd.pdf(ofn)



x = robjects.IntVector(range(10)) 

y = r.rnorm(10)

r.layout(r.matrix(robjects.IntVector([1,2,3,2]), nrow=2, ncol=2)) 

r.plot(r.runif(10), y, xlab="runif", ylab="foo/bar", col="red")

gd.dev_off()

Finally, there are R's special operators like %in% and %*%. I didn't find this yet, so try wrapping it:



r(''' 

is_in <- function(value, container) { 

  value %in% container

  }

''')



is_in = g['is_in']

L = robjects.StrVector(list('abcde'))

print is_in('a',L)

print is_in('f',L)



[1] TRUE



[1] FALSE



r('''m_mult <- function(m1,m0) { m1 %*% m0 }''')

m_mult = g['m_mult']

r('''mtoi <- function(m) { as.integer(m) }''')

mtoi = g['mtoi']



# an old friend

m0 = r.matrix(robjects.IntVector([0,1,1,1]), nrow=2)

m = m0

for i in range(10):

    print i + 3, 

    m = m_mult(m,m0)

    print mtoi(m)



3 [1] 1 1 1 2



4 [1] 1 2 2 3



5 [1] 2 3 3 5



6 [1] 3 5 5 8



7 [1]  5  8  8 13



8 [1]  8 13 13 21



9 [1] 13 21 21 34



10 [1] 21 34 34 55



11 [1] 34 55 55 89



12 [1]  55  89  89 144

Explanation of the example code

2011-08-18T05:33:00.005-04:00

This is a brief post to explain the R (Bioconductor) code that we used to draw a heat map for gene expression data (here and here). Again, first we load the example data:



library('ALL') 

data('ALL') 

library('limma')

The second group of commands is:



eset <- ALL[,ALL$mol.biol %in% c('BCR/ABL','ALL1/AF4')] 

f <- factor(as.character(eset$mol.biol)) 

design <- model.matrix(~f) 

fit <- eBayes(lmFit(eset,design)) 

sel <- p.adjust(fit$p.value[,2]) < 0.05 

esetSel <- eset[sel,]

The first of these constructs a "selector." The values in ALL$mol.biol are factors; two of the possibilities are in the vector c('BCR/ABL','ALL1/AF4'). The result is a vector of booleans. To give a simpler example:



> vowels = c('a','e','i','o','u')

> letters[1:6]

[1] "a" "b" "c" "d" "e" "f"

> letters %in% vowels

 [1]  TRUE FALSE FALSE FALSE  TRUE FALSE FALSE

 [8] FALSE  TRUE FALSE FALSE FALSE FALSE FALSE

[15]  TRUE FALSE FALSE FALSE FALSE FALSE  TRUE

[22] FALSE FALSE FALSE FALSE FALSE

The vector of True/False values from ALL$mol.biol %in% c('BCR/ABL','ALL1/AF4' is used to select the appropriate columns from the data, but is acting on the complex ExpressionSet object, as explained last time. The result is reduced in size.



> dim(eset)

Features  Samples 

   12625       47

The second line is just a bit of coaxing. Because the mol.biol descriptions for the smaller eset were derived from a "factor" with 6 possibilities, we remove the other 4 before moving on to model construction:



> f = eset$mol.biol

> class(f)

[1] "factor"

> length(f)

[1] 47

> levels(f)

[1] "ALL1/AF4" "BCR/ABL"  "E2A/PBX1" "NEG"     

[5] "NUP-98"   "p15/p16" 

> f = factor(as.character(f))

> class(f)

[1] "factor"

> length(f)

[1] 47

> levels(f)

[1] "ALL1/AF4" "BCR/ABL"

The middle two lines is where the magic happens. I'm going to skip that for the moment. You can read about it in (Gentleman et al, PMID 15461798). The result is in the variable `fit`. One of its sub-objects is a list of p-values:



> class(fit$p.value)

[1] "matrix"

> fit$p.value[1:3,]

           (Intercept)   fBCR/ABL

1000_at   2.713661e-58 0.05437095

1001_at   3.001944e-43 0.24423133

1002_f_at 2.273600e-46 0.10091357

> colnames(fit$p.value)

[1] "(Intercept)" "fBCR/ABL"   

> pv = fit$p.value[,2]

> class(pv)

[1] "numeric"

> pv[1:4]

   1000_at    1001_at  1002_f_at  1003_s_at 

0.05437095 0.24423133 0.10091357 0.30303775

p.adjust is a correction for the large number of comparisons. You can read about it by doing ?p.adjust. I wasn't clear on which is the default method, but it's not Bonferroni.



> r1 = p.adjust(pv)

> r2 = p.adjust(pv,method='bonferroni')

> all(r1 == r2)

[1] FALSE

So we just construct a selector for those rows with p < 0.05, and again, whittle down the ExpressionSet to something manageable:



> sel <- p.adjust(pv) < 0.05 

> esetSel <- eset[sel,]

> dim(esetSel)

Features  Samples 

     165       47

In the third part, we draw the heatmap.



color.map <- function(s) { 

if (s=='ALL1/AF4') 'red' else 'blue'} 

patient.colors <- unlist(lapply(esetSel$mol.biol,color.map)) 

heatmap(exprs(esetSel), 

col=topo.colors(100), 

ColSideColors=patient.colors)

In the first two lines we construct a function called color.map, and then use it to construct a vector of patient.colors. This is for the colored bars at the top of the plot.



> color.map <- function(s) { 

+ if (s=='ALL1/AF4') 'red' else 'blue'} 

>

We use lapply to feed values to our function, but the result is a list. So we use unlist to extract what we want.



> values = esetSel$mol.biol

> values[1:4]

[1] BCR/ABL  BCR/ABL  ALL1/AF4 BCR/ABL 

6 Levels: ALL1/AF4 BCR/ABL E2A/PBX1 ... p15/p16

> colors = lapply(values,color.map)

> class(colors)

[1] "list"

> colors[1:2]

[[1]]

[1] "blue"



[[2]]

[1] "blue"



> patient.colors = unlist(colors)

> patient.colors[1:4]

[1] "blue" "blue" "red"  "blue"

> class(patient.colors)

[1] "character"

We draw the heatmap with the command:



heatmap(exprs(esetSel)

All the rest is eye candy. In particular, we use the topo.colors. Without that, we get the default.

With the topo.colors, we would get colors as in the graphic at the top of the post.



heatmap(exprs(esetSel),col=topo.colors(100))

The last argument to heatmap gives us the bars on the columns to identify the sample by translocation.

We've explained everything but the actual statistical model fitting. That's for another time.

Dissecting an ExpressionSet object

2011-08-17T20:36:00.003-04:00

Last time, we were in R using some Bioconductor code. If we do this:



library('ALL')

data('ALL')

then the variable ALL is available



> ALL

ExpressionSet (storageMode: lockedEnvironment)

assayData: 12625 features, 128 samples 

  element names: exprs 

protocolData: none

phenoData

  sampleNames: 01005 01010 ... LAL4 (128 total)

  varLabels: cod diagnosis ... date last seen (21 total)

  varMetadata: labelDescription

featureData: none

experimentData: use 'experimentData(object)'

  pubMedIds: 14684422 16243790 

Annotation: hgu95av2

(note: your output may vary slightly depending on the R and Bioconductor version.)

Obviously the object that ALL refers to is complicated. Its components are called "slots" and they are accessible with the syntax: object_name@slot_name. For example:



pData = ALL@phenoData

You can find out which slots are available by parsing the above output or just doing:



> slotNames(ALL)

[1] "assayData"         "phenoData"        

[3] "featureData"       "experimentData"   

[5] "annotation"        "protocolData"     

[7] ".__classVersion__"

In turn the slots are also complex objects. Technically, a phenoData object is:



> class(ALL@phenoData)

[1] "AnnotatedDataFrame"

attr(,"package")

[1] "Biobase"

More important, its sub-objects are available with the syntax slot_name$sub-object_name. But first we must discover which names are available:



> pData = ALL@phenoData

> varLabels(pData)

 [1] "cod"            "diagnosis"      "sex"           

 [4] "age"            "BT"             "remission"     

 [7] "CR"             "date.cr"        "t(4;11)"       

[10] "t(9;22)"        "cyto.normal"    "citog"         

[13] "mol.biol"       "fusion protein" "mdr"           

[16] "kinet"          "ccr"            "relapse"       

[19] "transplant"     "f.u"            "date last seen"

So, for example:



> pData$sex[1:5]

[1] M M F M M

Levels: F M



> pData$mol.biol[1:5]

[1] BCR/ABL  NEG      BCR/ABL  ALL1/AF4 NEG     

Levels: ALL1/AF4 BCR/ABL E2A/PBX1 NEG NUP-98 p15/p16



> levels(pData$mol.biol)

[1] "ALL1/AF4" "BCR/ABL"  "E2A/PBX1" "NEG"     

[5] "NUP-98"   "p15/p16"

At this level, the objects are reasonably simple.



> class(pData$mol.biol)

[1] "factor"

> class(pData$sex)

[1] "factor"

Don't forget the other method names:



> sampleNames(pData)[1:5]

[1] "01005" "01010" "03002" "04006" "04007"



> featureNames(ALL)[1:8]

[1] "1000_at"   "1001_at"   "1002_f_at" "1003_s_at"

[5] "1004_at"   "1005_at"   "1006_at"   "1007_s_at"



> head(varMetadata(pData))

                                                         labelDescription

cod                                                            Patient ID

diagnosis                                               Date of diagnosis

sex                                                 Gender of the patient

age                                           Age of the patient at entry

BT                             does the patient have B-cell or T-cell ALL

remission  Complete remission(CR), refractory(REF) or NA. Derived from CR

And, I should point out that (without actually knowing it at the time) I shadowed a real function in the namespace by assigning to pData above.

If we reboot R and do:



library('ALL')

data('ALL')

we can use a different method for access:



> p = pData(ALL)

> class(p)

[1] "data.frame"

> names(p)

 [1] "cod"            "diagnosis"      "sex"           

 [4] "age"            "BT"             "remission"     

 [7] "CR"             "date.cr"        "t(4;11)"       

[10] "t(9;22)"        "cyto.normal"    "citog"         

[13] "mol.biol"       "fusion protein" "mdr"           

[16] "kinet"          "ccr"            "relapse"       

[19] "transplant"     "f.u"            "date last seen"

> p$sex[1:5]

[1] M M F M M

Levels: F M

The very first slot holds the gene expression data.



> aData = ALL@assayData

> class(aData)

[1] "environment"

I'm not really sure what that class is. And, unfortunately, it is not easy to find which sub-object names are available for this slot:



> names(aData)

NULL

> elements(aData)

Error: could not find function "elements"

> elementNames(aData)

Error: could not find function "elementNames"

except we can probably guess that exprs is one (based on the output of the call `ALL` above:



> e = aData$exprs

> class(e)

[1] "matrix"

> e[1:3,1:3]

             01005    01010    03002

1000_at   7.597323 7.479445 7.567593

1001_at   5.046194 4.932537 4.799294

1002_f_at 3.900466 4.208155 3.886169

> dim(e)

[1] 12625   128

> rownames(e)[1:8]

[1] "1000_at"   "1001_at"   "1002_f_at" "1003_s_at"

[5] "1004_at"   "1005_at"   "1006_at"   "1007_s_at"

> colnames(e)[1:3]

[1] "01005" "01010" "03002"

-------------------------------------------------------
So far, so good.

Now, obviously we can select rows and columns of a matrix like ALL@assayData$exprs, which we've assigned to the variable e:



> e[1:2,1:3]

           01005    01010    03002

1000_at 7.597323 7.479445 7.567593

1001_at 5.046194 4.932537 4.799294

What the Bioconductor guys have done is make this assignment valid for the ExpressionSet object itself:



> ALL[1:2,1:3]

ExpressionSet (storageMode: lockedEnvironment)

assayData: 2 features, 3 samples 

  element names: exprs 

protocolData: none

phenoData

  sampleNames: 01005 01010 03002

  varLabels: cod diagnosis ... date last seen (21 total)

  varMetadata: labelDescription

featureData: none

experimentData: use 'experimentData(object)'

  pubMedIds: 14684422 16243790 

Annotation: hgu95av2

That's either very cool or very confusing depending on where you sit.

Finally, the annotation. This requires another package:



ann = ALL@annotation

> ann

[1] "hgu95av2"



library('annaffy')

probeids <- featureNames(ALL)

symbols <- aafSymbol(probeids, 'hgu95av2')



> symbols <- aafSymbol(probeids, 'hgu95av2')

[1] "You are missing hgu95av2 looking to see if it is available."

Warning: unable to access index for repository http://brainarray.mbni.med.umich.edu/bioc/bin/macosx/leopard/contrib/2.11

Package hgu95av2 is available for download, would you like to install? [y/n] 

y

Warning: unable to access index for repository http://brainarray.mbni.med.umich.edu/bioc/bin/macosx/leopard/contrib/2.11

trying URL 'http://bioconductor.org/packages/2.6/data/annotation/bin/macosx/leopard/contrib/2.11/hgu95av2_2.2.0.tgz'

Content type 'application/x-gzip' length 9340297 bytes (8.9 Mb)

opened URL

==================================================

downloaded 8.9 Mb





The downloaded packages are in

	/var/folders/3g/3gHxNStTGZaWkGngOL39-++++TI/-Tmp-//RtmptUUJr1/downloaded_packages

Loading required package: hgu95av2



******* Deprecation warning *******:



The package 'hgu95av2' is deprecated and will not be supported in the future. 



Instead we strongly reccomend that you should start using the 'hgu95av2.db' package. 



We hope you will enjoy these new packages.  If you still have questions, you can search 

(http://dir.gmane.org/gmane.science.biology.informatics.conductor) or ask the mailing list 

(http://bioconductor.org/docs/mailList.html).

To repeat:



library('annaffy')

probeids <- featureNames(ALL)

symbols <- aafSymbol(probeids, 'hgu95av2')

> class(symbols)

[1] "aafList"

> genenames = getText(symbols)

> genenames[1:3]

[1] "MAPK3"   "TIE1"    "CYP2C19"

> probeids[1:3]

[1] "1000_at"   "1001_at"   "1002_f_at"

We'll write this all down on a couple post-it notes and put it on the monitor for reference.

A first heat map using Bioconductor

2011-08-17T17:15:00.005-04:00

I don't recall what the exact Google search target was that turned up a post about making a heatmap of gene expression data, but it really changed my life. It lead me to these two pages describing how to do it in R (here) and Python (here), and that motivated me to get into R, as much as I've been able to in spite of all I don't like about it.

I thought I'd work on this a little bit again; it's been 3 years or so since I gave a couple of talks about this, so I'm rusty, but it'll be fun. I would recommend that you go through the original posts by Peter Cock first. I'll try to add something of my own.

Of course, it's all based on the Bioconductor project in R, and there is a very nice write-up in the original paper (Gentleman et al, PMID 15461798). I even bought the book (our library sucks), and it's pretty good too, especially if you are serious about the field.

To begin with, to motivate things, we will recreate the graphic. We need a couple of R packages from Bioconductor to do this. You could get them from the Package Installer, but I was having trouble and so I tried this:



source("http://bioconductor.org/biocLite.R")

biocLite('Biobase')

biocLite('ALL')

biocLite('limma')

biocLite('annaffy')

Note: all these lines would actually look like this in R:



> source("http://bioconductor.org/biocLite.R")

For most of the commands here, I've left out the prompt so you can more easily copy and paste them. And I'm not showing the info that's printed as the packages are downloaded and installed.

We're going to just quickly do the plots, then later I'll explain more about it. First we need to load two modules and some data:



library('ALL')

data('ALL')

library('limma')

The first command will output:

Loading required package: Biobase



Welcome to Bioconductor



  Vignettes contain introductory material. To view, type

  'openVignette()'. To cite Bioconductor, see

  'citation("Biobase")' and for packages 'citation(pkgname)'.

Next we do the analysis (I know this won't make much sense yet, but we're just making sure we get the graphic without an error):



eset <- ALL[,ALL$mol.biol %in% c('BCR/ABL','ALL1/AF4')]

f <- factor(as.character(eset$mol.biol))

design <- model.matrix(~f)

fit <- eBayes(lmFit(eset,design))

sel <- p.adjust(fit$p.value[,2]) < 0.05

esetSel <- eset[sel,]

And in the third part we plot the data:



color.map <- function(s) {

  if (s=='ALL1/AF4') 'red' else 'blue'}

patient.colors <- unlist(lapply(esetSel$mol.biol,color.map))

heatmap(exprs(esetSel),

  col=topo.colors(100),

  ColSideColors=patient.colors)

The graphic is at the top of the post.

An overview of what we've done. We filtered for all samples with one of two chromosomal translocations. In the plots, the columns from the 'ALL1/AF4' samples are red, those with 'BCR/ABL' are blue. We now have 47 samples.



> dim(eset)

Features  Samples 

   12625       47

Each sample has 12625 "features"---gene expression measurements. We've done some fancy statistics and then filtered for those with p < 0.05 (correcting for the large number of measurements). In the end, we have 165 genes of interest to plot.



> dim(esetSel)

Features  Samples 

     165       47

By default, heatmap does clustering and draws dendrograms, but you can turn this off. Let's turn off the column clustering only, and since we haven't substituted the real gene names, let's turn that off too:



heatmap(exprs(esetSel),

  col=topo.colors(100),

  Colv=NA,

  ColSideColors=patient.colors,

  labRow=NA)

More later on this. Our goal is to understand more about "ExpressionSet" objects and the kinds of manipulation that are possible, the second is to see how much of this we can do from Python (using rpy2), so that the code is transparent and makes intuitive sense.

Zenoss

2011-08-17T04:25:00.004-04:00

I found this question on SF about network monitoring for OS X. One of answers described an open source solution called Zenoss.

From the looks of things it's very powerful. Here is a slightly old post about using it on OS X.

I downloaded and installed Zenoss. Unfortunately, it made my system feel unresponsive, even when the App was not running. (There are a bunch of process that run at startup, and eat up CPU).

Luckily, there is an uninstaller. I couldn't figure out how to get there in the Finder, but I just did:

open -a /usr/local/zenoss/uninstall.app/

It's probably no surprise that the uninstaller took more than 5 min to run, but it seems at least to have done its job and now everything is nice and snappy again.

High standards

2011-08-16T08:00:00.000-04:00

The nice thing about standards is that there are so many to choose from. And if you really don't like all the standards you just have to wait another year until the one arises you are looking for.

-- A. Tanenbaum, "Introduction to Computer Networks", quoted here

Trying Ubuntu Linux (6)

2011-08-15T08:30:00.003-04:00

Continuing with my Ubuntu Linux (11.0.4) server running in a VirtualBox VM under Mac OS X Lion, and making its capabilities available to other machines on my LAN. In two previous posts I described my search for a solution (here), and the actual solution (here), which used port forwarding with VirtualBox networking in NAT mode.

There is at least one other possible solution, which is to do the networking in bridged mode. To quote the docs:

With bridged networking, VirtualBox uses a device driver on your host system that filters data from your physical network adapter. This driver is therefore called a "net filter" driver. This allows VirtualBox to intercept data from the physical network and inject data into it, effectively creating a new network interface in software. When a guest is using such a new software interface, it looks to the host system as though the guest were physically connected to the interface using a network cable: the host can send data to the guest through that interface and receive data from it. This means that you can set up routing or bridging between the guest and the rest of your network.

To begin with I did something (that may not be necessary), I turned off the NAT rule from before (in OS X):

VBoxManage modifyvm Ubuntu --natpf1 delete "server"

Then I set networking to bridged in the VirtualBox Settings:

And restarted Apache (in Ubuntu):

sudo /etc/init.d/apache2 restart

ifconfig gives me the IP address that we have obtained from the DHCP server in Airport Extreme:

If I use the AirPort Utility I can see the corresponding MAC and IP addresses.

Now if I go to another machine on the network and point Safari at that IP, it works:

[ UPDATE: In the screenshot we specified port 8080, and that works because of what we did last time. But you can also do the ipadress only and it works too. ]

Pretty code

2011-08-14T15:07:00.013-04:00

[ UPDATE2: In an earlier version of this post, I missed one HTML tag in the code below, as a result the example wouldn't run as described when copy/pasted. Sorry. ]
I've been looking into solutions for posting code on the blog. Since the beginning, I've used <table> tags like this:

 <div style="overflow-x: scroll "> 
<table bgcolor="#ffffb0" border="0" width="100%" padding="4"> 
<tbody><tr><td><pre style=" hidden;font-family:monaco;"> 
my code here 
</pre></table></div>

With results looking like a bit like this, but larger :)

I borrowed the technique from Peter Norvig. He's moved on, but I hadn't.

In the last few days Blogger changed something that inserts newlines into such tables. I found another guy who noticed, but no solution. (Perhaps if I used the "new" editor?).

So, I googled around for syntax highlighting. Here is a listing of a bunch of methods. The one I decided to try is google-code-prettify. Like all the methods, it uses javascript and some css to do its work. Instructions are in the readme. The footprint in the html is pretty small, and I like the minimalist style.

To test it, I downloaded the distribution (see the link), and then put a file with the following html in the same directory as the source:

 
<?xml version="1.0" encoding="UTF-8" ?> 
<!DOCTYPE html> 
<html ... > 
. 
<head ... > 
. 
<link href="prettify.css" type="text/css" rel="stylesheet" /> 
<script type="text/javascript" src="prettify.js"></script> 

</head> 
<body onload="prettyPrint()"> 

<pre class="prettyprint" lang-python> 
def f(n, s='abc'): 
    print 'Hello world!' 
    for i in range(n): 
        # this is a loop 
        print s 
</pre> 

<pre class="prettyprint"> 
#include <math.h> 
#include "Python.h" 

int process (const char *c, double f) { 
    int i = 0; 
    if (c[0] == 'c') { 
        i = floor(f); 
    } 
    return i; 
} 
</pre> 


</body> 
</html>

and then drag-n-dropped it on Safari:

One small problem is this note:

Include the script and stylesheets in your document (you will need to make sure the css and js file are on your server, and adjust the paths in the script and link tag)

Server, what server?

Seriously, I'd be happy to do it, but I don't want to pay Comcast $$ for the privilege.

However, there is a very nice answer on Stack Overflow (here) that shows how to embed the relatively simple css instructions into the html document. Then you just load the code from

 
http://google-code-prettify.googlecode.com/svn/trunk/src/prettify.js 
http://google-code-prettify.googlecode.com/svn/trunk/src/lang-css.js

I tested it in a simple html document and it worked. Now I'll have to try to incorporate all this into a Blogger template.

[ UPDATE: I see it choked on the <math.h>. I'll have to remember that. ]

Trying Ubuntu Linux (5)

2011-08-14T10:43:00.005-04:00

The project I had yesterday (here) was to set up my server (Apache) in Linux under VirtualBox so that it is visible from the host (OS X Lion). I solved that eventually by sticking with NAT mode and issuing this command in Terminal (with Ubuntu powered down):

VBoxManage modifyvm Ubuntu --natpf1 "server,tcp,,8080,,8080"

This is exactly what the docs say to do. The extra twist was to modify /etc/apache2/ports.conf to add:

NameVirtualHost *:8080

Listen 8080

and to modify /etc/apache2/sites-available/default. I basically duplicated the whole <VirtualHost entry and changed the second part to read:

<VirtualHost *:8080>
..

And now it works. From the host in Safari or using curl in Terminal:

http://127.0.0.1:8080/cgi-bin/test.py



SERVER_SOFTWARE Apache/2.2.17 (Ubuntu)

SCRIPT_NAME /cgi-bin/test.py

SERVER_SIGNATURE

Apache/2.2.17 (Ubuntu) Server at localhost Port 8080

REQUEST_METHOD GET

SERVER_PROTOCOL HTTP/1.1

QUERY_STRING

..

HTTP_USER_AGENT Mozilla/5.0 (Macintosh; Intel Mac OS X 10_7) AppleWebKit/534.48.3 (KHTML, like Gecko) Version/5.1 Safari/534.48.3

HTTP_CONNECTION keep-alive

SERVER_NAME localhost

..

SERVER_PORT 8080

..

SERVER_ADMIN webmaster@localhost

HTTP_HOST localhost:8080

REQUEST_URI /cgi-bin/test.py

HTTP_ACCEPT text/html,application/xhtml+xml,application/xml;q=0.9,*/*;q=0.8

GATEWAY_INTERFACE CGI/1.1

..

HTTP_ACCEPT_LANGUAGE en-us

HTTP_ACCEPT_ENCODING gzip, deflate

The other thing I've done today is to install EMBOSS and PyCogent.

sudo apt-get install emboss

How easy is that! The rebase files are at

http://rebase.neb.com/rebase/rebase.files.html

We want # 5 and # 31 as described. chmod for the EMBOSS directories so we can write there:

sudo chmod 777 -R /usr/share/EMBOSS

rebaseextract

tell it where the files are, and it'll do its thing.

head -5 /usr/share/EMBOSS/data/REBASE/embossre.enz > cat



# REBASE enzyme patterns for EMBOSS

#

# Format:

# namepatternlenncutsbluntc1c2c3c4

#

put a sequence file DA19.txt on the Desktop

remap -sequence DA19.txt



DA19





                                                SgeI

                                                | PcsI

                                     SgeI       | SgeI 

                   Cac8I             | NspI     | | SgeI     HindIII

                   | PcsI Cac8I      | SgeI     | | | SgeI   |   Cac8I

                   \ \    \          \ \        \ \ \ \      \   \

          GACGAACGCTGGCGGCGTGCTTAACACATGCAAGTCGAACGGAGCGAATGAAAGCTTGCT



..

[ UPDATE: Something wrong with the formatting here. The HindIII site, AAGCTT, is at the far right. ]

Get Cython:

sudo apt-get install cython

Get Pycogent from http://sourceforge.net/projects/pycogent/:

tar -xzvf ~/Downloads/PyCogent-1.5.1.tgz

python setup.py build

sudo python setup.py install

I tested Pocogent by downloading sequences and making a phylogenetic tree as described here (I had to install a MUSCLE binary first). I pasted the source into four scripts.

Not exactly as before, but it looks reasonable.

Trying Ubuntu Linux (4)

2011-08-13T08:44:00.003-04:00

As discussed in previous posts (here, here, here), I have Ubuntu Linux 11.0.4 running under VirtualBox 4.1.0 on my Powerbook with the host running OS X Lion 10.7 (in a second partition, actually). I set up and tested Apache Server 2.2 on the guest. Now I'd like to access the guest server from the host machine, and later the other machines on my Wi-Fi network. I haven't yet succeeded, so yesterday I asked a question on StackOverflow, which got moved by someone to a related site, "serverfault" (here).

I got a couple of good suggestions, but didn't yet solve my problem. However, there are a bunch of "related questions" in the sidebar that look promising. Time to be (even more) systematic.

The VirtualBox manual (here) has a section: Configuring port forwarding with NAT.

So it's pretty clear that we could stay with NAT. In that case we need to do "port forwarding," although other people have suggested using "bridged mode." In the NAT method, we simply instruct the VM (as we would a router) that packets arriving on certain ports should be "forwarded to the guest, on the same or a different port." The manual says to do this from the command line in the host (the first line is their template, the second my implementation):



VBoxManage modifyvm "VM name" --natpf1 "guestssh,tcp,,2222,,22"

VBoxManage modifyvm Ubuntu --natpf1 "apache,tcp,,8888,,80"

There are six fields in the last argument: a name, which is "purely descriptive", the protocol for forwarding, and then two sets of an ip address and port. So, this example "forwards all TCP traffic arriving on the localhost interface (127.0.0.1) via port 2222 to port 22 in the guest":



VBoxManage modifyvm "VM name" --natpf1 "guestssh,tcp,127.0.0.1,2222,,22"

From the docs

Forwarding host ports < 1024 impossible:
On Unix-based hosts (e.g. Linux, Solaris, Mac OS X) it is not possible to bind to ports below 1024 from applications that are not run by root. As a result, if you try to configure such a port forwarding, the VM will refuse to start.

So that's why we used 8888. Under Port Forwarding we forward (in the second example) from the host port 8888 to the guest port 80.

I used the GUI to set this.

I believe we should not have to do the command line version. Just to be safe, I remove the forwarding rule I set when playing around, by doing this from the command line:



VBoxManage modifyvm Ubuntu --natpf1 delete "apache"

And then repeated setting up forwarding in the GUI (Settings > Network). To be sure that everybody has got the word, we:

Quit Ubuntu
Restart VirtualBox
Restart Ubuntu
Restart Apache with:

sudo /etc/init.d/apache2 restart

But it doesn't work. From the guest, I point Firefox at localhost or 127.0.0.1 and I can see the index page or run my scripts. Note that this works even though the Apache restart gave this message:

apache2: Could not reliably determine the server's fully qualified domain name, using 127.0.1.1 for ServerName

But from Safari on the host with the same ip address I get "Safari can't connect to the server" and if I add the port (127.0.0.1:8888) it just hangs. 127.0.1.1 also doesn't work.

Now, what could be wrong? From serverfault (here):

Iain says to try bridged mode.

I tested that early on but it didn't make any difference. Also, the docs clearly say that NAT should work. I'll have to test some more, but first..

Eric Fortis says:

"then from the host PC browser access the IP of the Ubuntu virtual machine, instead of 127.0.0.1. ifconfig will show you your IP."



te@VB:~$ ifconfig

eth0      Link encap:Ethernet  HWaddr 08:00:27:88:33:a4  

          inet addr:10.0.2.15  Bcast:10.0.2.255  Mask:255.255.255.0

          inet6 addr: fe80::a00:27ff:fe88:33a4/64 Scope:Link

          UP BROADCAST RUNNING MULTICAST  MTU:1500  Metric:1

          RX packets:122 errors:0 dropped:0 overruns:0 frame:0

          TX packets:182 errors:0 dropped:0 overruns:0 carrier:0

          collisions:0 txqueuelen:1000 

          RX bytes:43296 (43.2 KB)  TX bytes:23263 (23.2 KB)



lo        Link encap:Local Loopback  

          inet addr:127.0.0.1  Mask:255.0.0.0

          inet6 addr: ::1/128 Scope:Host

          UP LOOPBACK RUNNING  MTU:16436  Metric:1

          RX packets:18 errors:0 dropped:0 overruns:0 frame:0

          TX packets:18 errors:0 dropped:0 overruns:0 carrier:0

          collisions:0 txqueuelen:0 

          RX bytes:1642 (1.6 KB)  TX bytes:1642 (1.6 KB)

If I'm reading this correctly, the second entry is for "local loopback", that is, it doesn't go through the virtual network card. And the first one uses eth0. So that might make a difference (Ethernet v. Wi-Fi)

Also, the ip address is 10.0.2.15, which I think should not be visible beyond the "router"---the VM.

The third comment is from anthonysomerset:

is apache in your vhost configured to listen on the correct ports? this will work in bridged or nat mode, also check any firewall rules in the guest and that network access works out of the guest as well.

So:
- this will work in bridged or NAT
- make sure the Listen directive is correct
- check the firewall
- check that network access works out of the guest

The guest is set up to listen on port 80. That is the standard port set up in Apache. I wrote a Python script to filter comments (filter.py) and do:

python filter.py /etc/apache2/apache2.conf



Include mods-enabled/*.load

Include mods-enabled/*.conf

Include httpd.conf

Include ports.conf

..

Include conf.d/

Include sites-enabled/

So there are lots of places to look for possible conflicts! Luckily, httpd.conf is empty. ports.conf:

..

NameVirtualHost *:80

Listen 80

..

I'm not quite sure about the first one, but the Listen part is correct.

Firewall rules are here (at least, some are):

te@VB:/etc/apache2/sites-available$ python ~/Desktop/filter.py default

<VirtualHost *:80>

	ServerAdmin webmaster@localhost

	DocumentRoot /var/www

	<Directory />

		Options FollowSymLinks

		AllowOverride None

	</Directory>

	<Directory /var/www/>

		Options Indexes FollowSymLinks MultiViews

		AllowOverride None

		Order allow,deny

		Allow from all

		# Allow from 127.0.0.0/255.0.0.0 ::1/128

	</Directory>

..

Skipping the script stuff, this has been edited from my previous rule restricting access to go back to Allow from all.

About checking network access from the guest: if he means accessing a server on the host machine, I haven't tried that yet. But of course, Firefox works.

So, out of all this, the only thing I can see is that I don't know what this stuff is about.

NameVirtualHost *:80

<VirtualHost *:80>

I tried switching from Wi-Fi to Ethernet in the host, but that didn't help.

I'm going to have to look through all the other posts related to this topic. But it seems like it's getting to be too hard, after all, I used VirtualBox to make things easy.

Keeping up appearances

2011-08-13T07:36:00.001-04:00

Something's changed in blogger that messes with my code examples. The last normal one is from Sat August 6, and the first unusual one is from last Thu August 11. It looks the same as ever in Preview mode, but after publishing, all the code text has acquired an extra newline and is double-spaced.

Code is enclosed within tags:



<div style="overflow-x: scroll ">

<table bgcolor="#ffffb0" border="0" 

  width="100%" padding="4">

<tbody>

<tr>

<d>

<pre style=" hidden;font-family:monaco;">

..

</pre>

</table>

</div>

I liked this method because the output is monospaced, and the background colored according to the type of code (Python, R, output, etc).

I guess I'll have to look into CSS. If anyone has a good set of instructions to point me to, I'd be grateful.

Python C Extension

2011-08-12T08:30:00.005-04:00

Following Andrew Dalke's example, this is perhaps the simplest possible C extension for Python. For explanation of what's happening, see his page.

[ Something's gone wacko with this post---it's single spaced in preview, double after posting ]

x.c



#include <math.h>

#include "Python.h"



int process (const char *c, double f) {

    int i = 0;

    if (c[0] == 'c') {

        i = floor(f);

    }

    return i;

}



PyDoc_STRVAR(x__doc__,

"x module for process'ing stuff");



PyDoc_STRVAR(process__doc__,

"c,f -> do something with f if c == 'c'");



static PyObject *

py_process(PyObject *self, PyObject *args) {

    double f = 0;

    char *c;

    int i;

    if (!PyArg_ParseTuple(args, "sd:process", &c, &f))

    		return NULL; 

    i = process(c,f);

    return PyInt_FromLong((long) i);

}



static PyMethodDef x_methods[] = {

	{"process",  py_process, METH_VARARGS, process__doc__},

	{NULL, NULL}      /* sentinel */

};



PyMODINIT_FUNC

initx(void)

{

	Py_InitModule3("x", x_methods, x__doc__);

}

setup.py

from distutils.core import setup, Extension



setup(name="x", version="0.0",

	ext_modules = [Extension("x", ["x.c"])])

> python setup.py build

running build

running build_ext

building 'x' extension

creating build

creating build/temp.macosx-10.6-universal-2.6

gcc-4.2 -fno-strict-aliasing -fno-common -dynamic -DNDEBUG -g -fwrapv -Os -Wall -Wstrict-prototypes -DENABLE_DTRACE -arch i386 -arch ppc -arch x86_64 -pipe -I/System/Library/Frameworks/Python.framework/Versions/2.6/include/python2.6 -c x.c -o build/temp.macosx-10.6-universal-2.6/x.o

creating build/lib.macosx-10.6-universal-2.6

gcc-4.2 -Wl,-F. -bundle -undefined dynamic_lookup -arch i386 -arch ppc -arch x86_64 build/temp.macosx-10.6-universal-2.6/x.o -o build/lib.macosx-10.6-universal-2.6/x.so

You can control which architecture is build via setup.py with:

export ARCHFLAGS="-arch x86_64"

> ln -s build/lib.macosx-10.6-universal-2.6/x.so x.so

>

> python -c "import math;  import x;  print x.process('c',math.pi)"

3

> python -c "import math;  import x;  print x.process('a',math.pi)"

0

> python

Python 2.6.1 (r261:67515, Jun 24 2010, 21:47:49) 

[GCC 4.2.1 (Apple Inc. build 5646)] on darwin

Type "help", "copyright", "credits" or "license" for more information.

>>> import x

>>> x.__doc__

"x module for process'ing stuff"

>>> x.process.__doc__

"c,f -> do something with f if c == 'c'"

>>>

Looks good.

Note: on my 64-bit machines I'm getting a warning:

x.c: In function ‘process’:

x.c:7: warning: implicit conversion shortens 64-bit value into a 32-bit value

Converting i to long doesn't help. I'll try to look into that.

Trying Ubuntu Linux (3)

2011-08-12T08:00:00.001-04:00

I have an idea that someday I might try to teach students how to set up a web site (for example, EMBOSS on a local network would be nice). Python would be a perfect fit and it should be fun and extensible. I would probably want a neutral platform to do that (since most would be trained on Windows and I know nothing except OS X), so I looked into adding Apache Server to my Ubuntu 11.0.4 Desktop install (posts here and here).

There is a nice guide on the web from Ubuntu (here), and some simple instructions are here, and the Apache docs are here. I could have started by grabbing the Server version of Ubuntu, but instead I did:

sudo apt-get install apache2

As easy as that. Start and stop commands:

sudo /etc/init.d/apache2 start

sudo /etc/init.d/apache2 stop

sudo /etc/init.d/apache2 restart

To test it, just use Firefox and point it at localhost or 127.0.0.1. The index page is at /var/www/index.html (the "document root" is /var/www). If you want to get a little fancier, you can grab PHP:

sudo apt-get install php5 libapache2-mod-php5

sudo /etc/init.d/apache2 restart

And put this in /var/www/test.php:

<?php phpinfo(); ?>

Go to localhost/test.php and it should print a bunch of details about your setup.

I spent a lot of time working on access control (Apache docs here). I should probably have read the VirtualBox guide first (here):

A virtual machine with NAT enabled acts much like a real computer that connects to the Internet through a router. The "router", in this case, is the VirtualBox networking engine, which maps traffic from and to the virtual machine transparently. The disadvantage of NAT mode is that, much like a private network behind a router, the virtual machine is invisible and unreachable from the outside internet; you cannot run a server this way unless you set up port forwarding (described below).

Actually, we will want to do this later. A lot of time was wasted because I followed the "simple" instructions (e.g. Apache's) and couldn't figure out why stuff didn't work.

It turns out that although /etc/apache2/apache2.conf and /etc/apache2/httpd.conf have relevant settings there are a bunch of other files (in directories under /etc/apache2/). If you want to change access to the document root or scripts, you have to modify /etc/apache2/sites-available/default. Guess I should have read the first page (here) of the Ubuntu Server docs. It's right there!

Here is part of the file in the original form:

	ScriptAlias /cgi-bin/ /usr/lib/cgi-bin/

	<Directory "/usr/lib/cgi-bin">

		AllowOverride None

		Options +ExecCGI -MultiViews +SymLinksIfOwnerMatch

		Order allow,deny

		Allow from all

	</Directory>

I'm not sure what all of this does, but using ScriptAlias means that (Apache docs again):

The ScriptAlias directive tells Apache that a particular directory is set aside for CGI programs. Apache will assume that every file in this directory is a CGI program, and will attempt to execute it, when that particular resource is requested by a client.

I just modified /usr/lib/cgi-bin/ to be /home/te/cgi-bin/.

And, since we will ultimately care about access permissions, you should read (at least) this page of the docs, which explains that

		Order allow,deny

is explained as:

Allow,Deny
First, all Allow directives are evaluated; at least one must match, or the request is rejected. Next, all Deny directives are evaluated. If any matches, the request is rejected. Last, any requests which do not match an Allow or a Deny directive are denied by default.
Deny,Allow
First, all Deny directives are evaluated; if any match, the request is denied unless it also matches an Allow directive. Any requests which do not match any Allow or Deny directives are permitted.

I copied from further down in the file (about access to /usr/share/docs) to the ScriptAlias directive above:

        Order deny,allow

        Deny from all

        Allow from 127.0.0.0/255.0.0.0 ::1/128

The subnet masking stuff will need an explanation of its own.

So now, I can put the first script from the second part of this post (that prints data from os.environ in cgi-bin under my home directory, and it works.

Trying Ubuntu Linux (2)

2011-08-11T16:04:00.004-04:00

I'm installing software on Ubuntu 11.0.4 (in a Virtual Box VM running under OS X Lion) and my standard test is matplotlib. To begin with, Ubuntu 11.0 has python but does not have numpy. So, let's start there.

I found build instructions here, but didn't actually follow them. Instead, I took their advice and used apt-get. While I was still thinking I would build everything I did:

sudo apt-get install gfortran

sudo apt-get install python-dev

And then I went after BLAS, and LAPACK. Based on the library name on this page, I did:

sudo apt-get install libblas3gf

And using the library name from this page, I did:

sudo apt-get install liblapack3gf

Then, I just did:

sudo apt-get install python-numpy

And, amazingly enough, it works:

>>> import numpy

>>>

Now:

sudo apt-get install zlib-bin

sudo apt-get install libpng3

sudo apt-get install libfreetype6

sudo apt-get install python-matplotlib

Failed to fetch http://us.archive.ubuntu.com/ubuntu/pool/main/t/tcl8.5/tcl8.5_8.5.9-2_amd64.deb  Something wicked happened resolving 'us.archive.ubuntu.com:http' (-5 - No address associated with hostname)

E: Unable to fetch some archives, maybe run apt-get update or try with --fix-missing?

Oops..

sudo apt-get update

And now it gets tcl8.5

What else to do for matplotlib?
Nothing!

We'll try the example from here, and we get a png in the right place, which I just double click and get the graphic at the top. Here is the whole Desktop:

I really can't take any credit for this. I am ~~very impressed~~ speechless with the Ubuntu dev guys.

[ UPDATE: Did the same apt-get for scipy and the test at the bottom of this post works too! ]

Trying Ubuntu Linux (1)

2011-08-11T08:00:00.000-04:00

I thought I would play around with Linux a bit, trying to improve my Unix skills. I'm running it in a second partition on my Macbook that has OS X Lion installed (post here) and runs a Virtual Box VM.

I first heard about and tried Oracle's Virtual Box in connection with QIIME (post here).

I grabbed a VirtualBox version 4.1 binary installer for Mac OS X here.

Then I downloaded Ubuntu Linux 11.04 (now why is it called Natty Narwhal?) from here. There is another download page here but I think it has the same choices.

I ended up trying both Desktop "CD"s, because I was a little confused about the version I needed. There appears to be only 32-bit for i386 (ubuntu-11.04-desktop-i386.iso) but 64-bit for AMD (ubuntu-11.04-desktop-amd64.iso). However, the chip running the host OS is not relevant. The QIIME Linux is labeled AMD64 and is 64-bit and can run on my 64-bit machines under a VM. I suppose it would matter if I wanted to run Linux directly in the partition instead of in Virtual Box.

I have to confess, I've installed several different Linux versions (and then removed them) and the one described here a number of times, because I ran into problems with the install that I couldn't figure out how to back out of. But experimentation is relatively cheap, since it's just a VM. I don't have any worries about leftover crap on my disk. If worst comes to worst I can just erase the whole partition. Although it takes patience for the install steps.

Stock 64-bit Ubuntu Linux

So, I ran the Virtual Box installer, set the Name as Ubuntu, Operating System as Linux and the Version as Ubuntu 64. Then I set up the VM with all standard options (except 1 GB for RAM). Just click on through.

To actually install Linux using a .iso file, use the "First Run Wizard" by double-clicking the new machine in the list. For Select Installation Media: I navigated to the saved file. On one try I somehow got past the Wizard without properly installing and had to delete and re-configure the VM.

It takes quite a while to boot the first time, copying all the files over to the VM, I guess. Choose Erase disk and install Ubuntu (the disk size should be only the size of the VM disk). There are a bunch of dialogs to go through. After the install and reboot we get a message that Unity won't run on our hardware and so we're going to get Classic. Oh well.

Reboot..

Now, there are two things that seem really important to have but are not present in this stock install:

• a Shared Folder for file sharing between host and guest OS
• a Clipboard for copy/paste between host and guest OS

To have these features you must install what are called the "Guest Additions," which are additions to the Linux kernel. Some instructions are on this page:

After you have installed the Linux Guest system, there are a few additional packages needed for the Guest Additions. These packages enables you to create kernel modules. For Debian and Debian based distro's like Ubuntu, you need the following packages:
dkms
build-essential
linux-headers-generic

Find the Terminal (why not just drag it to the menu bar) and do:

sudo apt-get install dkms build-essential linux-headers-generic

You'll have to type in the command, since copy/paste doesn't work yet :)
When it's done:

Reboot to make the changes effective..

Next, "mount" the Guest Additions "iso": in the VB menu (in OS X) under Devices, choose Install Guest Additions. It shows up as a CD on the Desktop. Then from the dialog choose the default that automatically runs the shell script. Or hit Open AutoRun Prompt > OK and then in Terminal do:

media/VBOXADDITIONS{version}/VBoxLinuxAdditions.run

After it's done, reboot again..

At this point, copy/paste between host/guest should work. You'll have to click twice, once to switch focus to the VM and once to switch focus again to whatever application is receiving the paste.

There are two different key combos for Linux and neither is CMD. Not sure if this is standard Linux or it is just since the left CMD key is a special key for the VM. To copy and paste within Terminal on Linux use CTL-SHIFT-x and CTL-SHIFT-v rather than CMD (for Text Edit just use CTL-x and CTL-v).

Now we can try to get shared folders. In the Virtual Box window, from the Devices menu choose Shared Folders. Add a folder (navigate the host OS to the desired path). Also, make a folder under Linux that we'll use to share with, say /home/telliott/Desktop/share

In the dialog, choose Automount (and Make Permanent, see below). Even after the specification of a shared folder, you still have to mount it within Linux. The Qiime guys do:

sudo mount -a

Unfortunately, it doesn't seem to work for me. (Not in the Qiime install either). So then I tried this:

sudo mount -t vboxsf -o uid=1000,gid=1000 shared ~/Desktop/share

The first argument after vboxsf above is the name of the folder on the host OS (as defined by us in the Devices menu above). The second argument is the name of the sharepoint --- the path to a folder on the Guest.

It didn't work at first, because it turned out, I had failed to check Automount above. And, if you want to make a change to the Shared Folder setting, it seems that you have to reboot to be able to actually do it, otherwise there's a message about the resource being in use.

Also, I never did figure out how to change Make Permanent. Any change or attempt to delete in the Shared Folders dialog (even after a reboot) gives an error message: VERR_PERMISSION_DENIED. There's a note in the link about putting the mount command in

/etc/rc.local

but I don't see how that's going to work since it needs to run as root? Luckily, the VM allows me to suspend Linux in its current state. And now, thinking that I need to re-run that command every time upon booting, I try it and get:

sudo mount -t vboxsf -o uid=1000,gid=1000 shared ~/Desktop/share

/sbin/mount.vboxsf: mounting failed with the error: Invalid argument

After re-boot, the same mount fails! I am starting to get angry with this thing.

I gave up and set up Dropbox access from within the VM. It works great.

Now, to build some software.

PyMOL (1)

2011-08-06T17:07:00.010-04:00

If you're interested in Python and bioinformatics you've undoubtedly heard of Pymol.

If you're at a place that cares about molecular biology (my institution does not), you may even have a site license. I have a copy that I got for educational use about 3 years ago, and still have some examples that I made for a class I taught. That version has a Cocoa GUI.

There is also an open source version that is up on SourceForge (here). I thought it would be a good challenge to try to build the open source Pymol.

[ Update: I think the Cocoa GUI one is MacPymol here). I just noticed this on the legacy page. ]

The dependencies are listed in the README as:

PyMOL has the following external dependencies:
  1. OpenGL
  2. GLUT library for OpenGL (or freeglut)
  3. Python (v 2.1 or better) compiled with threads support
  4. libpng (can compiled without it, but image saves won't work)
  5. Tcl/Tk (for the external GUI)
  6. Numerical Python (technically optional - now using numpy instead)
  7. Python megawidgets (Pmw -- only required for the "external" GUI)
  8. FreeType2 (if you want nice fonts)

Most of these are already present in Lion, so I didn't think it would be too difficult. What I failed to realize is that the Tcl/Tk in Lion is apparently messed up (brief SO discussion here).

I downloaded pymol-v1.4.1.tar.bz2 from the main page at SourceForge, and did the standard distutils thing:

python setup.py build

The first error was:

layer0/ShaderMgr.c:173: error: ‘GLEW_OK’ undeclared (first use in this function)

The top Google hit shows it's related to OpenGL and glut.
If present, there should be a file glew.h (and likely glut.h).

> find /usr -iname *glut*
/usr/X11/include/GL/glut.h
/usr/X11/include/GL/glutf90.h
/usr/X11/lib/libglut.3.7.dylib
/usr/X11/lib/libglut.3.dylib
/usr/X11/lib/libglut.dylib
/usr/X11/lib/pkgconfig/glut.pc
> find /usr -iname *glew*
> find /usr -iname *GLEW*
>

No glew.

brew install glew

works and we have


$ find /usr/local/lib -iname *glew*
/usr/local/lib/libGLEW.1.5.8.dylib
/usr/local/lib/libGLEW.1.5.dylib
/usr/local/lib/libGLEW.a
/usr/local/lib/libGLEW.dylib
/usr/local/lib/pkgconfig/glew.pc

and we find "GLEW_OK defined on line 15273 of glew.h. Unfortunately, the build fails with a link error at the end:

ld: warning: ignoring file /usr/local/lib/libGLEW.dylib, file was built for unsupported file format which is not the architecture being linked (i386)
ld: duplicate symbol __CShaderMgr in build/temp.macosx-10.7-intel-2.7/layer2/RepCylBond.o and build/temp.macosx-10.7-intel-2.7/layer2/RepCartoon.o for architecture i386

> file /usr/local/lib/libGLEW.dylib
/usr/local/lib/libGLEW.dylib: Mach-O 64-bit dynamically linked shared library x86_64

So it seems that Homebrew has built x86_64, which can't work with what pymol is building (i386). What I should do is rebuild libGLEW.dylib with the correct architecture. I haven't accomplished that yet. Instead, I noticed a page that discusses this very same problem. One answer is in a different thread (here):

Kasper, not sure if this will help you, but I think the GLEW stuff was
added recently in V1.4.x and it seems to be still somewhat experimental.
May be to get you started you could switch to v1.3r1 from the repository
and try to embedd that version. V1.3.x should not contain the new GL
stuff.

So let's try 1.3. I grabbed pymol-1.3r2-src.tar.bz2 from here and it got past the previous problem but then dies with:

/usr/X11R6/include/ft2build.h:56:38: error: freetype/config/ftheader.h: No such file or directory

It's looking in /usr/X11R6 and can't find freetype which is in /usr/X11 (my post here).

So.. I modify the setup.py in pymol:

                     
    inc_dirs=["ov/src",
              "layer0","layer1","layer2",
              "layer3","layer4","layer5", 
              "/usr/X11R6/include",
              # my additions 
               "/usr/X11/include",
               "/usr/X11/include/freetype2",
               "/usr/local/include/GL",

              EXT+"/include",
              EXT+"/include/GL",
              EXT+"/include/freetype2",
   "modules/cealign/src",
   "modules/cealign/src/tnt",
              ]
    libs=[]
    pyogl_libs = []
    lib_dirs=[]
    def_macros=[("_PYMOL_MODULE",None),
                ("_PYMOL_LIBPNG",None),
                ("_PYMOL_FREETYPE",None),
                ]
    ext_comp_args=[]
    ext_link_args=[
               "-L/usr/X11R6/lib", "-lGL", "-lXxf86vm",
               "-L"+EXT+"/lib", "-lpng", "-lglut", "-lfreetype",
               # my additions
               "-lGLEW"
               ]

And now it builds.

But.. after doing ./pymol the GUI hangs. I think it's Python (IDLE?), which I never use.

However, I can get pymol to run from the command line with no GUI:

> pymol -c
 PyMOL(TM) Molecular Graphics System, Version 1.3.
..
 Command mode. No graphics front end.
 Detected 2 CPU cores.  Enabled multithreaded rendering.
 PyMOL: normal program termination.

I ran a short script and made a figure, it looks fine.

So then, we get wise words from Ned Deily (here), and see (here)

- You should not rely on the Apple-suppled Pythons if you want to use IDLE.

I thought about putting Tcl/Tk 8.4 in a different place, etc.

But finally I said (f*** it) and I reformatted my drive to add a new partition, and then used my USB flash drive with the Lion Installer (post here), and then I installed MacPorts and did

$ sudo port install pymol
--->  Computing dependencies for pymol
--->  Dependencies to be installed: freetype bzip2 zlib glew libpng mesa makedepend pkgconfig glib2 autoconf help2man gettext expat libiconv gperf ncurses ncursesw p5-locale-gettext perl5 perl5.12 perl5 perl5 m4 automake libtool xorg-xproto py27-libxml2 libxml2 libxml2 python27 db46 gdbm openssl python_select readline sqlite3 python27 xorg-dri2proto xorg-glproto xorg-libXfixes xorg-fixesproto xorg-libX11 xorg-bigreqsproto xorg-inputproto xorg-kbproto xorg-libXau xorg-libXdmcp xorg-libxcb xorg-libpthread-stubs xorg-xcb-proto xorg-util-macros xorg-xcmiscproto xorg-xextproto xorg-xf86bigfontproto xorg-xtrans xorg-libXi xorg-libXext xorg-libXmu xorg-libXt xorg-libsm xorg-libice py26-numeric python26 py26-pmw py26-tkinter tk Xft2 fontconfig xrender xorg-renderproto tcl xorg-libXScrnSaver xorg-scrnsaverproto xdpyinfo xorg-libXinerama xorg-xineramaproto xorg-libXtst xorg-recordproto xorg-libXxf86vm xorg-xf86vidmodeproto xorg-libdmx xorg-dmxproto

Incredible! Everything but the kitchen sink. Perhaps even more incredible, it worked. An hour or two later, we finish with:


--->  Attempting to fetch pymol-1.4_2.darwin_11.x86_64.tbz2 from http://packages.macports.org/pymol
--->  Fetching pymol
--->  Verifying checksum(s) for pymol
--->  Extracting pymol
--->  Applying patches to pymol
--->  Configuring pymol
--->  Building pymol
--->  Staging pymol into destroot
--->  Installing pymol @1.4_2
--->  Activating pymol @1.4_2
--->  Cleaning pymol
$

I got 1LMB.pdb from the PDB, and I ran this script with pymol script.pml:

load ~/Desktop/1LMB.PDB
hide everything
select D1, chain 1
select D2, chain 2
select R1, chain 3
select R2, chain 4
deselect all
show sticks, D1
show sticks, D2
show cartoon, R1
show cartoon, R2
color red, R1
zoom active, -10
move x, 4
bg_color white
ray 1200,1200
png ~/Desktop/x.png

and I got what's at the top of the post. Pretty cool! More later on this topic.

Exploring ctypes (5)

2011-08-05T11:51:00.002-04:00

I was a little snarky about how the types docs start with all Windows examples, but once you get to here they are really excellent. So I won't say too much more about them except to reproduce a fun example:

>>> from ctypes import *
>>> class cell(Structure):
...     pass
... 
>>> cell._fields_ = [('name', c_char_p),
...                  ('next', POINTER(cell))]
>>> 
>>> c1 = cell()
>>> c1.name = 'foo'
>>> c2 = cell()
>>> c2.name = 'bar'
>>> c1.next = pointer(c2)
>>> c2.next = pointer(c1)
>>> p = c1
>>> for i in range(8):
...     print p.name,
...     p = p.next[0]
... 
foo bar foo bar foo bar foo bar

This linked list is a classic example of something that needs a "forward declaration" of the cell type in C.

One other small point, from the docs:

Either the LoadLibrary() method of the dll loaders should be used, or you should load the library by creating an instance of CDLL by calling the constructor:

We used the second method before. To use the first with the example from the other day (here):

>>> import ctypes
>>> libadd = ctypes.cdll.LoadLibrary('add1.so')
>>> libadd.f1(2)
3
>>> libadd.f1(2,3)
3
>>> libadd.f1()
1858614497
>>> libadd.f1.argtypes = [ctypes.c_int]
>>> libadd.f1(2,3)
3
>>> libadd.f1()
Traceback (most recent call last):
  File "", line 1, in 
TypeError: this function takes at least 1 argument (0 given)
>>>

You can get some error checking with the argtypes call, but not everything.

A developer of Cython gave a nice summary of why you might want to start with that in the SO post cited the other day. But for my small purposes, I think ctypes will probably be everything I'll need.

Note to self: .pth

2011-08-05T08:08:00.002-04:00

I think I finally got .pth files sorted out. I first ran into this with easy_install (here). The Python library docs are here.

A path configuration file is a file whose name has the form package.pth and exists in one of the four directories mentioned above; its contents are additional items (one per line) to be added to sys.path. Non-existing items are never added to sys.path, but no check is made that the item refers to a directory (rather than a file). No item is added to sys.path more than once. Blank lines and lines beginning with # are skipped. Lines starting with import (followed by space or tab) are executed.

The four directories are:

sys.prefix + lib/site-packages
sys.prefix + lib/site-python
sys.exec_prefix + lib/site-packages
sys.exec_prefix + lib/site-python

So, not all directories on the path are searched for .pth files. I made a nested directory:

/Users/telliott/Desktop/temp/x

and put a .pth file in site-packages:

> cat /Library/Python/2.6/site-packages/x.pth
/Users/telliott/Desktop/temp/x

To test it:

> python -c 'import sys;  L = [p for p in sys.path if "Desk" in p];  print "\n".join(L)'
/Users/telliott/Desktop/temp/x

Of course, I could have just added a path to PYTHONPATH in ~/.bash_profile:

export PYTHONPATH=$PYTHONPATH:$HOME/Desktop/temp

> python -c 'import sys;  L = [p for p in sys.path if "Desk" in p];  print "\n".join(L)'
/Users/telliott/Desktop/temp
/Users/telliott/Desktop/temp/x

You need to restart Terminal to see the effect of changes to PYTHONPATH.

It's worth emphasizing (as pointed out in the quote) that .pth files can contain code. Here is easy_install.pth:

> cat /Library/Python/2.6/site-packages/easy-install.pth
import sys; sys.__plen = len(sys.path)
./Sphinx-0.6.3-py2.6.egg
./docutils-0.6-py2.6.egg
./Jinja2-2.2.1-py2.6.egg
./Pygments-1.1.1-py2.6.egg
./rpy2-2.1.7_20101104-py2.6-macosx-10.6-universal.egg
./numpy-1.5.0-py2.6-macosx-10.6-universal.egg
./Cython-0.14-py2.6-macosx-10.6-universal.egg
./PIL-1.1.7-py2.6-macosx-10.6-universal.egg
import sys; new=sys.path[sys.__plen:]; del sys.path[sys.__plen:]; p=getattr(sys,'__egginsert',0); sys.path[p:p]=new; sys.__egginsert = p+len(new)

It seems roundabout, but easy_install.pth is first adding these paths to sys.path, and then moving them to the position in the path specified by sys.__egginsert.

SciPy on Lion (1)

2011-08-04T19:27:00.004-04:00

I was looking for new things to do today, and realized I haven't really done much with SciPy. Their page about installing it is here. It seems like the developers are tired of dealing with people who can't get it installed and recommend that you download a binary or at least, install the Python.org version of Python first.

Apple ships its own version of Python with OS X. However, we strongly recommend installing the official Python distribution.

Of course, I wasn't gonna do that. There are other recommendations: that you really ought to get readline (somehow better than what comes with OS X), that you should get and use gcc-4.0 (only gcc-4.2 on Lion that I could find), that you should set LDFLAGS and FFLAGS and so on. I ignored it all.

> svn co http://svn.scipy.org/svn/scipy/trunk scipy
> cd scipy
> python setup.py build

The build fails with

unable to execute swig: No such file or directory

> brew install swig
Error: No available formula for swig

SWIG installation (download here) was a bit of a hiccup:

configure: error: in `/Users/telliott/Desktop/swig-2.0.4':
configure: error: 
        Cannot find pcre-config script from PCRE (Perl Compatible Regular Expressions)
        library package. This dependency is needed for configure to complete,
        Either:
        - Install the PCRE developer package on your system (preferred approach).
        - Download the PCRE source tarball, build and install on your system
          as you would for any package built from source distribution.
        - Use the Tools/pcre-build.sh script to build PCRE just for SWIG to statically
          link against. Run 'Tools/pcre-build.sh --help' for instructions.
          (quite easy and does not require privileges to install PCRE on your system)
        - Use configure --without-pcre to disable regular expressions support in SWIG
          (not recommended).

./configure --without-pcre
make
sudo make install

Back to SciPy and python setup.py build and problems as described here. I'm not too swift with this svn/git stuff so rather than find the patched version I just patched the files:

add #include <complex.h> to each of:

scipy/linalg/src/fblaswrap_veclib_c.c
scipy/lib/blas/fblaswrap_veclib_c.c.src
scipy/sparse/linalg/eigen/arpack/ARPACK/FWRAPPERS/veclib_cabi_c.c

And it seems to have built:

> python
Python 2.7.1 (r271:86832, Jun 16 2011, 16:59:05) 
[GCC 4.2.1 (Based on Apple Inc. build 5658) (LLVM build 2335.15.00)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> from scipy import *
>>>

I'll look into to testing it tomorrow.

Exploring ctypes (4)

2011-08-04T15:20:00.003-04:00

Introducing the subject of ctypes, I want to point you to a neat example from Andrew Dalke (second half of post here) that is a little beyond me, in many respects. But I wanted to give a shout out, and also say that I got it to run with a minor tweak.

• in the example file integrators.c, it should read
double integrate_gq(probe_func_t f) {
..
double integrate_rectangle(probe_func_t f) {

• the sym link for the built library should be (on 10.7):
ln -s build/lib.macosx-10.7-intel-2.7/integrators.so

And here is the output I got running his examples by cut and paste. First the build:

> python setup.py build
running build
running build_ext
building 'integrators' extension
llvm-gcc-4.2 -fno-strict-aliasing -fno-common -dynamic -g -Os -pipe -fno-common -fno-strict-aliasing -fwrapv -mno-fused-madd -DENABLE_DTRACE -DMACOSX -DNDEBUG -Wall -Wstrict-prototypes -Wshorten-64-to-32 -DNDEBUG -g -fwrapv -Os -Wall -Wstrict-prototypes -DENABLE_DTRACE -arch i386 -arch x86_64 -pipe -I/System/Library/Frameworks/Python.framework/Versions/2.7/include/python2.7 -c integrators.c -o build/temp.macosx-10.7-intel-2.7/integrators.o
creating build/lib.macosx-10.7-intel-2.7
llvm-gcc-4.2 -Wl,-F. -bundle -undefined dynamic_lookup -Wl,-F. -arch i386 -arch x86_64 build/temp.macosx-10.7-intel-2.7/integrators.o -o build/lib.macosx-10.7-intel-2.7/integrators.so

And the exploration in Python. We load the library:

> python
Python 2.7.1 (r271:86832, Jun 16 2011, 16:59:05) 
[GCC 4.2.1 (Based on Apple Inc. build 5658) (LLVM build 2335.15.00)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import ctypes
>>> ctypes.CDLL("integrators.so")
<h;CDLL 'integrators.so', handle 7fceb0483980 at 1056d9a90>

We fail to get our C integrator function to accept the Python-defined "probe" function as argument:

>>> integrators = _
>>> integrators.integrate_rectangle
<_FuncPtr object at 0x105759050>
>>> def flat_probe(x):
...     return 1.0
... 
>>> integrators.integrate_rectangle(flat_probe)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ctypes.ArgumentError: argument 1: <type 'exceptions.TypeError'>: Don't know how to convert parameter 1

[ About _ in Python: here and here ]

A simple cytypes.CFUNCTYPE definition does the trick:

>>> PROBE_FUNC = ctypes.CFUNCTYPE(ctypes.c_double, ctypes.c_double)
>>> 
>>> integrators.integrate_rectangle.argtypes = [PROBE_FUNC]
>>> integrators.integrate_rectangle.restype = ctypes.c_double
>>> integrators.integrate_rectangle(PROBE_FUNC(flat_probe)) 
2.0
>>> integrators.integrate_gq.argtypes=[PROBE_FUNC]
>>> integrators.integrate_gq.restype = ctypes.c_double
>>> integrators.integrate_gq(PROBE_FUNC(flat_probe))
2.0

And the original example

>>> import math
>>> def probe_function(x):
...     return math.cos(x)*(x+1)
... 
>>> integrators.integrate_gq(PROBE_FUNC(probe_function))
1.6829416883812192
>>> integrators.integrate_rectangle(PROBE_FUNC(probe_function))
1.7006012905844656
>>>

Where I'm going:

• ctypes library documentation
• learn more about distutils

Getting that build from these four lines is simply amazing:

setup.py

from distutils.core import setup, Extension

setup(name="integrators", version="0.0",
 ext_modules = [Extension("integrators", ["integrators.c"])])

And for reference, in case the page ever goes away here are his code listings:

integrators.h

#ifndef DALKE_INTEGRATORS_H
#define DALKE_INTEGRATORS_H

/* Two simple numerical integrators */

typedef double (*probe_func_t)(double);

double integrate_gq(probe_func_t f);
double integrate_rectangle(probe_func_t f);

#endif

integrators.c

/* Two numerical integrators */
#include "integrators.h"

/* 4 point Gaussian quadrature integration */
double integrate_gq(probe_func_t f) {
 return 
  0.652145155*(f(+0.339981044) + f(-0.339981044)) +
  0.347854845*(f(+0.861136312) + f(-0.861136312)); 
}

/* 4 point rectangle rule */
double integrate_rectangle(probe_func_t f) {
 return (f(-0.75)+f(-0.25)+f(0.25)+f(0.75))/2;
}

Exploring ctypes (3)

2011-08-04T08:00:00.000-04:00

Continuing with ctypes---here is an example of scoring a genome for potential sites for a DNA-binding protein according to the methods we talked about here.

We have a single function in a C module that will be compiled into a shared dynamic library. Then we load that library from Python, using ctypes, and call the function. The heavy lifting is done by the function score. The arguments are:

• the DNA sequence as a string (a str on the Python side)
• an array of floats which holds the scoring table
(in order acgt for each position in the site)
• the length of the site (n = 2 in the example)
• an array to hold the result for each position in the sequence

other variables:

• loop counters / index variables i, j, k
• the current character, c, so we can print it
• the current score for this char, d, also so we can print it
• the current score for this position in the sequence, r

The logic is straightforward and elementary. (Although perhaps there are faster approaches that people might suggest?)

mylib.c

#include <stdio.h>
#include <string.h>

void score(const char* dna, double p[], int n, double result[]) {
    char c;
    int i,j,k,N;
    double d,r;
    
    N = strlen(dna)-n+1;
    for (i=0; i<N; i++){
        r = 0;
        for (j=0; j<n; j++){
            c = dna[i+j];
            switch (c) {
                case 'a': { k=0; break; }
                case 'c': { k=1; break; }
                case 'g': { k=2; break; }
                case 't': { k=3; break; }
            }
            printf("%c ", c);
            d = p[j*4 + k];
            printf("%3.3f ",d);
            r += d;
        }
        printf("%3.3f\n", r);
        result[i] = r;
    }
}

We build the shared library using clang (or gcc, a quick test showed no speed difference):

clang -g -Wall -c mylib.c
clang -dynamiclib -current_version 1.0  mylib.o  -o mylib.dylib

The Python script is a little more unusual. The array initialization was just adapted from the docs (here), I don't yet understand how it works. (It seems like a pointer would be better).

Here is the second of those:

NFloats = ctypes.c_double * N
result = [0.0] * N
result = NFloats(*result)

And I was surprised that we don't have to worry about the types for arguments to our score function---that's why those lines are commented out at the top.

For this toy example, the DNA sequence is 14 nt, and the scoring table is given in the list L. Here is the script:

script.py

import ctypes
pre = '/Users/telliott/Desktop'
mylib = ctypes.CDLL(pre + '/mylib.dylib', ctypes.RTLD_GLOBAL)
#mylib.score.argtypes = [ctypes.c_char_p]
                        #ctypes.c_float_p,
                        #ctypes.c_int]
                        
dna = 'actgtcgactcgag'
L = [ 0.567,  -1.603,  -0.2245, 0.3605,
     -0.1175, -0.4655, -0.5326, 0.6898 ]
     
EightFloats = ctypes.c_double * len(L)
ff = EightFloats(*L)
n = int(len(L)/4)
N = len(dna) - n + 1

NFloats = ctypes.c_double * N
result = [0.0] * N
result = NFloats(*result)

mylib.score(dna, ff, n, result)
for f in result[:6]:  print round(f,3)

Here's the output. Each position is evaluated on a separate line. The last six values printed are the results received back in Python.

> python script.py
a 0.567 c -0.466 0.101
c -1.603 t 0.690 -0.913
t 0.360 g -0.533 -0.172
g -0.225 t 0.690 0.465
t 0.360 c -0.466 -0.105
c -1.603 g -0.533 -2.136
g -0.225 a -0.117 -0.342
a 0.567 c -0.466 0.101
c -1.603 t 0.690 -0.913
t 0.360 c -0.466 -0.105
c -1.603 g -0.533 -2.136
g -0.225 a -0.117 -0.342
a 0.567 g -0.533 0.034
0.101
-0.913
-0.172
0.465
-0.105
-2.136

Not bad for 50 lines of code and a little more than an hour of work! If you want to see it go fast, comment out the printf statements in the C code, and change:

dna = 'actgtcgactcgag'*400000

Now the DNA length is 5.6E6, and it finishes in about 5 seconds. It's impressive.

Exploring ctypes (2)

2011-08-03T20:39:00.003-04:00

I used an example from before (here) to make a dylib and then called the functions using ctypes. It seems very easy:

add1.c

int f1(int x) { return x + 1; }

add2.c

int f2(int x) { return x + 2; }

> gcc -g -Wall -c add*.c
> gcc -dynamiclib -current_version 1.0  add*.o  -o libadd.dylib
> file libadd.dylib
libadd.dylib: Mach-O 64-bit dynamically linked shared library x86_64

> python
Python 2.7.1 (r271:86832, Jun 16 2011, 16:59:05) 
[GCC 4.2.1 (Based on Apple Inc. build 5658) (LLVM build 2335.15.00)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import ctypes
>>> pre = '/Users/telliott/Desktop'
>>> libadd = ctypes.CDLL(pre + '/libadd.dylib', ctypes.RTLD_GLOBAL)
>>> x = libadd.f1(1)
>>> x
2
>>> y = libadd.f2(x)
>>> y
4
>>>

I guess one question is how to alter / determine the search path for ctypes in loading libraries. libc is in /usr/lib, but libpng is in /usr/X11/lib, as we discussed the other day (here):

>>> libc = ctypes.CDLL('libc.dylib', ctypes.RTLD_GLOBAL)
>>> libpng = ctypes.CDLL('libpng.dylib', ctypes.RTLD_GLOBAL)
Traceback (most recent call last):
  File "", line 1, in 
  File "/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/ctypes/__init__.py", line 353, in __init__
    self._handle = _dlopen(self._name, mode)
OSError: dlopen(libpng.dylib, 10): image not found
>>> pre = '/usr/X11/lib/'
>>> libpng = ctypes.CDLL(pre + 'libpng.dylib', ctypes.RTLD_GLOBAL)
>>>

>>> import ctypes.util
>>> ctypes.util.find_library('libpng.dylib')
>>> ctypes.util.find_library('png')
>>> ctypes.util.find_library('AGL')
'/System/Library/Frameworks/AGL.framework/AGL'

Obviously, find_library doesn't search /usr/X11/lib. Wonder if there is a way to tell it what we want?

After we look at Andrew's more sophisticated examples, we'll try something relevant to bioinformatics

Exploring ctypes (1)

2011-08-03T08:30:00.002-04:00

A perennial question in Python is what to do when Python isn't fast enough, and a related one is how to work with C libraries from Python. As discussed in this SO post, there are several possibilities, here I want to start working with ctypes.

This is now a Python built-in module. I was a bit put off when I looked at the official docs, since all the early examples are windows and there is almost no mention of OS X. (However, I will have to make sure I understand what's there later), I'd start instead with Andrew Dalke's quick tutorial (here). So the first part of this is my code duplicating a few of his examples.

Getting started is literally one line, preceded by the import, and followed by the function call itself. That's wonderful!

>>> import ctypes
>>> libc = ctypes.CDLL("libc.dylib", ctypes.RTLD_GLOBAL)
>>> print libc.strlen("Hello world!")
12
>>>

It's pretty self-explanatory, if you're at all familiar with C, except for the constant RTLD_GLOBAL. Under the hood, the library libc uses the functions dlopen and friends defined in clfcn.h. (We touched on these once before, in the last example of this post).

The prototype for the C call to load a library is:

void *dlopen(const char *file, int mode);

RTLD_GLOBAL is one of several enumerated modes, and the value we provide is passed in as an argument to dlopen. You can read about dlopen here.

Here is printf:

>>> libc.printf("the number %d and string %s\n", 10, 'spam')
the number 10 and string spam
30

I'm not sure why this works, because in the standard C library, printf is in libstdio (ref)

The function printf returns the number of characters printed, or a negative value to indicate an error. The next example introduces the idea that frequently we need to tell ctypes the type of input arguments or return values.

>>> libc.atoi('3')
3
>>> libc.atof('3.14')
3
>>> libc.atof.restype = ctypes.c_double
>>> libc.atof('3.14')
3.14

The assumed return type is int, and in order to get this example to work we have to specify that the returned value is a double. Here's the prototype:

double atof (const char *str);

Argument types are specified like this:

>>> libc.atan2.argtypes = [ctypes.c_double, ctypes.c_double]
>>> libc.atan2.restype = ctypes.c_double
>>> libc.atan2(3.0, 4.0)
0.6435011087932844
>>> import math
>>> math.atan2(3.0, 4.0)
0.6435011087932844

Finally, here's another function, this one from the math library:

>>> libmath = ctypes.CDLL("libm.dylib", ctypes.RTLD_GLOBAL)
>>> libmath.sqrt.restype = ctypes.c_double
>>> libmath.sqrt(2)
nan
>>> libmath.sqrt(2.0)
Traceback (most recent call last):
  File "", line 1, in 
ctypes.ArgumentError: argument 1: : Don't know how to convert parameter 1
>>>
>>> libmath.sqrt.argtypes = [ctypes.c_double]
>>> libmath.sqrt(2)
1.4142135623730951

And here is a glowing testimonial about controlling a camera via USB with code written in Python.

It'd be interesting to see what we could do on that in OS X.

Still to do:

• finish Andrew's tutorial (lots more there)
• read the official docs
• try to apply it to the problem of finding binding sites (here)

e again

2011-08-03T08:00:00.001-04:00

I came across a neat quote from Martin Gardner:

The number e is as ubiquitous as pi, turning up everywhere, and especially in probability theory. If you randomly select real numbers between 0 and 1, and continue until their sum exceeds 1, the expected number of choices is e.

-Martin Gardner
note in Calculus Made Easy, p. 153

Of course, this begs for a Python simulation:

import random

def mean(L):
    return float(sum(L))/len(L)

def oneTrial():
    counter = 1
    f = random.random()
    while f <= 1:
        counter +=1
        f += random.random()
    return counter

N = int(1E7)
L = [oneTrial() for i in range(N)]
print mean(L)

> python script.py 
2.7183671

A derivation would be nice, but that'll have to wait. I wonder about the variance.

pkg-config

2011-08-03T05:18:00.005-04:00

A reasonably brief post about pkg-config, which we used to help build matplotlib the other day (here).

The web page is here and there is a guide here.

As it says on the first page:

pkg-config is a helper tool used when compiling applications and libraries. It helps you insert the correct compiler options on the command line so an application can use gcc -o test test.c `pkg-config --libs --cflags glib-2.0` for instance, rather than hard-coding values on where to find glib (or other libraries). It is language-agnostic, so it can be used for defining the location of documentation tools, for instance.

And from the man page:

DESCRIPTION
The pkg-config program is used to retrieve
information about installed libraries in the
system. It is typically used to compile and
link against one or more libraries. Here is a
typical usage scenario in a Makefile:

program: program.c
cc program.c $(pkg-config --cflags --libs gnomeui)

pkg-config retrieves information about pack-
ages from special metadata files. These files
are named after the package, and has a .pc
extension. On most systems, pkg-config looks
in and
for these files. It will additionally look
in the colon-separated (on Windows, semicolon-
separated) list of directories specified by
the PKG_CONFIG_PATH environment variable.

The package name specified on the pkg-config
command line is defined to be the name of the
metadata file, minus the .pc extension. If a
library can install multiple versions simulta-
neously, it must give each version its own
name (for example, GTK 1.2 might have the
package name "gtk+" while GTK 2.0 has
"gtk+-2.0").

What's happening is that the matplotlib setup.py script uses pkg-config if it's present. The binary is in /usr/local/bin/pkg-config. It looks for a subdirectory of /usr/X11/lib (and /usr/local/lib) called pkgconfig that constains .pc files for packages that specify metadata about themselves to take advantage of pkgconfig functionality:

> cat /usr/X11/lib/pkgconfig/libpng.pc
prefix=/usr/X11
exec_prefix=${prefix}
libdir=${exec_prefix}/lib
includedir=${prefix}/include/libpng12

Name: libpng
Description: Loads and saves PNG files
Version: 1.2.44
Libs: -L${libdir} -lpng12
Libs.private: -lz 
Cflags: -I${includedir}

We used Homebrew only to grab pkg-config. You could get it yourself: 0.26 is the newest, but requires glib-2.0; 0.25 bundles glib with pkg-config. 0.25 builds with no dependencies using the standard ./configure; make; make install. Download page for version 0.25 and 0.26 is here. Or of course, you can use Homebrew like we did.

Not clear why it doesn't ship with OS X. The license is GPLv2. But then, so is BASH.

Some minor issues:

• when I tried using pkg-config on my iMac (Snow Leopard, littered with MacPorts and other stuff), matplotlib seemed to build and install but died with:

> python x.py
Traceback (most recent call last):
  File "x.py", line 1, in 
    import matplotlib.pyplot as plt
  File "/Library/Python/2.6/site-packages/matplotlib/pyplot.py", line 95, in 
    new_figure_manager, draw_if_interactive, _show = pylab_setup()
  File "/Library/Python/2.6/site-packages/matplotlib/backends/__init__.py", line 25, in pylab_setup
    globals(),locals(),[backend_name])
  File "/Library/Python/2.6/site-packages/matplotlib/backends/backend_macosx.py", line 243, in 
    class TimerMac(_macosx.Timer, TimerBase):
AttributeError: 'module' object has no attribute 'Timer'

For reasons that I don't understand, the version of the file in the build directory

/Users/telliott/Software/matplotlib-1.0.1/build/lib.macosx-10.6-universal-2.6/matplotlib/backends/backend_macosx.py

is not the same as the one in site-packages

/Library/Python/2.6/site-packages/matplotlib/backends/backend_macosx.py

Simply replacing the site-packages version fixed the problem:

> sudo cp /Users/telliott/Software/matplotlib-1.0.1/build/lib.macosx-10.6-universal-2.6/matplotlib/backends/backend_macosx.py /Library/Python/2.6/site-packages/matplotlib/backends/backend_macosx.py
> python -c "import matplotlib.pyplot as plt"
>

• the Homebrew "Formula" in /usr/local/Library/Formula/pkg-config.rb tells the installed pkg-config which paths to search for lib.pc files:

        /usr/local/lib/pkgconfig
        /usr/lib/pkgconfig
        /usr/X11/lib/pkgconfig

I'm not quite sure how that would work for a manual install.

• the libpng in /usr/local/lib that I installed shadowed the one in /usr/X11/lib so to use the X11 version, I removed the files from /usr/local

• The notation libpng.pcfileLibs.private: -lz indicates that libpng does depend on zlib, but it's already been linked against it (hence private). I'm not sure where it is, though.

matplotlib on OS X Lion--revised

2011-07-31T15:25:00.006-04:00

In the previous post, I repeated an installation of matplotlib (now for OS X Lion) following the same method that I've linked to a number of times on the blog. It worked this time as well, but lead me to consideration of other possible methods. I found one on the web here. So that's what this post is about.

BTW, whatever you do, do not follow the instructions that the matplotlib developers provide. You do not need another Python, including MacPython or the Endthought distribution or anything else..

To recap, matplotlib lists dependencies (in the make.osx file) of: libpng, libfreetype, and zlib. I read somewhere on the matplotlib site today (but can't find the link now) that zlib is not a required dependency. These days, the other two are actually provided by Apple:

> ls -al /usr/X11/lib/libpng*
-rwxr-xr-x  1 root  wheel  296864 Jul 29 16:11 /usr/X11/lib/libpng.3.dylib
lrwxr-xr-x  1 root  wheel      14 Jul 29 16:11 /usr/X11/lib/libpng.dylib -> libpng15.dylib
-rwxr-xr-x  1 root  wheel  294160 Jul 29 16:11 /usr/X11/lib/libpng12.0.dylib
-rwxr-xr-x  1 root  wheel  305008 Jul 29 16:11 /usr/X11/lib/libpng15.15.dylib
lrwxr-xr-x  1 root  wheel      17 Jul 29 16:11 /usr/X11/lib/libpng15.dylib -> libpng15.15.dylib

> ls -al /usr/X11/lib/libfreetype*
-rwxr-xr-x  1 root  wheel  1060416 Jul 29 16:11 /usr/X11/lib/libfreetype.6.dylib
lrwxr-xr-x  1 root  wheel       19 Jul 29 16:11 /usr/X11/lib/libfreetype.dylib -> libfreetype.6.dylib

I'm not sure at the moment whether X11 came with Xcode (as it used to) or was present in the original Lion install.

In any case, I spent an hour or two trying to figure out how to use the build commands from the makefile that comes with matplotlib, but point at these libraries. In the process I found what appears to be a bug:

/usr/X11/include/png.h:666: error: forward declaration of ‘struct png_info_def’

(and see Prashant's answer here), but we're going to use a different approach so it doesn't matter. I'm not sure how our method solved this in the end. As mentioned, the approach is to use Homebrew, in particular, something called pkgconfig. I used a Ruby script that downloads and installs Homebrew, as described here. The install gave this warning:

Warning: The following *evil* dylibs exist in /usr/local/lib
They may break builds or worse. You should consider deleting them:
/usr/local/lib/libfreetype.6.dylib
/usr/local/lib/libpng.3.dylib
/usr/local/lib/libpng12.0.dylib
/usr/local/lib/libz.1.2.5.dylib

These are, of course, the libraries we just installed in order to get matplotlib to build.

That's what got me working on the other method more seriously. I don't think the danger is all that great, but clearly it is better to use libraries that are (i) still available from the maintainers and (ii) have been vetted at least to some extent by other folks including Apple. AFAIK there aren't any known security issues with the versions we installed previously. I ripped them out anyway (they seemed to confuse pkgconfig).

If we ask Homebrew for the libraries, it just points us to the ones we saw in /usr/X11:

> brew search libpng
Apple distributes libpng with OS X, you can find it in /usr/X11/lib.
However not all build scripts look here, so you may need to call ENV.x11
in your formula's install function.

> brew search freetype
Apple distributes freetype with OS X, you can find it in /usr/X11/lib.
However not all build scripts look here, so you may need to call ENV.x11
in your formula's install function.

Following the instructions in the blog post I first installed pkgconfig:

> sudo brew install pkgconfig
Cowardly refusing to `sudo brew'
> brew install pkgconfig

!!

That's all the Homebrew we need. This is followed by:

cd ~/Desktop
git clone git://github.com/matplotlib/matplotlib.git
cd matplotlib/
python setup.py build
sudo python setup.py install

And we can see that I really did overwrite the first matplotlib install, and that we're actually using the libraries from /usr/X11, by first doing

export DYLD_PRINT_LIBRARIES=1

Then when we run a script that imports matplotlib.pyplot, the Terminal shows (among much else) this:

dyld: loaded: /usr/X11/lib/libfreetype.6.dylib
..
dyld: loaded: /usr/X11/lib/libpng15.15.dylib

So that's what I'd recommend and it seems to be working fine. This simple script works exactly as you'd expect.

import matplotlib.pyplot as plt
Y = [1,4,9,16]
plt.scatter(range(len(Y)),Y,s=250,color='r')
plt.savefig('example.png')

So the next thing to do is to figure out pkgconfig and Homebrew work their magic!

P.S. You will still need to make and edit ~/.matplotlibas discussed last time.

matplotlib on OS X Lion--old

2011-07-31T09:45:00.006-04:00

[ UPDATE: I'm going to leave this post up, for the record, but I found a better way, and that's in the following post. ]

Just a note to say that I got matplotlib installed with only a few issues to solve. As before, I relied on Gavin Huttley's instructions (here). Although the matplotlib instructions say you should install a different Python, I've found that just leads to confusion down the road, and I want to try to use the system version only if I can get away with it. YMMV.

The sticking points for the install always have to do with the basic dependencies, which are listed at the top of the file make.osx in the matplotlib source (from here). The source is the last link at the bottom.


ZLIBVERSION=1.2.3
PNGVERSION=1.2.39
FREETYPEVERSION=2.3.11

As silly as it seems, although OS X provides excellent functionality for dealing with compression, png images and fonts, we need to install these libraries for matplotlib to use. Unfortunately, I'm not expert enough to alter matplotlib to get rid of these dependencies, so we're stuck with it.

libpng

I got libpng-1.2.39 from here. I used this version b/c that's what matplotlib asks for, despite the fact that the current release is 1.5.4 (here). Standard magic incantations:

./configure
make
sudo make install

freetype

curl -O http://ftp.twaren.net/Unix/NonGNU/freetype/freetype-2.3.7.tar.bz2

According to Gavin's notes, we need these compiler flags for the other two (but they interfere with the libpng build):

export MACOSX_DEPLOYMENT_TARGET=10.7
export CFLAGS="-arch i386 -arch x86_64"
export FFLAGS="-arch i386 -arch x86_64"

./configure
make
sudo make install

zlib

zlib is a compression library. The current version is 1.2.5. I tried this:

curl -O http://www.zlib.net/zlib-1.2.3.tar.gz

but as I've seen before, although you get a (small) file back that has the right icon, it is really a 404 Page Not Found, which doesn't do me any good. It seems the older versions have been disappeared off the web and I couldn't find them poking around on the site. What I should've done at this point is try the version they do offer, but instead I went to my other machine and grabbed 1.2.3. I guess that since I'd used it before to build zlib, I found I needed to clean it first:

make clean
make
sudo make install

matplotlib

Finally, I got the matplotlib source. I altered the file make.osx:

PREFIX=/usr/local
PYVERSION=2.7
PYTHON=python${PYVERSION}
ZLIBVERSION=1.2.3
PNGVERSION=1.2.39
FREETYPEVERSION=2.3.7
MACOSX_DEPLOYMENT_TARGET=10.7
OSX_SDK_VER=10.7
ARCH_FLAGS="-arch i386-arch x86_64"

make -f make.osx mpl_build
sudo python setup.py install

mkdir ~/.matplotlib
cp matplotlibrc.template ~/.matplotlib/matplotlibrc

Edit the above file to have:

backend      : MacOSX

> python
Python 2.7.1 (r271:86832, Jun 16 2011, 16:59:05) 
[GCC 4.2.1 (Based on Apple Inc. build 5658) (LLVM build 2335.15.00)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> from matplotlib import *
>>>

Yes!

[ UPDATE: I see I used freetype 2.3.7 rather than 2.3.11, by mistake.

And if you don't like my approach, you could use Homebrew and follow this post. I don't have any experience with this yet. ]

OS X Lion

2011-07-30T10:42:00.006-04:00

OS X Lion has been released for a week now, so I thought I'd try it out on one of my machines. A definitive review is here. I'll be too busy to really play with it for a while, or to do much programming, but I thought I should at least try. An additional motivation is that my laptop (an aluminum MacBook from about December 08), had developed an issue. It completely drains the battery (starting from full charge) in sleep mode, even when no apps are running. The battery hasn't been particularly stressed:


  Health Information:
  Cycle Count: 568
  Condition: Normal
  Battery Installed: Yes
  Amperage (mA): 0
  Voltage (mV): 12546

And, it doesn't drain when powered off. This problem appeared recently, and I am interested to see if installing Lion might fix it, but that hasn't been tested yet. I'll update when I have those results.

This post is just to document that I got Lion onto a USB flash drive and used it to do the install and it worked fine. Instructions from here. I did this because it's the only way to do a clean install. One odd thing: if the machine on which you are attempting to do the download can't run Lion, you won't be allowed to do the download at all! I had to bring my laptop to work (and get the IT guys to reauthorize it for our network).

Also, I have two Apple IDs---an old one and one that goes with my MobileMe account. When I set up Mail, the ID for the computer was automatically set to be the latter, but I originally did the Lion download with the former. We'll see if and when that causes a problem.

It seems clear that the Mac OS will (in time) become more iOS-like (locked down, user not allowed to risk hurting himself, even if he wants to) and less easy to hack around with. But for now, it's still a win over Linux, for me.

Python is 2.7:


> python
Python 2.7.1 (r271:86832, Jun 16 2011, 16:59:05) 
[GCC 4.2.1 (Based on Apple Inc. build 5658) (LLVM build 2335.15.00)] on darwin
Type "help", "copyright", "credits" or "license" for more information.

The default scrolling is like that on the iPhone---as if you're moving the content, rather than a scrollbar. After 2 or 3 hours, I changed the setting back to the old method. Maybe I will try to retrain my brain at a later date.

I also downloaded Xcode 4. It quit in the middle for no reason, but completed on the second try. Weirdly, in the middle of the install I got this alert panel:

I say "weirdly" because iTunes is not running. Starting and quitting iTunes had no effect, either. Finally I went in Terminal and looked for an iTunes-related process:

  163 ??         0:00.03 /Applications/iTunes.app/Contents/MacOS/iTunesHelper.app/Contents/MacOS/iTunesHelper -psn_0_57358

after I killed it, the install finished fine. There is a mention in the big Ars Technica review that apps don't actually quit (or perhaps not always), and may not show the little dot for an active app, but there was no process for iTunes so I don't think so. It's weird that a bug like that should still be present at this stage of the game.

And I like the new looks in Terminal:

[ UPDATE: Yep, the battery does not drain on sleep anymore. (Well, 2% in 2 hr). That's the good news. The bad news is that the extension I had to help autofill passwords when forms contain autocomplete="off" doesn't work with Safari any more (post here, extension here). I sure hope Apple didn't do that on purpose. ]

UPDATE 2: I just needed to check to the right box in the Prefs.

Verification

2011-07-15T07:34:00.003-04:00

I posted two derivations for Euler's famous equation (here and here):

e^ix = cos x + i sin x

This can be verified in a particularly simple way.
The series representation of the exponential function:

e^x =  1 +  x + x²/2! +  x³/3! + x⁴/4! + ..

is especially neat because each term in the series is the derivative of the term following, and the result of that is:

d/dx e^x = e^x

Which is, indeed, one definition of this function.
Substitution of ix leads to a simple shift in the pattern:

e^ix = 1 + ix - x²/2! - ix³/3! + x⁴/4! + ..

repeating with a period of 4. Powers of x like:

4n + 1 are multiplied by  i
4n + 2                   -1
4n + 3                   -i
4n                        1

But remembering the series for sine and cosine, and multiplying the former by i:

  sin x =  x -  x³/3! + x⁵/5! - x⁷/7! + ..
i sin x = ix - ix³/3! + ..
  cos x =  1 -  x²/2! + x⁴/4! - x⁶/6! + ..

Adding:

cos x + i sin x = 1 + ix - x²/2! - ix³/3!+ x⁴/4! + ..
  = e^ix

Note on trig substitution

2011-07-14T17:04:00.003-04:00

A fact we needed in a recent post deriving Euler's formula is the value of a particular integral:

∫ dx / √(x² + 1) = ln (x + √(x² + 1))

This gave me a lot more trouble than expected, so I thought I'd work through it here for reference. The easiest way is to differentiate the answer. Let:

u = x + √(x² + 1)
du = [1 + 2x / 2 √(x² + 1)] dx
   = [1 + x / √(x² + 1)] dx

y = ln(u)
dy = 1/u du = [1 / (x + √(x² + 1)) ] [ 1 + x / (√(x² + 1)] dx

The trick is to notice that when we find the common denominator for the part at the far right and then add the terms, we generate the denominator of the part to the left.

dy = [1 / (x + √(x² + 1) ] [ (√(x² + 1) + x) / (√(x² + 1)] dx
   = 1 / (√(x² + 1)

The integrand is normally written more generally as:

∫ 1 / √(x² + a²)

and that doesn't change anything.

The forward version starts with a trigonometric substitution.

x = a tan y
y = tan^-1 (x/a)

Remembering that

(u/v)' = (vu' - uv')/v²
tan' = (cos² + sin²)/cos² = 1/cos² = sec²

dx = a sec² y dy

The other part of the integral is:

1 / √(x² + a²)

x = a tan y
x² = a² tan² y
x² + a² = a² tan² y + a²
  = a² ( 1 + tan² y)

If we start with the standard identity:

sin² y + cos² y = 1

and divide by cos² y:

tan² y + 1 = sec² y
x² = a² ( 1 + tan² y) = a² sec² y

Thus the integral reduces to:

∫ (a sec² y) / (a sec y)  dy = ∫ sec y dy

This small integral is itself a bit tricky. The answer, which I found on a really nice site (here), is to multiply top and bottom by:

       sec y + tan y
sec y  -------------
       sec y + tan y

Let u = sec y + tan y

If you work through it you'll see that

du = (sec y tan y + sec² y) dy

Thus, this is

∫ du/u = ln(u) = ln (sec y + tan y) 
  = ln [ ((√(x² + a²) + x ) / a ]

As I said, not so easy as other trig substitutions.

Euler's Gem 2

2011-07-14T09:30:00.005-04:00

Here is a sketch of a second derivation of Euler's famous formula:

e^iθ = cosθ + i sinθ

as presented by William Dunham in his book Euler, The Master of Us All. First post here.

The first step is to recall a standard trigonometric substitution in calculus:

y = sin x
x = sin^-1 y

√(1 - y²) = cos x

We're interested in the integral:

∫ dy / √(1 - y²)

Substituting with x we see that:

dy = cos x dx

And the integral is

∫ (1/cos x) cos x dx = ∫ dx = x
x = ∫ dy / √(1 - y²)

Now Euler makes a complex change of variable:

y = iz
x = ∫ dy / √(1 - y²)
  = ∫ i dz / √(1-(iz)²)
  = i ∫ dz / √(1 + z²)
  = i ln [√(1 + z²) + z]

The last step is another standard result from calculus which I will assume for the time being (more here).

Undo the substitution:

z = y/i = sin x / i
z² = -sin² x
√(1 + z²) = √(1 - sin² x)
  = cos x

x = i ln (cos x + sin x / i)

We will use two identities involving i:

u / i = - i u
1 / (cos u - i sin u) = (cos u + i sin u)

(For the second one, see the previous post). Now:

x = i ln (cos x + sin x / i)
x = i ln (cos x - i sin x)
ix = - ln (cos x - i sin x)
  = ln [ 1 / (cos x - i sin x) ]
  = ln (cos x + i sin x)

Just eponentiate:

e^ix = cos x + i sin x

Wow, again!

Euler's gem

2011-07-13T17:42:00.004-04:00

Here is a sketch of the derivation of Euler's famous formula:

e^iθ = cosθ + i sinθ

as presented by William Dunham in his book Euler, The Master of Us All.

The first part of the proof is similar to when we used Euler's formula to derive other formulas for trig functions of sums and differences of angles (post), only backward. Start from the definition of i:

i = √-1

To begin with, having i allows us to factor new expressions:

1 = cos² s + sin² s
  = (cos s + i sin s)(cos s - i sin s)

(I'm going to use s and t, as before, rather than θ and φ).

This shows where the original idea of cos + i sin comes from. (Of course, we could just as well do sin + i cos, that would result in a different convention for the orientation of the complex plane).

Suppose we have two angles s and t, we can multiply and then use the formulas from before (obtained by the geometric proof):

(cos s + i sin s)(cos t + i sin t) =
  = (cos s cos t - sin s sin t) + i(sin s cos t + cos s sin t)
  = cos(s + t) + i sin(s + t)

Set s = t:

(cos s + i sin s)² = cos(2s) + i sin(2s)

In fact Euler showed it works for fractional n but I'll assume that part:

[1] (cos s + i sin s)ⁿ = cos(ns) + i sin(ns)
n >= 1

If we multiply the difference rather than the sum:

(cos s - i sin s)(cos t - i sin t) =
  = (cos s cos t - sin s sin t) - i (sin s cos t + cos s sin t)
  = cos(s + t) - i sin(s + t)

Again, with s = t we have:

(cos s - i sin s)² = cos(2s) - i sin(2s)
[2] (cos s - i sin s)ⁿ = cos(ns) - i sin(ns)

Adding [1] and [2] we have:

2 cos(ns) = (cos s + i sin s)ⁿ + (cos s - i sin s)ⁿ

The middle part of the proof is where the magic happens. Let:

s = x/n

n -> ∞
s -> 0
cos s -> 1
sin s -> s

So..

cos x = cos ns = 1/2 [(cos s + i sin s)ⁿ + (cos s - i sin s)ⁿ]
cos x = 1/2 [(1 + is)ⁿ + (1 - is)ⁿ]
cos x = 1/2 [(1 + ix/n)ⁿ + (1 - ix/n)ⁿ]

But..

e^ix = (1 + ix/n)ⁿ 
as n -> ∞

cos x = 1/2 [e^ix + e^-ix]

By very similar manipulation to what's in the first part we can also handle the sine:

2i sin(ns) = (cos s + i sin s)ⁿ - (cos s - i sin s)ⁿ

We will obtain:

sin x = 1/(2i) [e^ix - e^-ix]

Now it's just a matter of addition:

cos x + i sin x = 1/2 [e^ix + e^-ix + e^ix - e^-ix]
  = 1/2 [e^ix + e^ix]
  = e^ix

Wow!

Rotating a hyperbola: general case

2011-06-27T07:00:00.002-04:00

Just to finish up quickly with the analysis of hyperbolas from the other day (here and here), following Stewart (pdf), we derived the following equations for x and y in terms of a coordinate system based on u and v (orthogonal, rotated through angle θ):

x = u cosθ - v sinθ
y = u sinθ + v cosθ

The general formula for a quadratic is:

Ax² + Bxy + Cy² + Dx + Ey + F = 0

When we transform to the new coordinates, we get:

Ax² = A[u² cos²θ - 2uv sinθ cosθ + v² sin²θ]
Bxy = B[u² cosθ sinθ + uv cos²θ - uv sin²θ - v² sinθ cosθ
Cy² = C[u² sin²θ + 2uv sinθ cosθ + v² cos²θ

Gathering the terms in uv we obtain the coefficients:

(C-A)(2sinθ cosθ) + B(cos²θ - sin²θ)

Remember the double angle formulas:

sin(s+t) = sin s cos t + cos s sin t
sin(2θ) = 2 sinθ cosθ

cos(s+t = cos s cos t - sin s sin t
cos(2θ) = cos²θ - sin²θ

So we obtain:

(C-A) sin(2θ) + B cos(2θ)

for the coefficients of xy. These must equal zero for all the xy terms to disappear. Thus:

tan(2θ) = B / (C-A)

This approach runs into a problem if C = A, as it does for our example:

xy = 1

But we can just invert the step at the end:

cot(2θ) = (C-A)/B = 0

The cotangent is zero when the cosine is zero, e.g. for

2θ = π/2

Thus, if

θ = π/4

then all the xy terms vanish, as we found before.

Rotational transformation: geometrical view

2011-06-26T13:34:00.006-04:00

A bit more about rotational transformations (see the previous post for the setup to this). We derived equations describing the coordinates of a point at x,y with respect to a new coordinate system in u,v that is rotated counter-clockwise by an angle θ:

x = u cos(θ) - v sin(θ)
y = u sin(θ) + v cos(θ)

We can also solve these equations for u and v. (It helps to know the answer). Notice that if we multiply the first equation by cos(θ) and the second one by sin(θ) and add, we'll end up with u by itself on the right-hand side:

x cos(θ) = u cos²(θ) - v sin(θ) cos(θ)
y sin(θ) = u sin²(θ) + v cos(θ) sin(θ)

x cos(θ) + y sin(θ) = u cos²(θ) + u sin²(θ) 
x cos(θ) + y sin(θ) = u

Similarly, if we multiply the first equation by sin(θ) and the second one by cos(θ) and subtract the first from the second, we end up with v by itself:

x sin(θ) = u sin(θ) cos(θ) - v sin²(θ)
y cos(θ) = u sin(θ) cos(θ) + v cos²(θ)

- x sin(θ) + y cos(&theta) = v cos²(θ) + v sin²(θ)
- x sin(θ) + y cos(&theta) = v

I came up with geometric explanations for these relationships, the one for x and y in terms of u and v is at the top of the post, and below is the same drawing but explaining u and v in terms of x and y.

A quicker way to the same place is to start by considering u and v to be the original coordinate system, and rotate clockwise through the angle θ. Then θ is negative, sin(θ) = -sin(-θ) and the cosine stays the same. We end up switching the signs of the sine terms, as we found.

Different views, same hyperbola

2011-06-26T10:56:00.009-04:00

In the algebra study guide that my son is using (ALEKS) they set up a hyperbola in two "standard" forms as:

(x-h)²/a² - (y-k)²/b² = 1
(y-k)²/b² - (x-h)²/a² = 1

or more generally:

(x-h)²/a² - (y-k)²/b² = +/- 1

Let's simplify our lives and consider an example centered at the origin:

x²/a² - y²/b² = 1

Then they divide the world into those hyperbolas that open up and down versus those that open up left and right.

By looking at the equation, I think we can agree that if x² equals 0 we've got a problem, since no value of y can satisfy the equation, and in fact x must not be less than 1. Hence, we can already predict that this system opens left and right, as the plot shows (a² = b² = 2):

The ALEKS review goes on to explain that for this type of equation, the vertices of the hyperbola (points of closest approach to the origin at h,k) are equal to

h +/- a, k

and the asymptotes have slopes equal to

+/- b/a

This all works great. The problem I had was that the simplest parabola I can think of is:

xy = 1

which doesn't fit the system.

The answer to my confusion is that any parabola may be rotated around its origin in the xy-plane. In 2 selected orientations its equation will have only x² and y² terms (when the vertices are on the x or the y-axis), whereas in 2 other selected orientations it may have only xy terms (when the vertices are on y= +/- x. The rest of the time it contains both.

This is explained in an excerpt from Stewart's Calculus I found here.

Consider a point P in the plane with coordinates x,y and distance from the origin r. Now, let's establish a new coordinate system u,v which is rotated counter-clockwise by the angle θ. A vector from the origin through P is rotated an angle φ with respect to the u axis and φ + θ with respect to the x-axis. Here's a screenshot from the pdf:

(Note, he uses capital X and Y for the second set of coordinates).

In the u,v system, the point P has coordinates

u = r cos(φ)
v = r sin(φ)

In the x,y system the coordinates are:

x = r cos(φ + θ)
y = r sin(φ + θ)

We remember (my post here) that:

sin(s+t) = sin s cos t + cos s sin t
cos(s+t) = cos s cos t - sin s sin t

Hence:

x = r cos(φ) cos(θ) - r sin(φ) sin(θ)
  = u cos(θ) - v sin(θ)

y = r sin(φ) cos(θ) + r cos(φ) sin(θ)
  = v cos(&theta) + u sin(θ)
  = u sin(&theta) + v cos(θ)

Now, consider the hyperbola:

xy = 1

Rewrite this in terms of u,v and multiply, giving 4 terms

+ u cos(θ) u sin(θ)
+ u cos(θ) v cos(θ)
- u sin(θ) v sin(θ)
- v sin(θ) v cos(θ)

If θ equals π/4, then

cos(θ) = sin(θ) = 1/√2

The two middle terms drop out and leave us with:

u²/2 - v²/2 = 1

If the middle terms don't cancel, we're left with a mixture including some fraction of
u²/2, v²/2, and u times v.

Note: we rotated the coordinate system counter-clockwise, which has the effect of rotating the plot clockwise, when the coordinate system is viewed in standard orientation.

One more thing about 3x3 determinants

2011-06-25T08:22:00.004-04:00

While we're on the topic of determinants (see yesterday's post), rather than pounding on arithmetic drills, the other thing I would show algebra students is that you can use any row or column to set up the computation. That can be useful if one choice contains zeros.

The standard approach uses the top row:

a b c
d e f
g h i

And the determinant is:

a(ei-hf) - b(di-fg) + c(dh-eg)

Where the sign of the second term is negative by the checkerboard rule:

+ - +
- + -
+ - +

Multiplying out we obtain:

aei - afh - bdi + bfg + cdh - ceg

The 6 terms contain three components, each taken from a different row and column. For example, the components of bdi are from:

b = (1,2)
d = (2,1)
i = (3,3)

The checkerboard rule makes the sign come out correctly.

If we're working with the top row or the middle column and so processing b x (di - fg) or d x (bi - ch), we'll need the minus sign; whereas if we're obtaining this term from i x (ae - bd) we already have a minus sign.

Let's try using the last row. We have:

a b c
d e f
g h i

g(bf - ce) - h(af - cd) + i(ae - bd)
bfg - cdg - afh + cdh + aei - bdi

Compare with the first example to see that all the terms are present.

Can we do it by the diagonal? Try ceg:

c(dh - eg) - e(ai - cg) + g(bf - ce)

Nope. Some terms are correct, but some are duplicates.