The Singular Value Decomposition and Congressional Voting

I am teaching a class about the SVD (Singular Value Decomposition) of a matrix this week.  I was inspired by a nice article of Carla Martin and Mason Porter, "The extraordinary SVD", to compute the SVD of the voting record of the 112th Congress (House of Representatives) to show to my class.

If you're interested in how this is done, here is the Sage code I used as a Sage worksheet. 

Here is the projection onto the first two singular vectors:

Review of "Sage: Beginner's Guide" by Craig Finch

I was asked by the publisher to review "Sage: Beginner's Guide" by Craig Finch. They sent me a free copy; I have no other conflicts of interest. I am generally biased towards Sage itself, as an avid user and minor developer.

On Amazon you can browse the table of contents, which gives a pretty good idea of the strengths of the book, namely basic computation and plotting, numerical calculations, and data analysis. The focus was an excellent choice considering what is already available. The current free Sage Tutorial is oriented much more towards pure mathematicians. There is a Numerical Computing With Sage as part of the standard documentation, but at the moment its quite short and nowhere near as helpful as Finch's book.

I liked the style of the book a lot. There are many code examples that illustrate how to accomplish concrete tasks, along with good explanations of what they are doing. Many of these are things that are unfortunately far from obvious to a beginner (or even intermediate) Sage user. Despite using Sage heavily for the last five years, I learned some new things. The book is particularly strong in showing how to use Numpy, Scipy, and Matplotlib. Sage wraps a lot of the functionality of these projects, but if you want to do something that isn't included in the standard interfaces it can be quite mystifying.

Chapter 9, "Learning Advanced Python Programming", might have been a little ambitious. There's nothing wrong with it, but its too short to provide enough. Fortunately there are a lot of good books, some of them free, that cover Python programming in much more depth. I would have preferred some of this space and effort to be devoted to using Cython and the @interact command, which are covered very briefly in Chapter 10.

I teach several classes using Sage and I will definitely advertise this text as a useful optional supplement (I consider it a little too expensive to add on as a mandatory second text). It would be nice if some institutions considered using Sage instead of its commercial competitors such as Maple, Matlab, and Mathematica - you could probably give every student a copy of this book for the money saved from license fees!

The only thing I disliked about the book was the quality of the illustrations. Sage output that was in LaTeX was not typeset, but instead looks as if a PNG was copied from a screenshot. Some of the examples would have benefited from being in color. The quality of the plots is also somewhat poor. This is not too big a deal if one is following along with Sage, since you can reproduce the figures. None of them are bad enough to obscure the content.

Overall this is a very impressive and useful introduction to Sage that should help any beginning user a great deal.

A family of Mobius transformations

I was trying to sketch the behavior of a Mobius transformation in my complex analysis course today. Its hard to convey on the blackboard, so I tried making a video which shows a homotopy of the image of the unit disk, from the identity to (1+z)/(1-z) and back again.

Plotting the zeta function

I lectured a tiny bit on the Riemann zeta function for the first time in my complex analysis course, which inspired me to make the following plot. Colors are the argument of the Riemann zeta function, brightness is proportional to magnitude. The default brightness map gives very low contrast, so I modified the magnitudes. To help with seeing the magnitudes a contour map with exponentially spaced contours is overlaid.

def xy_to_zeta_size(x,y):
    return abs(zeta(N(x+I*y)))
cvals = [e^i for i in srange(-7,1,.25)]
cp = contour_plot(xy_to_zeta_size,(-6,3),(-3,3),contours = cvals, fill=False, plot_points = 201)
rzeta(z) = zeta(z)/norm(zeta(z))^(.25)
rzf = fast_callable(rzeta,domain=CDF)
cparg = complex_plot(rzf,(-6,3),(-3,3))
show(cparg+cp,figsize=[18,12])

Youtube videos of n-body choreographies

For a presentation I'm doing tomorrow, I made some animations of n-body choreographies (n equal masses which share a common trajectory). I put these up on youtube since that seems to make more sense than having them mangled by Blogger or hosting them myself:

Twenty-one bodies
Eight bodies
Three bodies on a figure eight

The videos are done using Sage, the Tachyon raytracer, and ffmpeg. The orbit data is from Carlos Simo.

