Delay differential equations in Python

I wrote a very simple and user-friendly method, that I called ddeint, to solve delay differential equations (DDEs) in Python, using the ODE solving capabilities of the Python package Scipy. As usual the code is available at the end of the post :).

Example 1 : Sine

I start with an example whose exact solution is known so that I can check that the algorithm works as expected. We consider the following equation:

y(t) = sin(t) \ \ \ for\ \ \ t < 0

y'(t) = y(t-3\pi/2) \ \ \ for\ \ \ t \geq 0

The trick here is that sin(t-3\pi/2) = cos(t)=sin'(t) so the exact solution of this equation is actually the sine function.

from pylab import *

model = lambda Y,t : Y(t - 3*pi/2) # Model
tt = linspace(0,50,10000) # Time start, end, and number of points/steps
g=sin # Expression of Y(t) before the integration interval 
yy = ddeint(model,g,tt) # Solving

# PLOTTING
fig,ax=subplots(1)
ax.plot(tt,yy,c='r',label="$y(t)$")
ax.plot(tt,sin(tt),c='b',label="$sin(t)$")
ax.set_ylim(ymax=2) # make room for the legend
ax.legend()
show()

dde_sin

The resulting plot compares our solution (red) with the exact solution (blue). See how our result eventually detaches itself from the actual solution as a consequence of many successive approximations ? As DDEs tend to create chaotic behaviors, you can expect the error to explode very fast. As I am no DDE expert, I would recommend checking for convergence in all cases, i.e. increasing the time resolution and see how it affects the result. Keep in mind that the past values of Y(t) are computed by interpolating the values of Y found at the previous integration points, so the more points you ask for, the more precise your result.

Example 2 : Delayed negative feedback

You can select the parameters of your model at integration time, like in Scipy’s ODE and odeint. As an example, imagine a product with degradation rate r, and whose production rate is negatively linked to the quantity of this same product at the time (t-d):

y(t) = 0 \ \ \ for\ \ \ t < 0

y'(t) = \dfrac{1}{1+(\dfrac{y(t-d)}{K})^2} -ry(t) \ \ \ for\ \ \ t \geq 0

We have three parameters that we can choose freely. For K = 0.1, d = 5, r = 1, we obtain oscillations !

from pylab import *

# MODEL, WITH UNKNOWN PARAMETERS
model = lambda Y,t,k,d,r :  1/(1+(Y(t-d)/k)**2) - r*Y(t)

# HISTORY
g = lambda t:0 

# SOLVING
tt = linspace(0,50,10000)
yy = ddeint(model,g,tt,fargs=( 0.1 , 5 , 1 )) # K = 0.1, d = 5, r = 1

# PLOTTING
fig,ax=subplots(1)
ax.plot(tt,yy,lw=2)
show()

dde_negativefeedback

Example 3 : Lotka-Volterra system with delay

The variable Y can be a vector, which means that you can solve DDE systems of several variables. Here is a version of the famous Lotka-Volterra two-variables system, where we introduce some delay d. For d=0 we find the solution of a classical Lotka-Volterra system, and for d non-nul, the system undergoes an important amplification:

\big(x(t), y(t)\big) = (1,2) \ \ \ for\ \ t < 0, \ \ else

x'(t) = 0.5x(t)\big(1-y(t-d)\big)\\   y'(t) = -0.5y(t)\big(1-x(t-d)\big)

from pylab import *

def model(Y,t,d):
    x,y = Y(t)
    xd,yd = Y(t-d)
    return array([0.5*x*(1-yd), -0.5*y*(1-xd)])

g = lambda t : array([1,2])
tt = linspace(2,30,20000)
fig,ax=subplots(1)

for d in [0, 0.2]:
    yy = ddeint(model,g,tt,fargs=(d,))
    # WE PLOT X AGAINST Y
    ax.plot(yy[:,0],yy[:,1],lw=2,label='delay = %.01f'%d)

ax.legend()
show()

dde_lotka

Example 4 : A DDE with varying delay

This time the delay depends on the value of Y(t) !

y(t) = 1,\ \ \ t \leq 0

y'(t) = - y\big(t-3\cos(y(t))^2 \big),\ \ \ t > 0

from pylab import *
model = lambda Y,t:  -Y(t-3*cos(Y(t))**2)
tt = linspace(0,30,2000)
yy = ddeint(model, lambda t:1, tt)
fig,ax=subplots(1)
ax.plot(tt,yy,lw=2)
show()

dde_Ydependent

Code

Explanations

The code is written on top of Scipy’s ‘ode’ (ordinary differential equation) class, which accepts differential equations under the form

model(Y,t) = ” expression of Y'(t) ”

where Y, and the output Y', must be Numpy arrays (i.e. vectors).

For our needs, we need the input Y to be a function of time, more precisely a function that can compute Y(t) at any past or present t using the values of Y already computed. We also need Y(t) to return the value of some function g(t) if t is inferior to some time tc that marks the start of the integration.

To this end, I first implemented a class (ddeVar) of variables/functions which can be called at any time t: for $t<latex t_c$, it will return the value of g(t), and for t>tc, it will look for two already computed values Y_a and Y_b at times t_a<t<t_b, from which it will deduce Y(t) using a linear interpolation. Scipy offers many other kinds of interpolation, but these will be slower and won't support vectors for Y.

