Interactive coalescents

Genetic coalescents are interesting statistically; the variance in the time to coalescence is large, which is the kind of quantity I think human intuition has trouble with. So it helps a bit to be able to play with them (code can be found here):

Gfan in 3D

Having upgraded the Sage interface to gfan for version 0.3, I've been thinking about other ways to leverage Sage's capabilities in this respect. One thing I've been working on is a 3D Groebner fan representation. I have some working code for this now, which hopefully will end up in Sage-3.0 if I have time to polish it up. Below are a couple of screenshots of the 3D rendering of the Groebner fan of the ideal generated by (w^3-x, x^3-y, y^3-w, z^2-x-y-w):

restricted four-body problem

The Albouy-Chenciner equations for the restricted four-body problem seem to generally have 191 solutions. For masses m1=17/20, m2 = 19/20, and m3 = 1, there are 160 complex solutions in 6 variables (the mutual distances between the particles). Plotting each sextuplet as a polygon gives the following plot, followed by the configurations formed by the positive real solutions. Computed with Sage and phcpack.

phcpack, sage, and interact

I've been working on integrating Jan Verschelde's phcpack software with Sage. phcpack finds solutions to pretty nasty systems of multivariable polynomials by using polyhedral homotopy continuation. Sage can provide a nice frontend for this. Here's an interactive display of the complex solutions of the Albouy-Chenciner equations (from the paper "Le problème des n corps et les distances mutuelles", Inventiones Mathematicae 131 p.151-184, 1998) for the 4-vortex problem:

Color me Gfan

The latest version of Gfan has some new capabilities that I am excited to use for testing whether ideals are zero-dimensional. But first I have to rewrite the Sage interface to Gfan. I thought that I should try to give some Sage-added-value while I was at it, so I am converting Gfan's xfig output to Sage graphics and adding some color. Here's one result so far: a map of all the reduced Groebner bases for the 3-vortex problem, colored by maximum degree:

Sage 2.10.3

Sage 2.10.3 is out, with the first released version of the new interact command. As a spectator to the process, it looked like a tough fight between bugs and developers - which I think should be viewed as an entirely positive thing, since it is a consequence of improved QA practices.

Release notes are here, in case you need it the download page is here.

If I can stop playing with interact (which has already been useful to me in teaching and research after 1 day!!) I hope to contribute a wee bit to some upcoming releases. I am rewriting the gfan interface to make use of gfan 0.3. Also, a student and I are working on pretty graphics for path-tracking solutions from homotopy solvers of polynomial systems (phcpack). I'll have to revise the plan a bit to exploit interact.

Interact

William Stein and Co. have delivered once again, with a new "interact" command that looks amazing even in its beta form. You can almost smell Sage-3.0; it should be out before the roses are blooming here in Duluth. Among other tests, I used it to explore the CpG content of the human mitochondrion averaged over a variable length window:

Music

Okay, so its not about bioinformatics or Sage, but after a month of no posts I thought I should say something. In the near future I hope to put up some more technical stuff; right now I am working on updating the Gfan (Groebner basis fan software) component of Sage, which isn't very sexy.

So: some music picks. I have been very happy with Balkan Beat Box's eponymous first album. I highly recommend it. Also good recently was Dengue Fever's Escape From Dragon House.

coalescent

Given the title of this blog, I've meant to do some coalescent stuff in sage for a while. Now that I am about to teach about it in a class, I've finally found the motivation:
(It's cooler animated, but I can't figure out how to post a gif animation here.)

Sage in 3D

It appears that Jmol will be the backend for most 3D graphics in Sage. After some hard work by developers (that I had nothing to do with) some of that functionality was put in Sage in time for the joint math meetings in San Diego. William somehow managed to find a slew of 3D glasses so that folks could enjoy the stereoscopic view:

(photo by Robert Bradshaw)

It will be exciting to see what happens as the Jmol integration proceeds.

G+C content as a linked image map

An image map of the human mitochondrion, linked to its GenBank entry by region, showing the G+C content in a small window. Created in Sage using the biopython optional module and this source code.

three body problem central configurations