Here's a 2D picture of the twenty-body choreography.

Torus knots in Sage with Tachyon and @parallel

After working a bit more on a patch improving the Tachyon raytracer in Sage, I made a movie of a torus knot. The code below was parallelized to use 20 cores at a time on one of the Sage Foundation's lovely Sun Fire X4450s. Frames were animated with ffmpeg. Blogger might degrade the video, so I recommend downloading the original here.

pk = 3
qk = 23
res = 180
trad = .5
lab = 'tk_' + str(pk) + '_' + str(qk)
def tknot(theta):
    return [(2 + trad*cos(qk*theta/pk))*cos(theta), \\
(2 + trad*cos(qk*theta/pk))*sin(theta), trad*sin(qk*theta/pk)]
@parallel(ncpus = 20)
def tknotter(anum):
    th = anum*2.0*pi/res
    T = Tachyon(xres = 800, yres = 600, camera_center = (4,4,2.5), \\
raydepth=12, antialiasing=2)
    loc = tknot(th)
    nloc = [abs(x)/norm(vector(loc)) for x in loc]
    T.texture('s1',color=nloc,specular=1)
    T.texture('p1',color=(.5,.5,.5))
    T.sphere(loc,.05,'s1')
    ci = 0
    for sth in srange(0,th,2*pi/res):
        e1 = tknot(sth)
        e2 = tknot(sth+2*pi/res)
        ts = [abs(x)/norm(vector(e1)) for x in e1]
        T.texture('t'+str(ci), color=ts)
        T.fcylinder(e1, e2, .025,'t'+str(ci))
        ci += 1
    T.plane((0,0,-4),(0,0,1),'p1')
    T.light((6,6,11),.1,(1,1,1))
    T.light((5.9,5.9,11),.1,(1,1,1))
    T.save(DATA+lab+'%03d.png'%anum)
    return 1

qin = range(pk*res)
s = list(tknotter(qin))

for i in range(40):
    r2 = os.system('cp '+DATA+lab+'%03d.png '%(pk*res-1) + DATA+lab+'%03d.png'%(pk*res+i))

anum = qin[-1] + 41
th = anum*2*pi/res
T = Tachyon(xres = 800, yres = 600, camera_center = (5,5,3))
T.save(DATA+lab+'%03d.png'%anum)

The animation was then done with the following code; if you want to do something similar you would have to use the appropriate path to your ffmpeg executable.

import subprocess
subprocess.call('/home/mhampton/ffmpeg/ffmpeg -qmax 2 -i ' + \\
DATA + 'v%3d.png ./t2vest.mp4', shell=True, stdout=file('/dev/null','w'),stderr=file('/dev/null','w'))

Tachyon in Sage: depth of field perspective

Today I found a bug in the options to the tachyon raytracer in Sage. It should be possible to use alternate projections, including fisheye lens and perspective depth-of-field projections. I think I've fixed that, patch should get into sage-4.6 or so.

Here's a quick example:

T = Tachyon(xres=800,antialiasing=24, raydepth=12, projection = 'perspective_dof', focallength = '1.0', aperture = '.005')
T.light((0,5,7),1.0,(1,1,1))
T.texture('t1', opacity=1, specular = .3)
T.texture('t2', opacity=1, specular = .3, color = (0,0,1))
T.texture('t3', opacity = 1, specular = 1, color = (1,.8,1), diffuse=0.2)
T.plane((0,0,-1),(0,0,1),'t3')
ttlist = ['t1','t2']
tt = 't1'
T.cylinder((0,0,.1),(1,1/3,0),.05,'t3')
for q in srange(-3,100,.15):
    if tt == 't1':
        tt = 't2'
    else:
        tt = 't1'
    T.sphere((q,q/3+.3*sin(3*q),.1+.3*cos(3*q)), .1, tt)
T.show()

complex plot in sage

I've been amusing myself all day with doing some complex plotting in Sage. The default color scheme has 0=black, infinity=white, and red=real. A lot of the time, this doesn't work too well, but if you plot f/abs(f) the color will just depend on the complex argument. For example, if f = (z^5+1)/(z^5-1) we get:



Or a more baroque example, f = z*log(z)*exp((z^5+1)/(z^5-1)^(1/2)):