Such variables need to be updated every time a new value of Y(t) is computed, so I created a class 'dde' that inherits from Scipy's 'ode' class but overwrites its integration method so that our special function Y is updated after each integration step. Since 'ode' would feed the model with a vector Y (a Numpy array to be precise), which we don't want, we give to the integrator an interface function that takes a Numpy array Y as an argument, but immediately dumps it and calls the model with our special ddeVar variable Y (I hope that was clear 🙂 ).

Ok, here you are for the code

You will find the code and all the examples as an IPython notebook HERE (if you are a scientific pythonist and you don’t know about the IPython notebook, you are really missing something !). Just change the extension to .ipynb to be able to open it. In case you just asked for the code:

# REQUIRES PACKAGES Numpy AND Scipy INSTALLED
import numpy as np
import scipy.integrate
import scipy.interpolate
 
class ddeVar:
    """ special function-like variables for the integration of DDEs """
    
    
    def __init__(self,g,tc=0):
        """ g(t) = expression of Y(t) for t<tc """
        
        self.g = g
        self.tc= tc
        # We must fill the interpolator with 2 points minimum
        self.itpr = scipy.interpolate.interp1d(
            np.array([tc-1,tc]), # X
            np.array([self.g(tc),self.g(tc)]).T, # Y
            kind='linear', bounds_error=False,
            fill_value = self.g(tc))
            
            
    def update(self,t,Y):
        """ Add one new (ti,yi) to the interpolator """
        
        self.itpr.x = np.hstack([self.itpr.x, [t]])
        Y2 = Y if (Y.size==1) else np.array([Y]).T
        self.itpr.y = np.hstack([self.itpr.y, Y2])
        self.itpr.fill_value = Y
        
        
    def __call__(self,t=0):
        """ Y(t) will return the instance's value at time t """
        
        return (self.g(t) if (t<=self.tc) else self.itpr(t))
 


class dde(scipy.integrate.ode):
    """ Overwrites a few functions of scipy.integrate.ode"""
    
    
    def __init__(self,f,jac=None):
		
        def f2(t,y,args):
            return f(self.Y,t,*args)
        scipy.integrate.ode.__init__(self,f2,jac)
        self.set_f_params(None)
        
 
    def integrate(self, t, step=0, relax=0):
		
        scipy.integrate.ode.integrate(self,t,step,relax)
        self.Y.update(self.t,self.y)
        return self.y
        
 
    def set_initial_value(self,Y):
		
        self.Y = Y #!!! Y will be modified during integration
        scipy.integrate.ode.set_initial_value(self, Y(Y.tc), Y.tc)
 
 
 
def ddeint(func,g,tt,fargs=None):
    """ similar to scipy.integrate.odeint. Solves the DDE system
        defined by func at the times tt with 'history function' g
        and potential additional arguments for the model, fargs
    """
    
    dde_ = dde(func)
    dde_.set_initial_value(ddeVar(g,tt[0]))
    dde_.set_f_params(fargs if fargs else [])
    return np.array([g(tt[0])]+[dde_.integrate(dde_.t + dt)
                                 for dt in np.diff(tt)])

Other implementations

If you need a faster or more reliable implementation, have a look at the packages pyDDE and pydelay, which seem both very serious but are less friendly in their syntax.

Typing keyboard + python = Musical instrument !

In this post I will show you how to do this :

I am not the first to do that, but most people who play their computer on Youtube use either very expensive programs, or programs that won’t run on your computer, or not-so-efficient programs with not-that-much possibilities of extension, or cheap programs with a big lag between pressing the key and actually hearing the note.

So here is a very small Python script which will run fine even on a basic netbook. If you are not familiar with Python, you should take online courses , it is really worth it 🙂 !
If you are faminiliar with Python, then you are welcome to improve the code on its Github page.

Transforming your keyboard into a mixtable

My original idea was to transform my computer keyboard into a mixtable, to make something like in this very awesome video:

So my first move was to make a program that would take a configuration file my_configuration.cf containing this:


q, dog.wav
w, cat.wav
e, tarzan.wav
r, scream.mp3

And then if you hit q you would hear a dog from the dog.wav soundfile, if you hit w you’d hear a cat, etc…
This is pretty easy to do with Python’s pygame package. Here is my code (inspired by similar stuff from the package Mingus):


import pygame as pg
import csv
import time


SAMPLE_WIDTH = 16
FPS = 44100
N_CHANNELS = 2
BUFFER = 2**9
    
