http://cnr.lwlss.net

Guerilla knitting in Dublin.

Ageing has been solved. We are all out of a job.

Our paper describing the use of a simple, well-parameterised cell population model to compare the usefulness of various cellular markers as indicators of cellular senescence is now available in Experimental Gerontology. We used CaliBayes (published here) to give distributed estimates of the increase in senescent fraction with time in passaged human fibroblasts.

DOI: 10.1016/j.exger.2010.01.018.

A few weeks ago, I had the very great pleasure of visiting the home of Prof. David Jones to take some photographs of him and some of his creations as promotional material for a documentary about him by Adrin Neatrour. Prof. Jones is probably most famous for writing the Daedalus columns in New Scientist and Nature, but is also well known for building unridable bicycles, investigating Napoleon's death, and building puzzlingly perpetual motion machines.

He is a fascinating, mischevious and inspirational scientist and engineer with a house packed to the brim with amazing tools and gadgets. Prof. Jones showed me a lathe where he had recently carved a set of billiard balls with centres of gravity which were not located at their centres to give to a relative as a christmas present!

Adrin's documentary is really enjoyable, I recommend it to anyone, and it will shown, for free at the Tyneside Cinema at 3pm on Tuesday 16th March 2010 as part of the Newcastle Science Festival 2010.

The documentary is distributed by Galloping Films.

I stumbled across the Python winsound module today. It is a standard library on Windows machines and it has a function for controlling the PC speaker. Harf.

Hear it play the music from Close Encounters of a Third Kind from a list of frequencies and durations:

import winsound
from winsound import *
music=((784,1),(880,1),(698,1),(349,1),(523,3))
for x in music:
    Beep(x[0],550*x[1])

A list of frequencies and durations is a simple, computational form of musical notation.

Our new article published today in Briefings in Bioinformatics, describing CaliBayes (Bayesian inference for parameter distributions in SBML models, including discrete stochastic models) and BASIS (a tool for simulating cohorts of stochastic simulations from an SBML model and storing the results). Both of these systems are accessible as webservices and there are various convenient Python, Java and R client libraries to make consuming these webservices extremely straightforward. Both systems make significant computing power at Newcastle University publicly available.

Free access to .pdf

The DOI entry for the article is now working:

10.1093/bib/bbp072

After reading Darren's post about the Mandelbrot set I was browsing through the Wikipedia entry on the Mandelbrot set which has a few amazing animations of really deep zooms on the fractal bounded by the set.

The good/bad thing about fractals is that as you magnify them, the pattern and resulting image that you see is in some way similar to patterns you saw at lower magnifications. So the zooms are in some sense self-similar.

It occurred to me that viewing random walks at different magnifications should have a similar property. You should not be able to tell what magnification you are looking at, since the average trajectory should be random, regardless of the scale of the trajectory. As long as you are away from the scale at which you can observe any discrete steps making up the walk, that is.

So, we should be able to make semi-infinite zooms of random walks, in which all magnifications look believably random and as if they were generated from the same process but without any repeating or recurring patterns. I think the random walks are quite pleasing to look at. It is a bit tricky to create this zoom though, since I can't think of any cunning way to step from one magnification to another by "filling in" the gaps in the walk that have just come into view. So I have to create a massive random walk image, and progressively resize and crop it to achieve the zoom. Also, the PyGame code I posted a while ago is too slow for generating a massive random walk. However, there's no need to see the walk being generated in this case, so we can just colour pixels black and anti-alias during image resize (using the Python Imaging Library) instead of doing it in real-time.

Deep zooming is cool anyway as you can see in the amazing Powers of Ten video by Charles and Ray Eames.

So, here's a zoom on a random walk that I made (approximately 100 times zoom, made on a 32-bit machine with 3.5Gb RAM):

The frames were generated with the Python code below. If you have a more powerful computer than me to hand, you can generate a deeper zoom! Just increase the width and height parameters until you run out of memory or patience. I patched the output images together and created a .gif animation using the amazing VirtualDub.

import PIL,random,math,os
from PIL import Image,ImageDraw
 
# Define screen size and initialize
size=width,height=3500,3500
plot=Image.new("L",(width,height),'white')
pix=plot.load()
 
# Output Frame size
wf,hf=250,250
# Zoom parameters
zmax,zno=1.0,700
 
# Start at centre of image
x,y=width/2,height/2
pix[x,y]=0
 
# Prepare to find the bounding box for the walk
xmin,ymin,xmax,ymax=x,y,x,y
pcount=0
 
# Start walking until we hit the edge
while 1<=x<width-1 and 1<=y<height-1:
        # Generate new steps
        xnew=x+random.randint(-1,1)
        ynew=y+random.randint(-1,1)
        x,y=xnew,ynew
        pix[x,y]=0
        pcount+=1
        if x<xmin:xmin=x
        if x>xmax:xmax=x
        if y<ymin:ymin=y
        if y>ymax:ymax=y
 
# Crop to the walk      
plot=plot.crop((xmin,ymin,xmax,ymax))
 
