Numerical derivatives in parallel

My last post was about a fairly esoteric new language feature in python: coroutines. I wasn't reading up on them because they're whiz-bang cool (okay, not only because they're whiz-bang cool) but because I have an actual problem that they might help with. For the triple system, I find I need to do very slow numerical optimization problems. That is, the objective function takes seconds to evaluate. It can't be parallelized internally, and numerical minimization is largely a sequential operation: evaluate the function, figure out which way is downhill, go there, repeat. So this code currently uses just one core, and is almost as slow as the highly parallel Markov-Chain Monte Carlo code that does hundreds of times as much work. So obviously it would be nice to find some parallelism. I have a way to do that: parallelize gradient evaluations.

The other reason for this post is to make a point about code robustness. Getting numerical derivatives right is a notoriously difficult process, and getting them both efficient and right is even harder. The package I have been using, MINUIT, has been used by the particle physicists for decades, and they have hammered out all sorts of corner cases that a naive implementation might miss. So I'm not going to use a naive implementation, not least because my problem is numerically difficult. This means the code is not simple, and I don't want to try to restructure it manually to be more concurrent; I want the computer to handle the concurrency for me.

Numerical coroutines

Python 3.5 has new syntax for coroutines. Like many of python3's new features (unicode strings? who uses non-English languages in astronomy?), it is not immediately obvious that coroutines are of any use to people like me who mostly write relatively simple imperative scripts to wrangle large amounts of data or long computations. But cores aren't getting any faster, so we're going to have to start writing parallel code if we're going to take advantage of new machines. Now, if you want to take a gigantic Fourier transform, or do algebra on enormous matrices, you are going to need a witch design you a parallel algorithm. But our code is often full of easy opportunities for parallelism that we don't take advantage of because it's awkward. So I keep my eye out for ways to take advantage of this easy parallelism. I have found one, in the combination of coroutines and futures, and I'd like to describe it with a worked example. Maybe it'll sell you on some of these cool new language features!

Parallel ipython notebooks and the H test

The CPU on the nice shiny new server I log in to is really not much faster than that of the ratty old laptop I have in front of me. The server has more memory, and more disk space, but CPU-wise it just has more not faster. For that matter even my laptop has two cores. So if I have some heavy-duty computing task, I'd better find a way to make it use multiple cores in parallel. Some tasks are just plain hard to parallelize (solving ordinary differential equations, for example), but it turns up fairly often that I'm doing something embarrassingly parallel: there's a for loop somewhere that just does the same thing to lots of different pieces of input. If only it were easy to hand each piece to a different core! Well, there are various tools for doing this sort of thing, but most of them apply to scripts or programs that run non-interactively. It turns out that ipython offers tools for interactive parallel computing. I'm going to explain how I use them, by working through a test problem, checking some statistics on a periodicity test (the H test).

Remote ipython notebooks

I use ipython notebooks a lot. For interactive tinkering, for calculations with units, for interactive control of long-running processes, for easy parallelization (details to come), and to generate plots for all my papers. When I can, I run the notebooks on a server rather than my rather wimpy laptop; apart from the obvious advantages of more cores, more disk space, and more RAM, this also means I can leave things running for ages. It's all very convenient, but it requires some setup. As with the ssh configuration, I've found myself explaining the setup to several different people recently, so I thought I'd put it here. (Still working on that cute-kittens post.)

SSH incantations

I use ssh all the time, for shell access on remote machines, to transfer files, to carry VNC sessions, and to use different institutional affiliations to get at paywalled things. Over the years I have developed a particular way of setting it up that minimizes the pain involved in doing all of the above things. Several times in the past months I have tried to explain the whole setup to someone, so I think that's a sign it's time to post it here. I'll try to post something about, I don't know, fluffy kittens, boating on the Thames, seaborgium hexacarbonyl, or something soon for those of you who for some reason expected an interesting blog.

The Fortress of the North

I live in Groningen. It's a city of about 200,000, about 50,000 of whom are students at the universities. So I sort of thought it might resemble Kitchener-Waterloo, the college town where I did my undergraduate degree, but in fact I've never been anywhere quite like this before.

The problem with "beginner's mind"

I should say up front that I know nothing about authentic Zen Buddhism. But one of the concepts that has made it into popular culture that I rather like is "beginner's mind": an expert knows what she thinks about things, while a beginner sees everything anew and judges it on what it is. It sounds kind of charming, and like a recipe for creativity. But applying it too literally can get you in trouble.

Concretely, ASTRON (where I work) arranged for a company to come in and teach programming courses to anyone who wanted them. They taught three: Introductory Python, Numeric Python, and Advanced Python. Now, I think that's a brilliant idea: we astronomers all spend most of our time writing code at one level or another, and nobody seems to have bothered to teach us how to do it well. So that ASTRON (actually NWO I think) cares is a really good sign. So of course I wanted to participate. I have lots to learn about writing good programs. I don't think I'm an expert programmer, or even a Python expert, but that latter is partly because I try not to think of myself as an expert in anything. I knew the Introductory Python course would not be productive — I write python code every day. And the Numerical Python, again, was a good idea, but having written, for example, the reshape function in numpy, and the spatial module in scipy, I figured that was probably not going to be too productive either. But the Advanced Python course sounded promising. And I wanted to show my support for the whole idea. So I signed up.

Well, the course was the last three days, and it was a pretty good course, but not at all what I needed. Which should not have been much of a surprise: when I stopped to do the math, I realize that I wrote my first python program eighteen years ago. (Good God.) Still, it was interesting to see how they ran the course, and to think about how I would run one (because if I end up somewhere like McGill that has basically nothing for physics grad students, I will run one, official or not).

Visualizing the new pet

I recently wrote about a new object I am studying: a millisecond pulsar with two white dwarf companions. There is lots more I want to say about it, but I think it would be nice to show what it looks like, or at least, to show a video I made trying to make visible what's going on:

Edited to note that Blogger's YouTube embedding is distinctly flaky; video is here. 

This video shows the orbital motions in the triple system. The orbits are drawn to scale, showing the actual motions of the two stars (red and yellow) and the pulsar (white). The first ten seconds are played relatively slowly, showing the motion around the inner orbit, then we speed up to see the motion around the outer orbit. For a sense of the time scale, an "MJD" is a modified Julian day, so a single day long. The larger left panel shows all three bodies, with trails marking the motion of the outer companion and the center of mass of the inner system. The inset in the top right zooms in on the inner system, showing the pulsar and the companion, with trails marking their orbits. The dots that appear on the orbits mark moments when we have observations of the system, color-coded by telescope; it should be clear that we have quite thorough coverage of both orbits. Each measurement tells us the distance to the pulsar to within a kilometer, so that we can measure the tiny deviations of these orbits from perfect Keplerian ellipses, allowing us to reconstruct the orbit.

There's a little more to it than that.

New pet: PSR J0337+17

I did my PhD thesis on PSR J1023+0038, a millisecond pulsar that is at a fascinating point in its evolution. (In fact there have been developments since the thesis was submitted; more about that later.) But during a moment of procrastination, I got involved with a new and fascinating system. The name, unedifying as usual, is PSR J0337+17, and it is unique in that the pulsar has not just one white dwarf companion but two.

Flywheel energy storage

In the quest for something better to run our cars on than gasoline, one of the proposals is flywheels. In fact, for a while there were flywheel-powered buses running in Switzerland and Belgium. On one level, it makes a lot of sense: you're storing energy as mechanical motion, and we're pretty good at transmitting mechanical motion from place to place. On another level it scares the living daylights out of me: a car in motion uses tens of kilowatts, so the car must be able to store hundreds of kilowatt-hours. If you let all those loose at once bad things will happen: 100 g of TNT going off releases about a hundred kilowatt-hours. Fortunately it's hard to get gasoline to do this, but a flywheel is just itching to dump all its energy. Batteries are a little scary too, to be honest. But anyway, that's all a digression: I want to talk about some really staggering examples of flywheel energy storage: pulsars and black holes.

New job, new country

I just started a new job, and in the way of academic life, I had to move. I put some things in storage, hopped on a plane with two suitcases, and hey presto, I live in the Netherlands now. The new job is with ASTRON, the Dutch institute for radio astronomy, and it's going really well. Living in a new country actually takes more getting used to.

For one thing, in spite of the fact that I live in arguably the most pedestrianized city in Europe, and in a densely-populated country with good transit systems and cities built when rapid transit meant corn-fed horses, I find myself tempted to carpool to work. See, the thing is, ASTRON was built in the middle of a national park in order to minimize interference with the radio telescope that was also being built. Since the fifties, though, radio astronomy has made a few strides, and the telescope is now used chiefly in education. The telescopes we do use are all off-site, so we're in the middle of the park for historical reasons. It means we're twelve kilometers from the nearest train station, and even the bus drops you off in the village of Lhee, a kilometer or so away from work. So a ride to work is tempting. Fortunately, this is the Netherlands, so I've been biking.

Differentiating the solutions to differential equations

A kind of problem that turns up quite often in physics and elsewhere is to find the solution of an ordinary differential equation, given some initial conditions. That is, you have a system with some state $x$ represented as a vector of real numbers, you have an initial state $x_0$, and you have a rule describing the evolution of the state: $$ \frac{dx}{dt} = F(x,t) $$ And your goal is to find $x(t)$. This standard problem has some standard solution techniques, some quite advanced - Runge-Kutta methods, symplectic methods, Hermite integrators. A few are implemented in scipy. But it sometimes happens that solving this problem is only part of a more complicated process, say of fitting, where it would be nice to have the derivatives of the solution with respect to the various initial conditions. It turns out this isn't too hard to work out, usually.

Weighted Poisson Uncertainties

I recently ran across a rather awkward mathematical problem. I'm trying to make a histogram of photon arrival phases, complete with an uncertainty on the number of photons in each bin. Normally this is done by just taking the square root of the number of photons, which is at least approximately right based on Poisson statistics. But in my problem — data from the Fermi space telescope, which is not very good at localizing low-energy gamma rays — the photons are weighted: for each photon I have a probability that it really came from the source. So the values in the histogram should be the total probability. But what should the uncertainty be? The short version is: the square root of the sum of the squares of the weights.

ETA: This is in the literature, without justification as far as I can tell. See below.

git tear-hair

I use git. It's fast, it's convenient, it gives you github, and mostly it works. But several times now I've managed to bollix a git repository, been faced with impenetrable messages, and been unable to continue using git unless I sort them out. This latest time I had enough stubbornness to figure out how to un-bollix my repository, and I wanted to record it for posterity (specifically, for people whose problem-solving technique involves Googling for the specific error message). The error I received was:

error: object file .git/objects/86/5d2dffe9a3d72917934ed9693c7167efb6d8d5 is empty
fatal: loose object 865d2dffe9a3d72917934ed9693c7167efb6d8d5 (stored in .git/objects/86/5d2dffe9a3d72917934ed9693c7167efb6d8d5) is corrupt

Read on for how it happened, my understanding of what it means, and how I fixed it. Or skip this even-more-technical-than-usual post; I'll try to post something with corals or kittens soon.

Justice

Available from Amazon

I recently read (listened to? What's the right way to say that for an audio book?) "Justice: What's the right thing to do?" by Michael J. Sandel. It's an interesting book, though rather heavier going than Nate Silver's breezy book on prediction. The basic question it addresses is, how do you decide what's right? Not just what's legal — you need some basis for deciding which laws are just. He talks through several approaches, including some ideas of his own. While he talks about several currently-debated topics, he refrains from stating his position on any of them, instead pointing out what the conflicting ideas of justice are that motivate the two sides.

Talking about radio astronomy

McGill has a modest observatory on the roof, with a 14-inch (optical) telescope. It was installed quite some time ago, then left to molder for a few years. This annoyed me, just on principle, so I talked my way into rehabilitating it as best I could, with my radio astronomer's skills. Even with it working again, it was just used by the occasional grad student who wanted to take her friends up to see Jupiter, or the surface of the moon, or whatever. Fortunately, we recently acquired an enthusiastic and organized post-doc who pushed hard and set up a public outreach program. Typical nights start with a public lecture, then we take folks outside to two portable telescopes and the bigger one on the roof to look at the stars. Of course, we can't predict the weather, so we have a few demos we can do inside - comet-making, liquid nitrogen, a muon detector, but really the appeal is going out and looking up at the sky. But especially in the summer, we need to entertain the public until it gets dark. Since it's once a month, we have a perennial need for speakers. So I volunteered, to give a one-hour talk on radio astronomy for the general public. I had never given a one-hour talk before, and radio astronomy doesn't produce as many pretty pictures as I would like, but I think it came out pretty well. And we had a videographer who recorded the whole thing, so if you're curious, here it is.

Prediction

The book can be bought from Amazon

I recently finished reading Nate Silver's book "The Signal and The Noise: Why Most Predictions Fail But Some Don't". Nate Silver is currently somewhat famous because his election predictions - he runs the blog fivethirtyeight for the New York Times - his election predictions were basically spot on. The book, timed to come out shortly before the election (in what I think was a sort of gamble to capitalize on a successful prediction) is about prediction, generally defined. Structured as a sort of collection of case studies, the book talks about all the ways prediction can go wrong, and about what to do to try to get it right. There's no nice one-line answer, but I think he gives some pretty good advice.

The Maya-Muon experiment

This week's Friday colloquium talk really had me feeling like I was in a science fiction show. The speaker is building machines to use cosmic-ray muons to see the interiors of still-buried Mayan ruins. Does this not sound like a line of technobabble from Stargate SG-1? But it's nearly feasible.

Higgs liveblog

A candidate Higgs event; see below.

The Higgs detection (?) announcement is being broadcast live. So I thought I'd blog my comments as it unfolds.

Fireworks

In honor of Canada Day (and for my neighbours to the south, the fourth of July) here's a video:

While this might look like a meteor, and in fact it is asteroidal material falling to Earth, it's actually some extremely expensive fireworks. It's the Hayabusa spacecraft returning to Earth after visiting asteroid 25143 Itokawa, and the material it brought back was safely encapsulated.