def minimix(config_file,mode = 'instrument'):
    """
    Opens an interface that lets you press the keys of your keyboard
    to plays the related soundfiles as sepcified in the provided
    configuration file.
    
    Args:
       config_file (str):  A file associating keyboard keys with
                           file names. Each line must be of the form
                           key , filename.
                       
       mode (str) :
            instrument -- causes autorepeat of samples as long
                            as the key is pressed.
            sustain -- causes autorepeat of the samples until its
                       key is pressed again.
            player -- striking the key causes the sound to play once. 
    
    Returns:
       a list of the (time,events).
    """
    
    repeat = 0 if (mode is 'player') else (-1)
    
    pg.mixer.pre_init(FPS,-SAMPLE_WIDTH,N_CHANNELS,BUFFER)
    pg.init()
    screen = pg.display.set_mode((640,480))
    
    
    
    
    ##### READ CONF FILE
    
    key2sound = {}
    key2file = {}
    config = csv.reader(open(config_file, 'rb'), delimiter=',')
    
    for key, soundfile in config:
        
        key,soundfile = key.strip(' '),soundfile.strip(' ')
        
        if key is not '#':
            
            key2file[key] = soundfile
            key2sound[key] = pg.mixer.Sound(soundfile)
    
    events_list = []
    currently_playing = {k : False for k in key2sound.iterkeys()}
    
    
    ##### MAIN LOOP
    
    while True:
        
        event =  pg.event.wait()
      
        if event.type in (pg.KEYDOWN,pg.KEYUP):
            key = pg.key.name(event.key)
              
            if key in key2sound:
                
                if event.type == pg.KEYDOWN:
                    
                    if (mode == 'sustain') and currently_playing[key]:
                        
                        key2sound[key].fadeout(20)
                        currently_playing[key] = False
                        
                    else:
                        
                        key2sound[key].play(repeat) 
                        currently_playing[key] = True
                    
                    events_list.append((time.time(),key2file[key]))
                    
                elif event.type == pg.KEYUP and (mode == 'instrument'):
                    
                    key2sound[key].stop()
                    currently_playing[key] = False
                    
                    events_list.append((time.time(),key2file[key]))
            
            elif event.key == pg.K_ESCAPE:
                
                break
    
    pg.quit()
    
    return events_list

Transforming your keyboard into a musical instrument

If instead of using various noises like cat and dog you use different notes from the same instrument, then you turned your computer into some kind of piano. The problem is that a set of soundfiles with all the notes of an instrument is difficult to find on the internet, so I wrote a script that makes as many notes as you want from just one sound by shifting its pitch up or down. It uses the audio processing program Soundstretch that you will need to install first :

import os

def shift_wav(wavfile,output,shifts,verbose=False):
    """
    Makes new sounds by shifting the pitch of a sound.
    Requires soundstrech installed.
    
    Args:
        wavfile : Name of the file containing the original sound
        output: name to use as a prefix for the output files and for
                the output folder name
        shifts (list of int): specifies of how many half-tones the pitch
                shifts should be. For instance [-2,-1,1,2] will produce
                4 files containing the sound 2 half-tones lower, one
                halftone lower, one halftone higher and two halftones
                higher.
    """
    
    folder = os.path.dirname(output)
        
    if not os.path.exists(folder):
        
        os.makedirs(folder)
        
    for i,s in enumerate(shifts):
        
        outputfile = '%s%02d.wav'%(output,i)
        
        command = 'soundstretch %s %s -pitch=%d'%(wavfile,outputfile,s)
        if verbose:
            print command
        os.system(command)

Going further

There is so much one could do to improve the program.

On the musical side, for instance, finding configurations of the keyboard that are particularly ergonomic. In the video above I used this configuration:

azerty

I called it typewriter because it enables you to play very fast things while moving your hands a minimum (the video is a bad example :P) . But maybe there is better to find !

Also, one could start listing every cool piece of music that can be played on a typing keyboard. I use 46 keys ( almost 4 octaves !), that makes a lot of possibilities !

On the programming side, there is a lot of little things I can think of, like automatizing scale changes, introducing nuances, designing nice interfaces (why not a guitar-hero-like game where you would actually be playing music on a playback ?), writing a script that would take some sheet music (in a nice format, like ABC, MIDI, lylipond) and return the list of the keys you should strike to play it.

I actually wrote a lot more code, for instance to make it easier to write configuration files, for sound processing, etc., but since it is not strictly necessary I am not reporting it here (I’ll certainly put a working version on GitHub, or such, some day).

Philosophy of the musical keyboard

Your typing keyboard is a real instrument. Of course it is not its primary use, but our voice’s primary purpose was not to sing, either. Now do the math : how many people out there have a piano at home ? And how many have a computer ? That gives you an idea of how many people would like to play the piano, cannot, but could play their computers instead.

So promoting computer-keyboardism is ultimately about bringing music to the masses. It is about providing everyone with an instrument that you can practice at home, in the train, at work, and that will be familiar to anyone everywhere in the world.

There is more : how many of you, pianist readers, have started the piano for seduction purposes ? (yeah, sure, me neither…) But public places with a piano on which you could show off your mad skills are getting pretty rare, aren’t they ? Especially since most bars have traded their good old piano for a TV. But think about all these places with a computer at hand ! Yep, time for you to develop a talent that will be useful in real life !

So practice, get good, be one of the first composers for tomorrow’s instrument, impress your friends and spread the good news ! If you are still reading me after so much gibberish , then do not hesitate : you just proved how little you value your free time, you are the right person for the task !

Extract data from graph pictures with Python

If you want to transform a picture of a graph into exploitable data (which is very useful in science if you want to exploit a figure from an article without bothering the authors), here is a minimalistic interface written in python with the following features:

  • Data extraction from picture files or from a picture in the clipboard.
  • Data extraction from rotated graphs or graphs shown with (moderate) perspective.
  • Advanced interface (left-click to select a point, right-click to deselect).
  • Stores the points’ coordinates in a python variable and in the clipboard (for use in another application).