# Rotate if necessary to make image wider than tall
if xmax-xmin<ymax-ymin:
        plot=plot.rotate(90)
(w,h)=plot.size
 
# Set minium zoom to just capture entire walk
zmin=float(wf)/float(w)
 
# For geometric zooming:
# zoom(z) = a*pow(phi,b*z)
# insist that zoom(0)=zmin, zoom(zno)=zmax
phi=10.0
a=zmin
b=math.log(zmax/zmin,phi)/zno
 
for z in xrange(0,zno+1):
        zoom=a*pow(phi,b*z)
        print "FRAME: ",z, zoom
        blank=Image.new("L",(wf,hf),'white')
        if z>0 and z%50==0:
                # Throw away pixels to achieve deeper zoom with set memory
                cz=float(hf)/float(hnew)
                xa=max(0,int(round(float(w)/2.0-cz*float(w)/2.0)))
                ya=max(0,int(round(float(h)/2.0-cz*float(h)/2.0)))
                xb=min(w,int(round(float(w)/2.0+cz*float(w)/2.0)))
                yb=min(h,int(round(float(h)/2.0+cz*float(h)/2.0)))
                plot=plot.crop((xa,ya,xb,yb))
                (w,h)=plot.size
        wnew=int(round(zoom*float(w)))
        hnew=int(round(zoom*float(h)))
        # Resize the frame
        small=plot.resize((wnew+1,hnew+1),Image.ANTIALIAS)
        # This crop is to trim grey lines which appear
        small=small.crop((1,1,wnew+1,hnew+1))
        xmin=int(round(float(wnew)/2.0-float(wf)/2.0))
        ymin=int(round(float(hnew)/2.0-float(hf)/2.0))
        xmax=int(round(float(wnew)/2.0+float(wf)/2.0))
        ymax=int(round(float(hnew)/2.0+float(hf)/2.0))
        if wnew<wf or hnew<hf:
                frame=Image.new("L",(wf,hf),'white')
                px=int(round(float(wf-wnew)/2.0))
                py=int(round(float(hf-hnew)/2.0))
                frame.paste(small,(px,py))
        else:
                frame=small.crop((xmin,ymin,xmax,ymax))
        fname="FrameL%05d.png"%z
        frame.save(fname)
 
print "Zoom: ",1.0/zmin
print "Steps: ",pcount

I went to a nice talk today by Nicolas Le Novere about some current developments in SBML, a standard language for describing biological models to computers.

During his talk he mentioned some of his work on modelling neurons and flashed up an image representing a molecule diffusing into a synapse. I think it was undergoing Brownian motion, which is basically a random walk in space (looked like 2-dimensional space to me).

In May I wrote some Python code using PyGame which simulates random walks, generating still images of them, but also displaying an animation of their generation (see below). It generates anti-aliased lines, so it still looks smooth even when taking tiny step sizes. Simply run it, and it will show you how it is choosing the random steps in the walk, and finally it will produce a .png image which looks like this:

You can change the resolution, step-sizes etc. yourself. I am not sure if the steps should be normally distributed or uniformly distributed to be true Brownian motion, but I can't see any difference in the patterns generated by eye. Anyway, it's a nice simple example of using PyGame for something non-game related.

import pygame, sys, numpy
from pygame.locals import *

# Define screen size
size=width,height=1000,1000
# It's much faster if we don't show the walk developing...
showwalk=True

# Initialise screen, set size
pygame.init()
screen = pygame.display.set_mode(size)
pygame.display.set_caption('A Random Walk')

# Fill background with white
background = pygame.Surface(size)
background = background.convert()
background.fill((255, 255, 255))

# Start at centre of image
x,y=width/2,height/2

def Exit():
        # Blit everything to the screen
        screen.blit(background, (0, 0))
        pygame.display.flip()
        # Save image
        fname="RandWalk"+str(numpy.random.randint(0,99999))+".png"
        pygame.image.save(background,fname)
        # Quit gracefully
        pygame.quit()
        sys.exit(0)

# Start walking until we hit the edge
while 0<=x<width and 0<=y<height:
        # Generate new steps
        xnew=x+numpy.random.uniform(-1.0,1.0)
        ynew=y+numpy.random.uniform(-1.0,1.0)
        # Let's try to draw an anti-aliased line
        pygame.draw.aaline(background, (0,0,0), (x,y), (xnew,ynew), True)
        x,y=xnew,ynew
        # Blit everything to the screen if we want to see development
        if showwalk:
                screen.blit(background, (0, 0))
                pygame.display.flip()
        # Exit gracefully if required by pressing q
        for event in pygame.event.get():
                if event.type == KEYDOWN and event.key==K_q:
                        Exit()
Exit()

Imagine you want a list of coordinates, following a spiral, on a discrete, regular grid. I can think of several examples where this could be useful for digital image construction, where the regular grid corresponds to pixel coordinates (integer only). It's easy to code up a square spiral like Ulam's spiral, since pixels are naturally arrayed on a rectangular grid, but using square spirals in image construction leaves obvious square corners. A circular spiral can be smoother and better.

