Revision as of 09:31, 3 May 2008

Contents

Sage and Cython: a brief introduction

Abstract

This is a quick introduction to Sage, a powerful new computational platform that builds on the strengths of Python. This article was directly inspired by Julius B. Lucks' "Python: All A Scientist Needs"; I recommend reading it first as it explains some of the attractions of Python and biopython. One thing I will highlight is using Cython[1] in Sage to make very fast code in an easy way.

Sage is a free and open-source project for computation of all sorts that uses Python as its primary language and "glue". One of the goals of Sage is to provide a viable free and open-source alternative to Matlab, Maple, and Mathematica. Sage unifies a great deal of open-source mathematical and statistical software; it includes biopython as an optional package and the statistics language R by default.

Sage notebook interface

The notebook starts by default in native Sage mode, which is almost identical to python. However, other environments are available. The screenshot below shows a Sage notebook in R mode, with some simple R commands executed [2].

The @interact command in the Sage notebook provides an easy way to make simple GUIs to explore data. In the example below, a user can enter the URL of a fasta-formatted protein file and a PROSITE-style regular expression. Using biopython and the "re" module of python we can search for and display matches to the pattern. For this screenshot, I used proteins from the malaria-causing Plasmodium falciparum and a fragment of the transthyretin pattern (Prosite PS00768). This is just for illustrative purposes - I do not claim any significance for these matches.

Here is another example of the interact command, along with some of Sage's 2D plotting. We take a human mitochondrial DNA sequence and plot the fraction of CG dinucleotides in a window of variable size. This fraction is divided by 16, to normalize it against the expected value in the case of independently, uniformly distributed nucleotides (obviously a simple model in this case).

Because Sage uses a browser as its primary GUI, you can continue working on a project remotely from anywhere in the world. It is also very easy to collaborate with others on a shared server. Similarly, you can make it easy for students to use Sage and work in groups through the notebook. I am teaching a class right now using Sage:

If I get a plea for help, I can dive right into a group's worksheet and help them figure it out.

Cython

Sage initially used an alternative to SWIG (described in Julius's article) called Pyrex to compile Python code to C when performance concerns demanded it. Because they needed to extend Pyrex in various ways, they created a friendly fork of Pyrex called "Cython". I believe it is fair to say that Cython is the easiest way to create C code in Python.

Lets look at a toy example of Cython usage. Consider the following python function, which we can imagine is the inner loop of some more complex operation:
<syntax type="python">
def slow_logistic(n,a,lamb = 3.82):

Time: CPU 6.17 s, Wall: 6.17 s
</syntax>
So the pure python version takes about 6 seconds for 10^6 iterations. If we tell Sage we'd like to use Cython, and declare variables to C types, we get quite a speedup:

So "Cythonizing" our loop sped it up by a factor of over 100. Often by converting a few crucial functions to Cython we get all the speed we need.

For a slightly more real-world example, suppose we want to interactively explore the behavior of a system of ODEs as we change parameters. The example below uses the "repressilator" system of Elowitz and Leibler [3], and a crude Euler's method (just for illustrative purposes!) to solve it. The behavior is not too bad for smaller runs but it gets a bit sluggish as we increase the simulation time.

<syntax type="python">
def m_der(m,p,a,a0,n):

return -m + a/(1+p^n) - a0

def p_der(m,p,b):

return -b*(p-m)

def der_vec(m_list,p_list,a,a0,n,b):

d_vec = []
for i in (0,1,2):
d_vec.append(m_der(m_list[i],p_list[(i+1)%3],a,a0,n))
for i in (0,1,2):
d_vec.append(p_der(m_list[i],p_list[i],b))
return d_vec

Note that the above code is still using quite a few python objects and methods, and could be sped up further with more effort.

Other notes

Sage includes scipy, numpy, and matplotlib for 2D plotting. For 3D plotting, Sage has Jmol and the tachyon raytracer.

Sage is also notable for its friendly community. In particular, if you get stuck you can usually get very fast support by writing to the sage-support list. There is also pretty good documentation available.

There is much more to Sage but hopefully this introduction has whet your appetite.