You can launch the interface with

points = pic2data()

This will either start a session using the picture from the clipboard, or , if there is none, wait for the clipboard to contain a picture. Alternatively you can use a picture from a file with

points = pic2data('graph.jpeg')

You will then be asked you to place the origin of the graph, as well as the coordinates of this origin (in case it it not (0,0)), and one reference point for each axis X and Y (i.e. points of these axis whose coordinates you know). Then you can select/deselect as many points of the curve as you want, and exit with the middle button.The list of selected points [(x1,y1),(x2,y2),…] is returned.

By default the program will consider that the graph is rectangular and parralel to the edges of the pictures (wich I will call straight in what follows). This will typically be the case for a graph from a scientific article. As a consequence the algorithm will automatically replace the reference point you chose for the X axis in order to put it at the same height as the origin, and it will replace the reference point for Y exactly above the origin. However if the graph on the picture is not straight, like in a photo, use the argument straight=False.

As an example, let us take a photo with a graph, like this one.

Fig. 1: Young Frederic Chopin disguised as Mozart.

As the graph is not straight we will use

points = pic2data('mozart.jpeg', straight = False)

Which gets you to that:

After placing the points and getting their coordinates one can redraw the plot with

from pylab import *
figure()
x,y = zip(*points)
plot(x,y,'o')
show()

And voilà !

Here is the code. Happy curving !

from urlparse import urlparse

import pygtk
import gtk
import tkSimpleDialog

import matplotlib.image as mpimg
import matplotlib.pyplot as plt

import numpy as np

def tellme(s):
    print s
    plt.title(s,fontsize=16)
    plt.draw()

def pic2data(source='clipboard',straight=True):
    """ GUI to get data from a XY graph image. Either provide the graph
        as a path to an image in 'source' or copy it to the clipboard.
    """
    
    ##### GET THE IMAGE
    
    clipboard = gtk.clipboard_get()
    
    if source=='clipboard':
        
        # This chunk tries the text content of the clipboard
        # and empties it if it is not a file path
        
        print "Waiting for an image in the clipboard..." 
        while not ( clipboard.wait_is_uris_available()
                    or clipboard.wait_is_image_available()):
            pass
            
        if clipboard.wait_is_uris_available(): # it's a path to a file !
            
             source = clipboard.wait_for_uris()[0]
             source = urlparse(source).path
             return pic2data(source)
        
        image = clipboard.wait_for_image().get_pixels_array()
        origin = 'upper'
    
    else: # source is a path to a file !
        
        image = mpimg.imread(source)
        origin = 'lower'

    ###### DISPLAY THE IMAGE
    
    plt.ion() # interactive mode !
    fig, ax = plt.subplots(1)
    imgplot = ax.imshow(image, origin=origin)
    fig.canvas.draw()
    plt.draw()
    
    ##### PROMPT THE AXES
    
    def promptPoint(text=None):
        
        if text is not None: tellme(text)
        return  np.array(plt.ginput(1,timeout=-1)[0])
    
    def askValue(text='',initialvalue=0.0):
        return tkSimpleDialog.askfloat(text, 'Value:',
                                         initialvalue=initialvalue)
    
    origin = promptPoint('Place the origin')
    origin_value = askValue('X origin',0),askValue('Y origin',0)
                                         
    Xref =  promptPoint('Place the X reference')
    Xref_value = askValue('X reference',1.0)
    
    Yref =  promptPoint('Place the Y reference')
    Yref_value = askValue('Y reference',1.0)
    
    if straight :
        
        Xref[1] = origin[1]
        Yref[0] = origin[0]
    
    ##### PROMPT THE POINTS
    
    selected_points = []
    
    tellme("Select your points !")
    print "Right-click or press 's' to select"
    print "Left-click or press 'del' to deselect"
    print "Middle-click or press 'Enter' to confirm"
    print "Note that the keyboard may not work."
    
    selected_points = plt.ginput(-1,timeout=-1)
    
    ##### RETURN THE POINTS COORDINATES
    
    #~ selected_points.sort() # sorts the points in increasing x order
    
    # compute the coordinates of the points in the user-defined system
    
    OXref = Xref - origin
    OYref = Yref - origin
    xScale =  (Xref_value - origin_value[0]) / np.linalg.norm(OXref)
    yScale =  (Yref_value - origin_value[1]) / np.linalg.norm(OYref)
    
    ux = OXref / np.linalg.norm(OXref)
    uy = OYref / np.linalg.norm(OYref)
    
    result = [(ux.dot(pt - origin) * xScale + origin_value[0],
               uy.dot(pt - origin) * yScale + origin_value[1])
               for pt in selected_points ]
    
    # copy the result to the clipboard
    
    
    clipboard.set_text('[' + '\n'.join([str(p) for p in result]) + ']')
    
    clipboard.store() # makes the data available to other applications
    
    plt.ioff()
    
    return result

Animate your 3D plots with Python’s Matplotlib

When you have a complicated 3D plot to show in a video or slideshow, it can be nice to animate it:

I obtained this surface with

import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import axes3d

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
X, Y, Z = axes3d.get_test_data(0.05)
s = ax.plot_surface(X, Y, Z, cmap=cm.jet)
plt.axis('off') # remove axes for visual appeal