All 99 (complex) solutions of the Albouy-Chenciner equations for central configurations of the Newtonian three-body problem, with equal masses. There are three "distances" for each solution; each set of distances are drawn as a triangle. The unit circle is included for a sense of scale. The solutions of primary interest are real, but understanding the complex solutions is important for analyzing these equations.

This was computed with Sage, with the bulk of the work was done with Jan Verschelde's PHCpack. At the moment, PHCpack is in sage as an optional package, so you would have to install that to get the code to work: source code.

spheres

Light reflecting in four tangent spheres....this relates to a math problem I abandoned a few years ago, maybe I'll work on it again someday. The sage code for the above:

t = Tachyon(camera_center=(0,-4,1), xres = 1200, yres = 800, raydepth = 12, aspectratio=.75, antialiasing = True)
t.light((0.02,0.012,0.001), 0.01, (1,0,0))
t.light((0,0,10), 0.01, (0,0,1))
t.texture('s', color = (.8,1,1), opacity = .9, specular = .95, diffuse = .3, ambient = 0.05)
t.texture('p', color = (0,0,1), opacity = 1, specular = .2)
t.sphere((-1,-.57735,-0.7071),1,'s')
t.sphere((1,-.57735,-0.7071),1,'s')
t.sphere((0,1.15465,-0.7071),1,'s')
t.sphere((0,0,0.9259),1,'s')
t.plane((0,0,-1.9259),(0,0,1),'p')
t.show()

Cobwebs

This is mostly a test of how code looks in a post. The sage functions below are some initial attempts to write basic dynamical systems plotting functions for educational purposes.

def cobweb(a_function, start, mask = 0, iterations = 20, xmin = 0, xmax = 1):
'''
Returns a graphics object of a plot of the function and a cobweb trajectory starting from the value start.

INPUT:
a_function: a function of one variable
start: the starting value of the iteration
mask: (optional) the number of initial iterates to ignore
iterations: (optional) the number of iterations to draw, following the masked iterations
xmin: (optional) the lower end of the plotted interval
xmax: (optional) the upper end of the plotted interval

EXAMPLES:
sage: f = lambda x: 3.9*x*(1-x)
sage: show(cobweb(f,.01,iterations=200), xmin = 0, xmax = 1, ymin=0)

Note: This is very slow with symbolic functions.
'''
basic_plot = plot(a_function, xmin = xmin, xmax = xmax)
id_plot = plot(lambda x: x, xmin = xmin, xmax = xmax)
iter_list = []
current = start
for i in range(mask):
current = a_function(current)
for i in range(iterations):
iter_list.append([current,a_function(current)])
current = a_function(current)
iter_list.append([current,current])
cobweb = line(iter_list)
return basic_plot + id_plot + cobweb

def orbit_diagram(a_function,parameter_interval, domain=[0,1], mask = 50, iterations = 200, param_num = 500.0):
'''
Returns a plot of the iterations of a function as a function of a parameter value.

INPUT:
a_function: a function of one variable
parameter_interval: a two-element list of the lowest and highest parameters to plot.
domain: (optional) a two-element list of the lowest and highest input values to iterate
mask: (optional) the number of initial iterates to ignore
iterations: (optional) the number of iterations to draw, following the masked iterations

EXAMPLES:
sage: f = lambda x,m: m*x*(1-x)
sage: show(orbit_diagram(f,[3.4,4], mask = 100, iterations = 500), xmin=3.4, ymin=0)

NOTES:
This is pretty crude so far.
'''
point_list = []
plen = RDF(parameter_interval[1] - parameter_interval[0])
seed = random()*(domain[1]-domain[0])+domain[0]
for i in srange(parameter_interval[0],parameter_interval[1],plen/param_num):
for x in range(mask):
seed = a_function(seed,i)
for x in srange(iterations):
seed = a_function(seed,i)
point_list.append((i,seed))
return point(point_list,pointsize=1,rgbcolor=(0,0,0))

A blog by any other name...

I've always hated the word 'blog', which I think is the only reason I haven't started one until now. But now with the sage projects new blog, I guess I'll give it a try.