Let's look at the Archimedes spiral in polar coordinates:

If we choose a = 0 and select b so that:

we can write the equations describing the spiral in Cartesian coordinates, in terms of , the distance between spiral loops:

To get an approximation of this smooth curve on a regular grid, we can think about trying to move along it in jumps of constant length , where the jump size is achieved by selecting a change in angle . Computing the distance between one point and the next we can write down this expression:

We can rearrange this expression into something simpler:

It would be nice to solve this for , then we would know how big an angle we need to jump to move exactly at a given time. Remember that as the radius gets bigger, the angle required to achieve a given displacement decreases. Unfortunately, I can't solve it analytically, due to the presence of the cosine term, however we can approximate cos() with the Taylor Series expansion up to 3rd order, giving us this expression:

This can easily be solved for small values of . If we define the perpendicular distance between pixel centres as 1, we know that the maximum distance between any two adjacent pixels is 1/. Therefore if we choose our parameters like this:

we would expect to move along the spiral, moving close to every pixel at least once. If we continue to step along the spiral like this indefinitely, and round every coordinate to the nearest integer value, we should hit every pixel in Archimedes spiral order!

This diagram shows the result of applying this algorithm:

The blue curve is a continuous Archimedes spiral, the numbered red points are the pixel approximations of the continuous curve and the red lines demonstrate the actual path taken between unique, discrete approximation points.

Below is a demonstration of what this sort of thing can look like. In this case I sorted the pixels of an original image by intensity and re-drew them along the approximate pixellated Archimedes spiral path:

The original image is here:

I've written some python code to generate these rasterized pixel spirals. You specify an image size (rectangle) and the coordinates of the spiral focus, and it generates a unique list of points lying inside that rectangle, in spiral order (as in the diagram above). It's quite fast, it can generate enough points to fill a 6 megapixel image in a couple of minutes.

from math import *
 
def uniquer(seq, idfun=None): # http://code.activestate.com/recipes/52560/
    # Generates a unique list of elements from seq
    # maintaining their original order
    if idfun is None:
        def idfun(x): return x
    seen = {}
    result = []
    for item in seq:
        marker = idfun(item)
        if marker in seen: continue
        seen[marker] = 1
        result.append(item)
    return result
 
def x(theta,alpha):
    return int(round(alpha*theta*cos(theta)/(2*pi)))
 
def y(theta,alpha):
    return int(round(alpha*theta*sin(theta)/(2*pi)))
 
def delta(alpha, theta, omega):
    return (2*pi*omega)/sqrt(pow(alpha,2)*(1+pow(theta,2)))
 
def dist(a,b):
    return sqrt(pow(a[0]-b[0],2)+pow(a[1]-b[1],2))
 
def rasterspiral(TL,BR,F):
    # Distance between spiral loops
    alphas=1/sqrt(2)
    # Distance along spiral at each step
    omegas=1/sqrt(2)
    # Starting theta
    thetas=0
 
    # Define top left pixel and bottom right pixel coordinates
    (tlx,tly)=TL
    (brx,bry)=BR
    # Spiral focus coordinates (must lie within rectangle above)
    (xf,yf)=F
 
    # Find which corner is farthest from the spiral focus
    corners=[[tlx,tly],[tlx,bry],[brx,bry],[brx,tlx]]
    distances=[dist([xf,yf],i) for i in corners]
    farthestind=max(xrange(len(distances)), key=distances.__getitem__)
    farthest=corners[farthestind]
 
    # Initialise the spiral at F
    xc,yc=xf,yf
    results=[(xc,yc)]
 
    # Keep spiralling until we've capture the farthest corner
    while [xc,yc]!=farthest:
        # Update the step in angle based on current angle etc.
        deltanew=delta(alphas,thetas,omegas)
        thetas=thetas+deltanew
        # Calculate the next coordinates
        xc,yc=xf+x(thetas,alphas),yf+y(thetas,alphas)
        # If these coordinates lie within the defined rectangle
        # then store them
        if tlx<=xc<=brx and tly<=yc<=bry:
            results.append((xc,yc))
 
    # Now eliminate duplicate copies of points on the spiral
    results=uniquer(results)
    return results
 
# Example of use:
test=rasterspiral((0,0),(50,100),(25,50))

bikehike is a UK site on which you can record run/hike/cycle routes online. Routes can be input by either clicking on a google map (or indeed using google routefinder) or by uploading a set of GPS coordinates as generated by a GPS-aware training watch for example. This is pretty much how the gmap-pedometer site works, BUT this is far more awesome since as well as using the full range of google map tools (satellite view, nice interface for zooming around, postcode search etc.) it also incorporates information from the ordnance survey, such as elevation data! Also OS maps tend to have many more nameplaces on them than google maps.

To see it in action, here is the route I run on Wednesdays. If you check the Show Elevation Data box, you can see the route elevation, and if you click on the Satellite button, you can see a browsable overhead photograph of the route. This is far cooler than the gmap-pedometer routes I had saved previously.

Now all I need is one of those GPS watches...