To animate it I created the function rotanimate that you can use like that:

import numpy as np
angles = np.linspace(0,360,21)[:-1] # A list of 20 angles between 0 and 360

# create an animated gif (20ms between frames)
rotanimate(ax, angles,'movie.gif',delay=20) 

# create a movie with 10 frames per seconds and 'quality' 2000
rotanimate(ax, angles,'movie.mp4',fps=10,bitrate=2000)

# create an ogv movie
rotanimate(ax, angles, 'movie.ogv',fps=10) 

Here is the source-code:

import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import axes3d
import os, sys
import numpy as np


##### TO CREATE A SERIES OF PICTURES

def make_views(ax,angles,elevation=None, width=4, height = 3,
                prefix='tmprot_',**kwargs):
    """
    Makes jpeg pictures of the given 3d ax, with different angles.
    Args:
        ax (3D axis): te ax
        angles (list): the list of angles (in degree) under which to
                       take the picture.
        width,height (float): size, in inches, of the output images.
        prefix (str): prefix for the files created. 
    
    Returns: the list of files created (for later removal)
    """
    
    files = []
    ax.figure.set_size_inches(width,height)
    
    for i,angle in enumerate(angles):
    
        ax.view_init(elev = elevation, azim=angle)
        fname = '%s%03d.jpeg'%(prefix,i)
        ax.figure.savefig(fname)
        files.append(fname)
    
    return files



##### TO TRANSFORM THE SERIES OF PICTURE INTO AN ANIMATION

def make_movie(files,output, fps=10,bitrate=1800,**kwargs):
    """
    Uses mencoder, produces a .mp4/.ogv/... movie from a list of
    picture files.
    """
    
    output_name, output_ext = os.path.splitext(output)
    command = { '.mp4' : 'mencoder "mf://%s" -mf fps=%d -o %s.mp4 -ovc lavc\
                         -lavcopts vcodec=msmpeg4v2:vbitrate=%d'
                         %(",".join(files),fps,output_name,bitrate)}
                         
    command['.ogv'] = command['.mp4'] + '; ffmpeg -i %s.mp4 -r %d %s'%(output_name,fps,output)
    
    print command[output_ext]
    output_ext = os.path.splitext(output)[1]
    os.system(command[output_ext])



def make_gif(files,output,delay=100, repeat=True,**kwargs):
    """
    Uses imageMagick to produce an animated .gif from a list of
    picture files.
    """
    
    loop = -1 if repeat else 0
    os.system('convert -delay %d -loop %d %s %s'
              %(delay,loop," ".join(files),output))




def make_strip(files,output,**kwargs):
    """
    Uses imageMagick to produce a .jpeg strip from a list of
    picture files.
    """
    
    os.system('montage -tile 1x -geometry +0+0 %s %s'%(" ".join(files),output))
    
    
    
##### MAIN FUNCTION

def rotanimate(ax, angles, output, **kwargs):
    """
    Produces an animation (.mp4,.ogv,.gif,.jpeg,.png) from a 3D plot on
    a 3D ax
    
    Args:
        ax (3D axis): the ax containing the plot of interest
        angles (list): the list of angles (in degree) under which to
                       show the plot.
        output : name of the output file. The extension determines the
                 kind of animation used.
        **kwargs:
            - width : in inches
            - heigth: in inches
            - framerate : frames per second
            - delay : delay between frames in milliseconds
            - repeat : True or False (.gif only)
    """ 
        
    output_ext = os.path.splitext(output)[1]

    files = make_views(ax,angles, **kwargs)
    
    D = { '.mp4' : make_movie,
          '.ogv' : make_movie,
          '.gif': make_gif ,
          '.jpeg': make_strip,
          '.png':make_strip}
          
    D[output_ext](files,output,**kwargs)
    
    for f in files:
        os.remove(f)
    

##### EXAMPLE

if __name__ == '__main__':

    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    X, Y, Z = axes3d.get_test_data(0.05)
    s = ax.plot_surface(X, Y, Z, cmap=cm.jet)
    plt.axis('off') # remove axes for visual appeal
    
    angles = np.linspace(0,360,21)[:-1] # Take 20 angles between 0 and 360

    # create an animated gif (20ms between frames)
    rotanimate(ax, angles,'movie.gif',delay=20) 

    # create a movie with 10 frames per seconds and 'quality' 2000
    rotanimate(ax, angles,'movie.mp4',fps=10,bitrate=2000)

    # create an ogv movie
    rotanimate(ax, angles, 'movie.ogv',fps=10) 

Sound up your Python !

Long simulations and optimisations are a plea to the modern modelist. Worst of all are the programs of average length (say, 5 to 15 minutes) that do not let you the time to start anything else during the runs, and have you lose your day before you know it. I wrote a few Python functions to ease the pain. The first is called ring() and produces… a ring (like an egg timer ! but if you have office-mates you may use a more neutral sound like somenone coughing, or scratching, etc…). You can use it like this:

def myLongFunction():
        for i in range(100000):
            print "I'm the world's most useless function !"

myLongFunction();ring()

In case you want your function to always ring at the end, you can also use my decorator @ringatend when you define the function :