Sometimes the default coloring is helpful, for example when showing the essential singularity at 0 of exp(1/z):

Plucker's Quartic Curve

A quartic curve in the plane has a maximum of 28 real bitangents. Plucker came up with the following sort of example that attains that bound, which I illustrated in a particular case with phcpack and Sage:

More polytopal animations

I was inspired by a still image from the Wikipedia entry on the 24-cell to make an animation of rotating and then stereographically projecting the 24-cell, showing the images of the original polytope edges. Math on Wikipedia is pretty spotty, but the polytope entries have really become top-notch in the last few years.
Here's a still from the animation:



and a medium and small resolution version of the animation:

small (300x300, 6 MB)
medium (512x512, 17 MB).

600-cell animations

Over the last year or so, when I have time, I've been trying to fine tune ray-traced animations of polytopes (using Sage, Tachyon, ffmpeg, and cddlib). My latest effort was to do a loop of Schlegel projections of the beautiful 600-cell. The 600-cell is one of the six regular polytopes in 4 dimensions; its usually considered the 4D analog of the icosahedron. Here's a still from my animation:



Since Blogger's re-encoding butchers my movies, I will just post links. Here are four sizes:

small (320x200, 3.5 MB)
medium (640x400, 7 MB)
large (1280x800, 13 MB)
large, high-res (1280x800, 33 MB).

24 cell movie

This is basically a test of my workflow for a more ambitious project, but its somewhat amusing in its own right. The original 48 MB version is here (much better quality).

Hearing a trigonometric identity

If you force a harmonic oscillator with natural frequency w_0 with a sinusoidal force of frequency w, the resulting steady-state is a linear combination of those two frequencies. If its something in the audible range, you hear either the two seperate frequencies or one frequency with a beat if they are close together. In other words, you hear either the right or left hand side of the identity shown below.



Sage source code available here. Ironically, because of my use of html I can't seem to show the source directly here.

Small frequency difference.

Smallish frequency difference (close to the Ma Bell tone).

Intermediate frequency difference.

Large frequency difference - you definitely hear it as a 2-note chord.

symmetry in chaos

Recently SIAM had a deal to get some books from their catalog cheap if you bought the new edition of "Symmetry in Chaos" by Field and Golubitsky. So I did, and couldn't help but try to reproduce some of their figures. I think they might have a typo in their parameters for their figure 2.3 (or I am making some mistake, quite likely), but after exploring a bit with other parameters I got the following, with iterates of

Polytopes in Sage, take 1

Here's a first attempt at a short video intro to polytopes in Sage. Sage 4.0, which should be released soon, has some new functionality that was added during Sage Days 15 (thanks to David Perkinson for reviewing some of that, and Rob Beezer for helping David learn the review ropes).

The video is about 4 MB and can be found here. I am not sure how the format will work on other computers so if anyone takes a look and it doesn't work let me know. Also, I am interested in getting constructive criticism on the video.

Slouching towards Adelaide

At some point when I was in grad school I became aware of some work on symmetric Venn diagrams. If you google this, you will find this link, which has been maintained but not changed too much since 1997. Other than that, there isn't a lot on the web and there is a particular lack of quantitative direction on how to construct the beautiful rotationally symmetric Venn diagrams such as Adelaide. This was named by Anthony Edwards after the city in which he discovered it. I have always wanted to go to Adelaide, it holds a strange attraction for me, so perhaps that is why that particular 7-set Venn diagram stuck in my head.

When I started working on my coloring book, I immediately thought of the Adelaide diagram but I didn't know how to construct it. After some mistakes today, I think I finally have it down. Here is a colored version:



Code (in Sage) for some version of this will be in the final coloring book.

Cayley Cubic

To further my polynomial education I've been thinking some lately about the Cayley cubic, which in the affine form I was using is given by x^2 + y^2 + z^2 + z(y^2-x^2) -1 = 0.

I made a little movie of some of the real solution set of the Cayley cubic. Here's one of the stills:

I think the quicktime version works better in browsers, or you can download a slightly higher-quality mp4.

Curves of Pursuit

When I visited Stellenbosch University in February, I noticed the math department had some nice posters on various mathematical topics that I thought could be implemented in Sage. The only one I did was "Curves of Pursuit", in which one has n points in a regular n-gon, and each point chases the next. The solutions are prettier when discrete time steps are taken. You can see the code at