@ringatend
def myLongFunction():
        for i in range(100000):
            print "I'm the world's most useless function !"

myLongFunction() # will ring at the end, no need to call 'ring()'

Now, if you really want to make your colleagues jealous (or crazy) you can ask Python to play some music in the background while the function runs, with my decorator @inmusic

@inmusic
def myLongFunction():
   .................

myLongFunction() # will play music during the function's run !

Believe me, simulations look way shorter with on some Beny Hill (but they are actually slower 🙂 )

You can also combine the decorators to get a function that will play some music, and ring at the end:

@ringatend
@inmusic
def myLongFunction():

Here is the code for all that. You will need to install Pygame to play the sounds. I have also included some functions to splash popups: popup and @popupatend:

rom Tkinter import *
import tkMessageBox
import time
import pygame

# RING

def ring(sound='./cough.wav'):
    
    pygame.mixer.init()
    sound = pygame.mixer.Sound(sound=)
    sound.play()
    time.sleep(sound.get_length())
    pygame.mixer.quit()



def ringatend(target):
    """ a decorator to make your functions rings when they are done """
    
    def f(*args,**kwargs):
        
        result = target(*args,**kwargs)
        ring()
        return result
        
    return f



# IN MUSIC !

def inmusic(target):
    """ a decorator to play your functions in music """
            
    def f(*args,**kwargs):
        
        pygame.mixer.init()
        sound = pygame.mixer.Sound('./music.wav')
        sound.play()
        result = target(*args,**kwargs)
        sound.stop()
        pygame.mixer.quit()
        return result
            
    return f



# POPUPS

def popup(txt='Finished'):
    root = Tk()
    tkMessageBox.showinfo('Finished!',txt)
    root.destroy()
    root.mainloop()



def popupatend(target):
    
    
    def f(*args,**kwargs):
        
        # Run the function
        result = target(*args,**kwargs)
        
        # Make a pretty text
        fname = target.func_name
        pretty_args = ','.join([str(a) for a in args])
        pretty_kwargs = ','.join(['%s=%s'%(str(k),str(v))
                       for k,v in kwargs.iteritems()])  
        descr ='%s(%s,%s)'%(fname,pretty_args, pretty_kwargs)
        descr = descr[:100] # in case it really is too long
        
        # popup description
        popup(txt='%s done !'%(descr))
    
        return result
    
    return f



### EXAMPLES

@ringatend
@inmusic
def myLongFunction():
        for i in range(100000):
            print "I'm the world's most useless function !"

myLongFunction()

@popupatend
def myLongFunction():
        for i in range(100000):
            print "I'm the most useless function !"

myLongFunction()

If you want to go further with the decorators, here are versions of the decorators with options that let you choose the sound or music you want to be played.



def ringatend2(soundFile):
    """ a decorator generator using a specified soundfile """
        
    def decorator(target):
        """ a decorator using the specified soundfile """
        
        def f(*args,**kwargs):
            
            result = target(*args,**kwargs)
            ring(soundFile)
            return result
            
        return f
    
    return decorator

def inmusic2(soundFile):
    """ a decorator generator using a specified soundfile """
        
    def decorator(target):
        """ a decorator to play your functions in music """
                
        def f(*args,**kwargs):
            
            pygame.mixer.init()
            sound = pygame.mixer.Sound(soundFile)
            sound.play()
            result = target(*args,**kwargs)
            sound.stop()
            pygame.mixer.quit()
            return result
                
        return f
    
    return decorator


### EXAMPLES

@ringatend2('./tarzan.wav')
@inmusic2('./diesIrae.wav')
def myLongFunction():
        for i in range(10000):
            print "The world's most useless function !"

A GUI for the exploration of functions with Python / Matplotlib

Sliders can be a great tools for lazy engineers and modelists who want to play around with the parameters of a program without recompiling it a thousand times.
Most scientific languages ( python-pylab, matlab, scilab…) support some basic sliders methods but it can be quite long to create a GUI. In this blog I present a light, user-friendly function I implemented for my everyday work.

Example with a numerical function

def volume(x,y,z):
	""" Volume of a box with width x, heigth y, and depth z """
	return x*y*z

intervals = [ { 'label' :  'width',  'valmin': 1 , 'valmax': 5 },
              { 'label' :  'height',  'valmin': 1 , 'valmax': 5 },
              { 'label' :  'depth',  'valmin': 1 , 'valmax': 5 } ] 

inputExplorer(volume,intervals)

The new value of volume is automatically computed each time you move the sliders. It is also computed if you press Enter (which can be useful to get several values, when your function is non-deterministic). In case you are studying a function that is long to evaluate, you can decide that it will only be evaluated when Enter is pressed, by adding “wait_for_enter = True” to the arguments list.

Example with a graphical function

You can also decide to use a function that doesn’t return any value, but updates a plot. Let us illustrate this with the Lotka-Volterra model, in which wolves eat rabbits and die, leading to periodical fluctuations of both populations. From wikipedia:

\frac{dx}{dt} = x(\alpha - \beta y)
\frac{dy}{dt} = - y(\gamma - \delta  x)

This can be simulated easily in python with scipy’s function odeint :

import numpy as np
from scipy.integrate import odeint

# model

def model(state,t, a,b,c,d):
	x,y = state
	return [ x*(a-b*y) , -y*(c - d*x) ] # dx/dt, dy/dt

# vector of times

t_vec = np.linspace(0,10,500) # 500 time points between 0 and 10

# initial conditions

x0 = 0.5
y0 = 1

# parameters

a = 1
b = 2
c = 1
d = 3

# simulate and return the values at each time t in t_vec

result = odeint(model, [x0,y0], ts, args = (a,b,c,d) ) 

# plot

import matplotlib.pyplot as plt
plt.plot(ts,result)
plt.show()

However, trying by hand several values of x0, y0, a, b, c, d, can be fastidious. So let us wrap the solving and the plotting in a function and feed it to inputExplorer. Let me rewrite the whole code from scratch :

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint

def model(state,t, a,b,c,d):
	x,y = state
	return [ x*(a-b*y) , -y*(c - d*x) ]

ts = np.linspace(0,10,500)

fig,ax = plt.subplots(1)

def plotDynamics(x0,y0,a,b,c,d):
	ax.clear()
	ax.plot(ts, odeint(model, [x0,y0], ts, args = (a,b,c,d)) )
	fig.canvas.draw()

sliders = [ { 'label' :  label,  'valmin': 1 , 'valmax': 5 }
		 for label in [ 'x0','y0','a','b','c','d' ] ]

inputExplorer(plotDynamics,sliders)

Code

Enjoy !

import matplotlib.pyplot as plt
from matplotlib.widgets import Slider, Button

def inputExplorer(f, sliders_properties, wait_for_validation = False):
    """ A light GUI to manually explore and tune the outputs of 
        a function.
        slider_properties is a list of dicts (arguments for Slider )
        whose keys are in ( label, valmin, valmax, valinit=0.5, 
        valfmt='%1.2f', closedmin=True, closedmax=True, slidermin=None,
        slidermax=None, dragging=True)
        
        def volume(x,y,z):
            return x*y*z
    
        intervals = [ { 'label' :  'width',  'valmin': 1 , 'valmax': 5 },
                  { 'label' :  'height',  'valmin': 1 , 'valmax': 5 },
                  { 'label' :  'depth',  'valmin': 1 , 'valmax': 5 } ]
        inputExplorer(volume,intervals)
    """
        
    nVars = len(sliders_properties)
    slider_width = 1.0/nVars
    print slider_width
    
    # CREATE THE CANVAS
    
    figure,ax = plt.subplots(1)
    figure.canvas.set_window_title( "Inputs for '%s'"%(f.func_name) )
    
    # choose an appropriate height
    
    width,height = figure.get_size_inches()
    height = min(0.5*nVars,8)
    figure.set_size_inches(width,height,forward = True)
    
    
    # hide the axis
    ax.set_frame_on(False)
    ax.get_xaxis().set_visible(False)
    ax.get_yaxis().set_visible(False)
    

    # CREATE THE SLIDERS
    
    sliders = []
    
    for i, properties in enumerate(sliders_properties):
        
        ax = plt.axes([0.1 , 0.95-0.9*(i+1)*slider_width,
                       0.8 , 0.8* slider_width])
        sliders.append( Slider(ax=ax, **properties) )
    
    
    # CREATE THE CALLBACK FUNCTIONS
    
    def on_changed(event) : 
        
        res = f(*(s.val for s in sliders))
        
        if res is not None:
            
            print res
    
    def on_key_press(event):
        
        if event.key is 'enter':
            
            on_changed(event)   
    
    figure.canvas.mpl_connect('key_press_event', on_key_press)
    
    # AUTOMATIC UPDATE ?
    
    if not wait_for_validation:
        
        for s in sliders :
            
            s.on_changed(on_changed)
    
    
    # DISPLAY THE SLIDERS
    
    plt.show()

E. coli’s regulation network from RegulonDB to Python

RegulonDB provides a continuously updated network of interactions of the bacterium E. coli. In this post I will show how easy it is to explore this network using Python with its package NetworkX. In case you are in a hurry, I put the whole script with enough comments at the end of this post.

Getting the Graph

First we create an empty graph:

import networkx as nx
eColiNetwork = nx.DiGraph(name = 'E. Coli Network')

Then we fetch the RegulonDB file containing the informations on the network. To get an internet file in Python we use urllib2:

from urllib2 import urlopen
url = "http://regulondb.ccg.unam.mx/data/network_tf_gene.txt"
regulonFile = urlopen(url)