http://wiki.sagemath.org/interact/graphics#CurvesofPursuit

Carl Witty does Sage on the G1

"the entire compile took 20943 minutes of real time (about
14.5 days)"

This is just heroic, borderline legendary: compiling Sage on the Android G1. There is nothing to be done but read the post from sage-devel and see the picture.

I have compiled and run Sage 3.2.3 on my T-Mobile G1 cell phone, and
large portions of it actually work.

300 files had failing doctests; this means that all doctests passed in
864 files. A lot of the failing doctests are with pexpect (maxima,
gap, etc.); I don't know why these fail. When I try simple things
with gap and maxima, they do work. Many more tests fail due to
timeouts; the 300-second timeout is far too short for the G1.
(Doctesting any file with doctests takes >50 seconds; doctesting a
file with no doctests takes >4 seconds. The entire doctest run took
270724 seconds, or a little over 3 days.)

Unfortunately, it's far too slow to be useful for anything (just
starting Sage takes about a minute, if there's enough free memory; but
since Android likes to keep the memory full all the time, there never
is enough free memory, and it takes much longer to start). So I don't
plan to do anything further with the project.

Of course, the sensible way to do this is to find some nice fast
computer and set up a cross-compiler. I didn't do this the sensible
way; I actually did the build on the cell phone. And of course, the
build would stop every once in a while due to a bug, which I would
then patch around. So, adding up the 14 chunks of compilation
involved, the entire compile took 20943 minutes of real time (about
14.5 days) and 9438 minutes (about 6.5 days) of CPU time. You'll note
that the CPU time adds up to less than half of the real time; most of
the rest of the time was spent swapping. The G1 has about 100MB of
RAM; I set up a swap file on my micro-SD card, to allow the build to
proceed at all. On several files, the compiler uses more than 300MB.
While compiling those files, the CPU is typically less than 1% active;
I'm pretty sure that there were files that took more than a day to
compile. The phone was essentially unusable for smartphone activities
(web browsing, etc.) while the build was running. In fact, the whole
smartphone UI crashed and restarted on a fairly regular basis during
the build (I'm guessing it has some sort of watchdog timer and reboots
itself if it detects that some operation is taking an unreasonably
long time), but the build just continued anyway (running inside a
screen session).

The ATLAS tuning process took about 3.2 days (with almost as much CPU
time as real time; that didn't swap much). Skimming through the logs,
I see reports of numbers like 3 to 6 MFLOPs.

The build was performed inside a Debian armel testing chroot; it's
pretty nice being able to run real Debian on my cell phone.
(Everything works; I can ssh in and out, etc., although when it's on
the cell phone network, inbound connections are problematic because
it's behind a NAT.)

To make the whole project even more pointless, since the phone is
running real Debian (inside the chroot), I can probably just wait a
few weeks or months for Tim Abbott to get all the portability issues
fixed in the sagemath Debian package, and then install Sage on the
phone with apt-get.

Here's the processor:

cwitty@localhost:/tmp$ cat /proc/cpuinfo
Processor : ARMv6-compatible processor rev 2 (v6l)
BogoMIPS : 383.38
Features : swp half thumb fastmult edsp java
CPU implementer : 0x41
CPU architecture: 6TEJ
CPU variant : 0x1
CPU part : 0xb36
CPU revision : 2
Cache type : write-back
Cache clean : cp15 c7 ops
Cache lockdown : format C
Cache format : Harvard
I size : 32768
I assoc : 4
I line length : 32
I sets : 256
D size : 32768
D assoc : 4
D line length : 32
D sets : 256

Hardware : trout
Revision : 0080
Serial : 0000000000000000

Carl

Live from Cape Town

I arrived at AIMS (the African Institute for Mathematical Sciences) yesterday. Its a great place, right on the Muizenberg beach near Cape Town - quite a thrill after being in Duluth this time of year.

I have volunteered to teach a short course on bioinformatics algorithms here, using biopython and Sage as the computational platform. It will be very interesting to see how it goes over. I am using the Plasmodium species which cause malaria as a source of examples and projects.

Below is a satellite picture of the cape peninsula and a view of Muizenberg that includes the main AIMS building.