Now we will read this file line by line to fill our graph with nodes (genes) and edges (interactions). The file contains uninformative lines (comments begining with ‘#’ or empty lines) and lines of the form ‘gene1 gene2 sign …’ where the sign is ‘+’ or ‘-‘ according to whether the interaction is an activation or an inhibition.

for line in regulonFile :

	# ignore comments and empty lines

	if not line.startswith(('\n','\t','#')):

		# A line is of the form "gene1	gene2 sign ..."

		g1,g2, sign = line.split('\t')[:3]

		# standardize gene names ( all characters to lower case )

		g1,g2 = g1.lower(),g2.lower()

		# add the edge to the network

		eColiNetwork.add_edge(g1,g2)
		eColiNetwork[g1][g2]['sign'] = sign

Easy as pie with Python ! You can then save this graph for offline use :

nx.write_gpickle(eColiNetwork,'EcoliNetwork.gn') # To save
eColiNetwork = nx.read_gpickle('EcoliNetwork.gn') # To load

Study and Statistics

It is very easy to get informations on a NetworkX graph. Let us start with some basic statistics:

nx.info(eColiNetwork)
>>> Name: E. Coli Network
    Type: DiGraph
    Number of nodes: 1696
    Number of edges: 3809
    Average in degree:   2.2459
    Average out degree:   2.2459

Note that there are only 1700 genes out of the circa 4000 genes of E. coli. Moreover it doesn’t have so much edges, so we can expect it not to be fully connected, i.e. maybe you cannot always go from one gene to another in this graph. Let us check that out.

# undirected version of E. coli's network
undirected_net = eColiNetwork.to_undirected()
 
connected_components = nx.connected_components(undirected_net)

print 'Number of connected components : ', \
		len(connected_components )

print 'Length of the connected components : ', \
		[ len(c) for c in connected_components]

 

>>> Number of connected components :  22
>>> Length of the connected components :  [1603, 13, 7, 6, 6, 6, 6, 5, 5, 5, 4, 4, 4, 3, 3, 3, 3, 3, 2, 2, 2, 1]

Indeed we have 21 small satellite networks and one big central component. Let’s learn more about the latter:

main_component = nx.connected_component_subgraphs(undirected_net)[0]

print "Average shortest path in E. coli's main connected component : ", \
	   nx.average_shortest_path_length(main_component)

print "Diameter of E. coli's main connected component : ", \
	   nx.diameter(main_component)

 

>>>  Average shortest path in E. coli's main connected component :  3.61111539459
>>>  Diameter of E. coli's main connected component :  9

We can also get information on a specific gene, like its regulators (which are its predecessors in the graph, i.e. the genes pointing at it) and the genes that it regulates (its successors). Let’s have a look at acs and allr :

acs_regulators = eColiNetwork.predecessors('acs')
print "Genes regulating acs : ", " ".join( acs_regulators)

allr_targets = eColiNetwork.successors('allr')
print "Genes regulated by allr : ", " ".join(allr_targets)
>>> Genes regulating acs :  fis crp ihf
>>> Genes regulated by allr :  allb alla hyi gcl glxr alls ybby ybbw glxk

Sketching a Subgraph

Because it is so easy, let us draw some networks with NetworkX. Displaying the whole graph, or even its main component, would be ugly, as they are too big, so we choose the second connected component:

component = nx.connected_component_subgraphs(undirected_net)[1]

We use pylab for its plotting capabilities :

from pylab import *

figure()
nx.draw(component, node_size = 1000)
show()

And voilà !

E. coli subnetwork

E. coli subnetwork

Complete Code

Have fun 🙂


# (F)utilities to explore the network of Escherichia coli
# as given on RegulonDB


##### FETCH THE NETWORK ON THE INTERNET ################################


import networkx as nx
from urllib2 import urlopen

# Fetch E. coli's network on Regulon DB as a .txt file

url = "http://regulondb.ccg.unam.mx/data/network_tf_gene.txt"
regulonFile = urlopen(url)

# Create an (empty) directed graph

eColiNetwork = nx.DiGraph(name = 'E. Coli Network')

# Read the file and fill our network

for line in regulonFile :

	# ignore comments and empty lines

	if not line.startswith(('\n','\t','#')):

		# A line is of the form "gene1	gene2	+ "

		g1,g2, sign = line.split('\t')[:3]

		# standardize gene names ( all characters to lower case )

		g1,g2 = g1.lower(),g2.lower()

		# add the edge to the network

		eColiNetwork.add_edge(g1,g2)
		eColiNetwork[g1][g2]['sign'] = sign


########### SAVE AND LOAD ##############################################


import pickle

# Save the graph as EcoliNetwork.gn

with open('EcoliNetwork.gn', 'wb') as output:
	pickle.dump(eColiNetwork, output, pickle.HIGHEST_PROTOCOL)

# Load the graph

eColiNetwork = pickle.load(open('EcoliNetwork.gn', 'r'))


#### PRINT INFOS ON THE NETWORK ########################################


# undirected version of E. coli's network

undirected_net = eColiNetwork.to_undirected()

# Stats

print '-- General informations : \n', \
		nx.info(eColiNetwork)

connected_components = nx.connected_components(undirected_net)

print 'Number of connected components : ', \
		len(connected_components )

print 'Length of the connected components : ', \
		[ len(c) for c in connected_components]

main_component = nx.connected_component_subgraphs(undirected_net)[0]

print "-- Average shortest path in E. coli's main connected component : ", \
	   nx.average_shortest_path_length(main_component)

print "-- Diameter of E. coli's main connected component : ", \
	   nx.diameter(main_component)

# Infos on a specific gene :

print "-- Genes regulating acs : ", \
	   " ".join( eColiNetwork.predecessors('acs') )

print "-- Genes regulated by allr : ", \
	   " ".join( eColiNetwork.successors('allr') )


###########  D R A W  ##################################################


from pylab import *

figure()
component = nx.connected_component_subgraphs(undirected_net)[1]
nx.draw(component, node_size = 1000)
show()