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.

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.

NuSTAR about to launch

NASA is getting ready to launch a new X-ray satellite, NuSTAR. This satellite will observe harder X-rays than we are currently able to. More specifically, this satellite is able to produce high-resolution images of the hard X-rays that come from the decay of isotopes produced in supernovae; the ability to observe hard X-rays with high spectral resolution should also let us probe neutron star formation. Plus it may help us pin down that hard X-ray emission from anomalous X-ray pulsars that nobody understands.

If all this sounds hopelessly technical, sorry. It's still impressive to me that they're going to take this hundred million dollar spaceship, fill it with what amounts to high explosives, strap it underneath a discontinued Lockheed airplane, fly it up into the air near Kwajalein, then drop it and light the fuse. The rocket will then burn for a few hundred seconds, dropping several stages and winding up in low Earth orbit. Once there, the mast holding the mirrors will extend from its one-meter storage canister to its full ten-meter length, the satellite will extend its solar arrays, and scientists will begin debugging and calibrating the instruments. That's if everything goes well; no Pegasus launcher (which this is) has ever actually exploded, though one pair of satellites failed to exit its fairing. So there are some people chewing fingernails, and there's some real drama here. You can watch it online.

Sorry for the hurried nature of the post; I'm planning to watch it live and report on the results as soon as we know.

Edited to add: Success! At least so far. The satellite is in orbit and the solar panels work. The last big worry-point is the extension of the boom; the ten-meter beam that holds the mirrors in position relative to the cameras is currently folded into a one-meter canister. Supposedly, the scientists were reassured that "we can't give you any details, but we know how" to build such a boom, but it won't extend for another week. That'll be followed by a few weeks worth of testing and calibration. But so far so good.

Local Maxima and social annealing

Dan's Data has a new post/magazine article up, about local maxima. Frustratingly, there's no comment section, but I've been reading some interesting things that seem relevant, so I'll post them here. (TL;DR of Dan's post: socially we seem to get stuck in local maxima instead of finding the real best way to do things.)

Testing whether a signal is broad-band

Radio pulsars are generally broad-band sources — you can hear them over a very wide range of radio frequencies, from about a hundred megahertz up to, in a few cases, tens of gigahertz. Their emission does change with frequency, decreasing as a power law, but over reasonable bandwidths we expect to see their signal in all frequency channels.


Radio-frequency interference ("RFI"), on the other hand, is very frequently narrow-band, appearing in just one, or a few frequency channels. In fact, one way we try to manage RFI is by using really wide bandwidths, so that there's so much power from the pulsar signal that narrow-band RFI is drowned out. Unfortunately, all too often the RFI is so strong that even a narrow-band signal can dominate the total power. And since it's narrow-band, dedispersion doesn't smear it out like it would a broad-band interference spike. So narrow-band RFI is one of the kinds of interference that is particularly hard to sift out from pulsar candidates.

Just recently, on the arxiv, a paper came out that attempts to address the problem of testing whether a candidate is broad-band:

Multimoment Radio Transient Detection
Laura Spitler, Jim Cordes, Shami Chatterjee, Julia Stone
We present a multimoment technique for signal classification and apply it to the detection of fast radio transients in incoherently dedispersed data. Specifically, we define a spectral modulation index in terms of the fractional variation in intensity across a spectrum. A signal whose intensity is distributed evenly across the entire band has a much lower modulation index than a spectrum with the same intensity localized in a single channel. We are interested in broadband pulses and use the modulation index to excise narrowband radio frequency interference (RFI) by applying a modulation index threshold above which candidate events are removed. The technique is tested both with simulations and using data from sources of known radio pulses (RRAT J1928+15 and giant pulses from the Crab pulsar). We find that our technique is effective at eliminating not only narrowband RFI but also spurious signals from bright, real pulses that are dedispersed at incorrect dispersion measures. The method is generalized to coherent dedispersion, image cubes, and astrophysical narrowband signals that are steady in time. We suggest that the modulation index, along with other statistics using higher-order moments, should be incorporated into signal detection pipelines to characterize and classify signals.

This looks very promising, but there's some testing I wish they'd done. More on the subject below the jump.

Faster-than-light neutrinos: keeping time

There's been a recent announcement of evidence that muon neutrinos may travel faster than light. That would be really weird - in particular it would pose serious problems for Einstein's theories of relativity. But even the discoverers aren't ready to claim that; they just describe their results and say they're puzzling. Personally I think it's unlikely they're right, but figuring out why not may be very interesting.


Faster than Light? par CNRS

Fortunately for us, they have posted a preprint of their paper on arxiv.org. All images below are from that paper.

I can't say much about the particle physics, or the details of instrument calibration, but I can address one possible way people have suggested the result may be wrong: inaccuracies in the time standards at the two endpoints.

Least squares and timing noise, part 2

Simulated time series

In my previous post I described a new paper about fitting pulsar parameters in the presence of timing noise using a general least-squares method. It seems like a good approach, but I'd like to look at it in more detail. So: python to the rescue!

Least-squares and timing noise

Figure 4 from the paper: residuals and spectrum

Working with long-term pulsar timing data sets is a nuisance because of so-called "timing noise". Not only is this noise above and beyond the usual observational uncertainties, perhaps because it is torque noise, it tends to be strongest at low frequencies (it is very "red"). Often so much so that leakage from the very lowest frequencies dominates at all analysis frequencies. Various people, myself included, have tried various approaches for dealing with this noise, but a recent arxiv paper shows real promise:

Pulsar timing analysis in the presence of correlated noise

W. Coles, G. Hobbs, D. J. Champion, R. N. Manchester, J. P. W. Verbiest

Pulsar timing observations are usually analysed with least-square-fitting procedures under the assumption that the timing residuals are uncorrelated (statistically "white"). Pulsar observers are well aware that this assumption often breaks down and causes severe errors in estimating the parameters of the timing model and their uncertainties. Ad hoc methods for minimizing these errors have been developed, but we show that they are far from optimal. Compensation for temporal correlation can be done optimally if the covariance matrix of the residuals is known using a linear transformation that whitens both the residuals and the timing model. We adopt a transformation based on the Cholesky decomposition of the covariance matrix, but the transformation is not unique. We show how to estimate the covariance matrix with sufficient accuracy to optimize the pulsar timing analysis. We also show how to apply this procedure to estimate the spectrum of any time series with a steep red power-law spectrum, including those with irregular sampling and variable error bars, which are otherwise very difficult to analyse.

I'd like to look at it in more detail and try some of the techniques on test data.

How tempo2 does its fitting

Pulsars can be difficult objects to study: for example, their radio pulses can randomly change in brightness, turn off, turn on, change shape, and we really don't know why. Nevertheless there has been some excellent science done by studying those very radio pulses. The trick has mostly been to simply not care how bright they are or what shape they have and focus on when they arrive. Since this comes from the rotation of the pulsar, this tends to be very regular. After all, for a ball 10 km across, with more mass than the sun, smoothed to within a millimeter by its own gravity, it takes an awful lot to change how fast it's spinning. What's more, time is the quantity we can make the best measurements of - world time standards drift by something like microseconds over decades, which is something like one part in ten to the fourteen. So pulsar timing is a powerful technique, that can measure pulsar positions and distances, spin-down rates and braking indices, and binary orbits. The standard software has been tempo, which is written in FORTRAN and has certain limitations. A new tool has recently appeared, tempo2, written in C++ and boasting good handling of timing effects down to the nanosecond level. Unfortunately the documentation on this tool is so far somewhat limited, so I've been figuring out how it works. I'd like to describe it, as best I understand it, here. This particular post will talk about how a timing solution is fit to a set of pulse arrival times.

Millisecond pulsar X-rays

Neutron stars are so tiny and so dense that the natural temperature scale for them has them glowing in the X-rays. What's more, they serve as powerful accelerators of electrons, which then naturally produce X-rays. So it turns out that X-ray telescopes provide a very interesting view of pulsars. In light of this, when we discovered my pet source, J1023, we took an X-ray observation of it. It's taken us considerable time to analyze the results and put together a paper describing them, but I think the result will be a valuable contribution to the literature. (More importantly for me, it should make a chapter of my PhD thesis.) The result is the cumbersomely-named "X-ray Variability and Evidence for Pulsations from the Unique Radio Pulsar/X-ray Binary Transition Object FIRST J102347.6+003841".

Much of the paper is devoted to details of data analysis, which I will spare you. But I think the gist is interesting, and not too hard to summarize. 

Real interbinning

(This is another of my highly-technical "note to self" signal processing posts. I'll put up something less arcane soon.)


The Fourier transform is great for finding periodic signals: you take an FFT and a periodic signal looks like a peak in the output. Well, in an ideal world, that is; you only really get such a neat and tidy peak if the periodic signal is exactly at a Fourier frequency, which happens when it makes exactly an integral number of turns over the course of the data set. If the signal is somewhere between two Fourier frequencies, the power is spread over several Fourier output values. While it's possible to interpolate a very accurate value based on the ~32 nearest values, this can be expensive, and there's a shortcut called "interbinning"; it doesn't reconstruct the phase, but you just take a scaled average of the two neighbouring bins and get a decent approximation to the value at the midpoint between two independent Fourier frequencies.

My problem, for this post, is that the theory is all nicely worked out when going from the time domain to the frequency domain, but I want to do something analogous while going from the frequency domain back to the time domain, if that's possible. (I haven't done a literature search, or even looked very carefully at the frequency-domain interbinning papers; I thought this would be a good exercise for me.)

Glitches and flares

Recently on the arxiv: Searching for X-ray Variability in the Glitching Anomalous X-ray Pulsar 1E 1841-045 in Kes 73, by Zhu and Kaspi.

Pulsars normally spin down very regularly — like clockwork, as the saying goes, and many pulsars spin down as regularly as a good atomic clock. But some pulsars, once in a while, will suddenly start spinning more quickly. This sudden (instantaneous as far as we can measure) spin-up is called a "glitch", and its full explanation remains mysterious. Generally, all we see is that the pulsar is suddenly spinning faster: no heating of the crust, no sudden X-ray emission, no radiative changes at all, just a suddenly-faster pulsar.

Anomalous X-ray pulsars (AXPs) are one kind of "magnetar", pulsars whose magnetic field is so enormous that its decay powers the X-ray emission of the star. They exhibit many peculiar behaviours, and are a major field of study in pulsar research. AXPs will occasionally become much more active for a while: they become much brighter, they emit random blasts of X-rays, and they do other peculiar things. It seems as if they may glitch every time they become active like this. If we want to find some sort of relationship between glitches and these active periods, it would be valuable to know whether an active period happens every time AXPs glitch, or whether AXPs sometimes have "quiet" glitches, like normal pulsars. That's what this paper tries to answer.

RFI

I'm currently involved in another low-frequency pulsar survey with the GBT. This is a pointed survey, so we take two-minute pointings in an array tiling the whole northern part of the sky. This is an incredibly sensitive setup, so we will hopefully find a nice collection of new pulsars. Unfortunately, the same arrangement that makes us sensitive to pulsars makes us sensitive to lightning, electric motors, the local power lines, passing aircraft... all sorts of signals, generally outside our intended piece of sky and/or outside our intended frequency range, but also generally vastly more powerful than the pulsars we're looking for. So dealing with this junk — generically called RFI (Radio Frequency Interference) — is an important part of our survey strategy. My topic today is not so much dealing with the mild RFI we see all the time as dealing with what we see every now and then: a tremendous signal comes booming in, overwhelming everything and ruining the data (as seen on the right).

First of all a brief summary of how our system works: The telescope has a receiver at the primary focus, basically a pair of crossed dipoles. This converts the electromagnetic waves to signals in the wires, where they're amplified, encoded in analog on light and sent down an optical fiber to the control room, where it's converted down to a standard frequency range. This is fed to GUPPI, a so-called "pulsar backend". This instrument does analog-to-digital conversion on the signal then uses a polyphase filterbank to efficiently split it into 4096 channels, then computes the power in each channel; this power is averaged over 81.92 microseconds, after which the set of 4096 powers is dumped to disk.

Pulsars are broadband signals, so we use the widest bandwidth our receiver can handle, about 80 MHz centered near 350 MHz. This gives us the best signal-to-noise from our pulsars, but it also means we can't stick to the few narrow bands reserved for radio astronomy. In fact, the region we're using is allocated for various kinds of radio operation, including "aeronautical radionavigation" around 330 MHz. The GBT is in the National Radio Quiet Zone, so with luck much of this spectrum might be clear, but if aircraft are really using it to navigate there's no keeping them out of our airspace. (And at least once when I've complained of RFI the GBT operator has noticed an airplane passing to the north, for what that's worth.) In any case, there is interference to deal with.

Old radioastronomy systems, in order to get the widest possible bandwidth, quantized the input signal to one bit, or three levels, on initial digitization. This has the serious drawback that any small change in the signal levels causes the digitizer to overload or underflow, trashing the output. Newer backends improve on this by using 8-bit digitization, so that you have more than twenty decibels difference between underflow and overflow; it also helps that they are sampling a wider bandwidth, so that it takes more power to substantially change the levels. So generally, when a noise signal is picked up by GUPPI, it is narrowband and produces excess power in one channel. But the polyphase filterbank is good at keeping it from spilling into neigbouring channels (much better than a simple FFT) as long as it doesn't overflow the system, so usually you see something like this:

This shows, in red, the average power in each channel, and in blue, the standard deviation in each channel; since what we are receiving is basically noise, the standard deviation in each channel is large. The overall profile is the sensitivity of our system as a function of frequency. You can see it's not flat over the 80 MHz we want and zero elsewhere, but some complicated bumpy shape; that's the joy of analog components. In fact, given that it's an 80 MHz band at only 350 MHz, we're doing well to have a filter even that good. The point of this plot, though, is those vertical lines. Those are channels where some narrowband signal is coming in, so strongly that it stands way out above the noise. But it's okay; we have four thousand channels, we can afford to throw away a handful, and we have tools to detect when to throw them away (clever tools that look not just for strong signals but also for weaker signals that are narrowband and periodic, since those are a real headache for pulsar searching). So this plot is of one of our good beams.

The beams I'm worried about are the ones that look like this:

Here what's happened is we have some narrowband signal that is so monstrous it's overloading our digitizer. The digitizer's saturation is a nonlinear effect, so it doesn't stay politely in one channel; instead it splatters all over the place, producing both ugly broad peaks as on the right, and (I think) the regularly-spaced vertical lines (though they might be something else). When this happens there's not much use to be had from the data; while it's possible we might detect a pulsar through all that muck, on the one hand we'd be vastly less sensitive, while on the other if the interference has any periodic component to it, it'll produce zillions of false-positive periodicity candidates to wade through. So we need to discard beams like this one.

Fortunately, since we're doing a pointed survey, if we find that a beam has been trashed by RFI, we don't have to just write it off and leave a gap in our survey coverage: we can just put the beam back on our to-observe list. No problem. All I need to do is come up with a scheme for detecting these bad beams automatically, since we often run our survey observations semi-unattended (once they're up and running, I often give the operator my cell phone number and ask them to call me if something like a wind stow or a snow dump comes up; one of my colleagues will often do this and then go to sleep).

Initially, we were planning to run the RFI detection and excision scripts as part of our survey processing. Unfortunately, pulsar survey processing is incredibly CPU-intensive - we are just now finishing off the processing for the drift-scan survey, whose datataking finished in mid-2007. Clearly if we wait until the beams are processed to decide which are RFI-laden, the survey will long be over before we find out which beams we should have reobserved. Fortunately, the RFI detection code isn't too compute-intensive; an eight-processor modern machine running jobs on all processors ought to be able to process eight two-minute beams in about fifteen minutes.

The RFI detection code isn't perfect, though: I found that it was able to detect overloaded channels, but a few overloaded channels are normal. When too many are overloaded we can discard the beam, but I found that there were beams where not very many channels were actually overloaded, but those few splattered into neighboring channels. So my prototype system does an additional test: it compares the bandpass of each beam with the instrument's normal bandpass. This sort of works, but it revealed some surprising things:


This is a different plot of the same sort of thing. The red curve is the bandpass of the system, estimated by taking the median of all the day's beams. The green curve is one of the horrible contaminated beams. But the blue curve is from a beam that is probably fine. But its bandpass looks decidedly different (enough so to trick an early generation of my script into marking it bad). It turns out that what has happened is that this is a beam shortly after a contaminated beam. While the RFI was going on, we ran a system "balance", which adjusts the gain/attenuation at many points along the signal path so that in each stage it's at a comfortable level. Since there was powerful RFI, the gains and attenuations all got readjusted, and all differently. Since different parts of the signal chain affect the bandpass differently, we got a different bandpass. (We also got a drastically reduced signal amplitude, but fortunately GUPPI gives us plenty of dynamic range, as I mentioned above.) So I have to be a little careful when comparing bandpasses so as not to reject minor changes like this. Fortunately the "splatter" from bad RFI is pretty obvious in the statistics, so now I think I have a working bad beams detector. I'm just waiting for feedback from my collaborators before putting it into service.


Some comments on the tools I used to build it: I wrote the code in python and numpy, of course. But key to the process was a module by the author of the RFI detection tool that let me construct python objects from the detector's statistics files. Given this, I used medians to construct a bandpass for each file (since medians are better than means at discarding crazy outliers, which is the whole point of the exercise). I then used masked arrays to flag any bad channels, and scale the result so its median is one. I repeat this process on many beams, then take a median (using the masked array median to nicely ignore any bogus data points). This gives me my system bandpass. Comparing an individual beam to the bandpass then proceeds by constructing a masked median of the file as before, scaling it so its median matches that of the system bandpass, and then counting the points where it is very different.

In all, it has proved invaluable to have the masked array tools; they just do the Right Thing with bad data, vastly simplifying my code.

Thresholds

I apologize for three highly-technical posts in a row; I'm trying to work something out and setting it down "on paper" as it were is helping. I promise I'll come up with a post about kittens or something soon.

Suppose you're searching for pulsars. You're going through the Fourier transform looking for peaks. Now suppose you've found one: how strong is it? Is it statistically significant? For that matter, is it better than any of the list of peaks you're already keeping track of?

To answer this question you need two pieces of information: how strong the background noise is, and how likely that noise would have just randomly produced a peak this high. There's some cleverness in estimating the background, since real signals don't have perfectly flat white noise backgrounds, but I'm going to leave that aside for the moment. My question for today is, what are the statistics of a coherent peak-based search?

I think they can reasonably be estimated by ignoring any oversampling and assuming that, in the absence of a signal, if you're using n harmonics, the profile consists of 2n independent Gaussians with standard deviation 1. So the question is, for a fixed false positive probability p, what threshold should we set? The answer is roughly the threshold for a single Gaussian to exceed p/2n.

This is a little inconvenient to work with, since it requires the error function and its inverse (or an approximation), but flops are free. I suppose calling the error function code might cause cache misses, but I don't imagine needing it that often. Specifically, my idea is this: a simple implementation might just set a threshold at the beginning and run through the whole FFT. But if the observation underwent bad enough RFI, you could find yourself with millions of candidate signals. Since you'd be fine-tuning and storing every one, this could slow things down a lot. My idea is instead to specify a maximum number of candidates - generously, maybe a few hundred - and keep the best ones. This means raising the threshold every time you get a new candidate once the list is full. This wouldn't require the error function except that you don't only want to look at candidates with the full 64 harmonics - you also need to consider those with fewer. And converting thresholds between different numbers of harmonics does require the error function.

Coherent harmonic summing

As I alluded to in a previous post, you can (and it is sometimes useful to) take a giant FFT of a time series, extract a series of harmonically-related coefficients, and with an inverse FFT, produce a "pulse profile" showing the data "folded" at the period of the fundamental. This is interesting in part because you can use the peak height to gauge the strength of the signal, taking advantage of the relative phases of the harmonics. My question today is, when you're extracting those coefficients, how accurate do you need to be?

If all you want is the power, then a first approximation would be to simply choose the nearest Fourier frequency. This, after all, is where your FFT naturally measures the coefficient. I'll call these "independent Fourier frequencies", and the spacing between them the "independent Fourier spacing" (IFS). It turns out that if you have a frequency exactly halfway between two independent Fourier frequencies, you lose a substantial amount of sensitivity: a nearly 36% loss of signal-to-noise. You can cut this back to something like 7.4% using an ultra-simple interpolation scheme called "interbinning".

What about coherent reconstruction? Well, now we need not just the amplitude but the phase. The amplitude is relatively easy to get close to since it's at a maximum at the point of interest, so that the first derivative is zero and you have little dependence on the frequency. The phase does not have an extremum at the correct frequency, so it may well be varying rapidly. In fact, for a signal that is present uniformly throughout the observation, the phase changes by pi units per IFS. So if we use a spacing of IFS/2, which would have been adequate for the power, our phase will be wrong by as much as forty-five degrees.

What we care about, though, is not the phase exactly, but the real part. More, there are really two questions here: how well do we need to interpolate the FFT to get reasonably accurate Fourier coefficients, and how closely must we space our inverse FFTs?

First the easy question: how well do we need to interpolate to get decent Fourier coefficients. There are techniques for doing really good-quality interpolation in FFTs - you can use sinc interpolation (based on the 32 nearest samples), or you can just pad the time series before taking your giant FFT. But this is somewhat expensive. Since flops are free as long as they only access memory you've recently read anyway, to speed things up you can always linearly interpolate between more-carefully interpolated samples. So here's a plot of the error introduced by that linear interpolation:



So if you do proper Fourier interpolation to a spacing of IFS/8, and linear interpolation beyond that, then your measured amplitude will be off by less than about 3%. Using IFS/4 still leaves the worst-case error at about 7%, about the same as interbinning.

Now for the tougher question: how far off can our estimate of the frequency be when we reconstruct our profile? This is crucial, since it determines how many inverse FFTs we need to do, and these are probably the most expensive part of the calculation. We can safely think of this in the time domain: we have a narrowly-peaked pulsar, and we're folding the incoming data at a slightly wrong frequency. How much does this lower the peak amplitude?

Well, if the frequency is wrong, over the course of the observation the pulsar's peak will drift in phase. The average profile will then be smeared out by the amount of drift. Exactly how much it will lose in peak height depends on the pulse shape. An error in frequency of one IFS will result in one full cycle of drift (in fact this is what defines an IFS). So if we are going up to the nth harmonic, then an error in (fundamental) frequency of IFS/n results in one full turn of drift for that harmonic, and we might as well not have bothered with that harmonic. But focusing on an individual harmonic can answer the question.

The loss in amplitude of a harmonic whose frequency is off by x is given by the sinc function; if we approximate the sinc function with its quadratic Taylor polynomial, we get a loss equal to x²π²/6, where x is the frequency error in IFS. Now, if we suppose that the profile is effectively a delta function, so that its n coefficients are all 1/n, then the total error is the sum over m of m²x²π²/6n or x²n(n+1)(2n+1)π²/36n.

What this works out to, once you clear up the messy math, is that if you sample the top harmonic at IFS/2, you lose about 14%; if you sample at IFS/4 you lose about 3.6%, and if the odd number doesn't make you queasy and you sample at IFS/3 you lose about 6.3%. (It turns out the number of harmonics is nearly irrelevant.)

So, in short, if you interpolate the FFT to IFS/8 and linearly interpolate beyond that, and you take an FFT every time you advance the top harmonic by its IFS/4. you'll lose no more than about 5% sensitivity at worst. If you cut your number of FFTs in half (which probably cuts your runtime in half), you lose at worst maybe 20% of your sensitivity, but you can probably make it up by setting a threshold 20% lower, then "tuning up" each candidate to find the best period.

Flops

In a recent astronomy talk, the speaker was discussing some sort of heavy-duty calculations, and how to make them faster. The way he put it, ultimately, was "flops are free". That is, CPUs are so fast now that they spend almost all their time waiting for data to be fetched from main memory: this means that the time it takes for a job to run is determined not by how many floating-point calculations it has to run but by how much data it must read from main memory.

This is kind of an astonishing statement, really: it says you can add floating-point calculations to your code at no extra cost. Of course such a statement can only be true for certain computing tasks on certain platforms, so I thought I'd test how true it is for a particular computing task I had in mind.

The task, you will probably not be surprised to discover, is pulsar searching. Specifically, I want to think about going through the Fourier transform of an observation looking for periodic signals.

Traditionally the way to do this is to take a giant FFT, square it to extract the power, and look for statistically-significant peaks. In fact, since pulsar signals are typically not just sine waves, one normally adds up the power from several harmonics and looks for statistically-significant peaks in the total power. There are lots of important details, but let's leave those aside for the moment. The idea I've been thinking about for a while now is based on a distinction.

When you estimate the total power in the harmonics, by a beautiful theorem in harmonic analysis, you are effectively estimating the root-mean-squared amplitude of the pulsar's folded pulse profile. But if you look at a profile that has one narrow peak, as many pulsars do, the peak-minus-mean amplitude can be much much more significant than the root-mean-squared amplitude. Back in the Fourier domain, when you have a single peak, not only is there power in many harmonics, but those harmonics line up in phase. The usual approach ignores the phases of the harmonics entirely. How much could be gained by paying attention to them?

Some informal experiments suggest that for a pulsar that is on only 1% of the time (which is not too rare among slow pulsars), this approach may offer something like a 40% improvement in sensitivity. Since that's the same as doubling the observation time, I think it's worth looking into.

As with everything, this improved sensitivity comes at a cost. In particular, all the codes we have to compute it are substantially slower than the code we have that uses the incoherent approach. So, thinking about how to speed things up, it occurred to me to look into just why the code was slow.

I'm sure there are brilliant and subtle tools to count cache misses and do nanosecond timing on running code, but I thought I'd take a more direct approach: just write dummy code that does only one thing and time it. In particular, what I want to compare is the time to extract 64 harmonics from a giant FFT to the time it takes to take an inverse FFT of those harmonics to reconstruct the profile.

I should say that my initial feeling was that the inverse FFT would dominate the time - after all, even an n log n algorithm requires a fair number of floating-point operations to produce a 256-point real inverse FFT. But the output is:

Array reading 1.71528 GB/s
Planning FFT
Running FFTs
FFT time 1.57851e-06 seconds
FFT rate 633510 /s
FFT rate/read rate: 0.176112

That is, yes, the FFTs are the slow step, but only by a factor of five. That is, reading all those coefficients in from memory takes a fifth as long as doing all those FFTs. (Incidentally, if these numbers seem low, the machine I'm running this on is a couple of years old. Edit: it's not the machine pictured above, which is even older, and which does most of my number-crunching.) Here's the code:

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#include <complex.h>
#include <fftw3.h>

const int repeats=20;
const int array_size=1<<26;
const int irfft_size=256;
const int fft_batch=1<<18;

int main(int argc, char*argv[]) {
  double t, tnew;
  float*array;
  double sum;
  double array_rate;
  int i,j,k;
  struct timeval start, end;
  fftwf_plan plan;
  fftwf_complex *in;
  float *out;

  array = (float*)malloc(array_size*sizeof(float));
  if (!array) {
      perror("Allocation failed");
      return 1;
  }
  sum = 0;

  t = 1e100;
  for (j=0;j<repeats;j++) {
      gettimeofday(&start,0);
      for (i=0;i<array_size;i++)
          sum += array[i];
      gettimeofday(&end,0);
      tnew = end.tv_sec-start.tv_sec+(end.tv_usec-start.tv_usec)/1e6;
      if (tnew<t) t=tnew;
  }
  printf("Array reading time %g seconds\n", t, sum);
  array_rate = array_size*sizeof(float)/t;
  printf("Array reading %g GB/s\n", array_rate/(1<<30));
  free(array);

  in = (fftwf_complex*) fftwf_malloc(sizeof(fftwf_complex)*(irfft_size/2+1));
  out = (float*) fftwf_malloc(sizeof(float)*irfft_size);
  printf("Planning FFT\n");
  plan = fftwf_plan_dft_c2r_1d(irfft_size, in, out, FFTW_MEASURE | FFTW_EXHAUSTIVE | FFTW_DESTROY_INPUT);
  printf("Running FFTs\n");

  t = 1e100;
  for (j=0;j<repeats;j++) {
      gettimeofday(&start,0);
      for (i=0;i<fft_batch;i++) {
          for (k=0;k<(irfft_size/2+1);k++) {
              in[k] = 0;
          }
          fftwf_execute(plan);
      }
      gettimeofday(&end,0);
      tnew = end.tv_sec-start.tv_sec+(end.tv_usec-start.tv_usec)/1e6;
      if (tnew<t) t=tnew;
  }
  printf("FFT time %g seconds\n", t/fft_batch);
  printf("FFT rate %g /s\n", fft_batch/t, sum);

  printf("FFT rate/read rate: %g\n", (fft_batch/t)/(array_rate/(2*sizeof(float)*(irfft_size/4))));

  fftwf_destroy_plan(plan);
  fftwf_free(in);
  fftwf_free(out);
  return 0;
}

Compiled with:

gcc -O9 -march=native -ffast-math -lfftw3f -lm timer.c

What this tells me is interesting. On the one hand, if I can do something really clever to make the FFTs faster, or avoid them, I can get some improvement, but sooner or later I'm going to hit the limit of memory loading. On the other hand, if I can't make them faster or fewer - and I probably can't - there's not much point sweating much over reducing the memory loads.

Anyway, the upshot of all this is: on modern hardware, you can do an awful lot of flops - a 256-point real FFT - for not much more than the cost of loading the data in from main memory. So if you have some clever mathematical trick to reduce your data size (interpolation, say) it may be worth implementing.

Plots

I have been doing some X-ray astronomy. In optical astronomy, spectroscopy is a very powerful tool: by looking for emission and absorption lines you can identify the elements present in a gas (helium was discovered this way, for example); the shapes of the lines can tell you about temperatures and velocities in the object, and the shape of the broadband spectrum can also tell you about the temperature and conditions in the emission region. In X-rays, things are more difficult, for a number of reasons. Unfortunately, lines are much rarer (at least when looking at neutron stars), telescope time is very scarce (since the telescopes must be in space), and there's always a shortage of photons. But X-ray spectroscopy still has the potential to tell you about temperatures, sizes, and compositions of neutron stars (for example). So that's what I've been working on.

The standard tool for X-ray spectroscopy is xspec, one of those pieces of scientific software that's had a great deal of cleverness built into it, very little of which has gone into making it easy to use. It could be worse - at least its interface is not stuck in the FORTRAN era, in fact it has a tcl interpreter built in (yack) - but its plotting in particular is pretty rudimentary, tending to produce monstrosities like this:



The worst part of this graph, apart from the fact that it's practically unreadable even for those with normal color vision, is that it's quite deceptive. It looks as if there's clear evidence for a bend in the spectrum just below 2 keV. But look at the plot below the jump for comparison.





The data points on this plot are identical, but drawing a single power-law through the whole thing makes it look like the data's completely straight. Ordinarily one would choose between the models based on the statistics, but we have so few photons (about 8000) that both models are perfectly adequate fits to the data. I suppose Occam's Razor tells me I should pick the simpler model, though whether this should be the simple but not particularly physical power-law or the physically plausible but more complicated power-law plus neutron-star atmosphere model isn't entirely clear to me.

I'll keep thinking about how to improve the graphics, but the problem is I have five data sets in which each data point has its own vertical and horizontal error bars, and the model gives slightly different predictions for each data set (since they use different instruments with slightly different responses). The plotting tools provided by xspec are also not very flexible (and I haven't found a good way to export the relevant data so I can use my own plotting tools).

What really bothered me, though, was how strongly each graph suggests its own interpretation. It would have been easy to look at one and assume that it told the whole story.

Light bending

Whether or not General Relativity is the correct theory of gravity on very large scales, it has passed all tests (many pulsar-based) when applied to planetary and solar system scales. One important feature of the theory is that in a gravitational field, light follows curved trajectories (technically geodesics are "as straight as possible", but they are curved in the practical sense). In familiar settings, this means tiny but detectable effects in laboratory experiments, or small but measurable deflection of stars near the Sun. But it turns out that pulsars are small enough and massive enough that light near their surfaces is bent a great deal - enough that you can actually see almost all of the surface of the pulsar at the same time.

This may sound bizarre, and to some degree it is, but it produces potentially measurable effects. A pulsar's radio emission is produced by plasma somewhere in its magnetosphere, and in fact we're not at all sure just where the emitting plasma is. But for many pulsars, their X-ray emission comes from "hot spots" on the surface, at the magnetic poles. For young pulsars, these hot spots arise because the magnetic field in the crust makes it much easier for heat to flow out where the field is vertical than where it's at an angle. For the very old millisecond pulsars, these hot spots arise from gigantic sparks in the magnetosphere blasting the surface with high-energy particles, heating it. In either case, we sometimes see X-ray pulsations with a thermal spectrum, and light-bending can explain some of the properties of these pulsations.

I did a few very simple simulations of the light-bending, and made some illustrative videos.


I made a set of three videos illustrating this effect. For a rotating pulsar, there are a number of geometric parameters, including the angle between the line of sight and the rotation axis, the angle between the rotation axis and the magnetic axis, and the size of the hot polar cap. I fixed values for all these. There is also the question of the physical size of the pulsar: the more compact and dense it is, the more light-bending we will see. I have generated three videos. The first shows the geometry with no light bending:



Below the actual animation I have a little plot showing a pulse profile, based on a very simple model (blackbody emission from the polar cap). If I make the pulsar more compact (R=3M), I get:



And more compact still (R=2.1M):



The size of each pulsar model is given in "geometrized units", where R=2M is the size of a black hole, the most extreme possible light bending. For a 1.4 solar mass neutron star, this is 4.1 km. Realistic neutron star models vary quite a lot in radius, from ~6 km to ~24 km (~3M to ~12M), so light bending will probably not be as strong as the neutron stars depicted here.

There are other effects to consider as well; neutron stars are expected to have atmospheres, and in fact their spectra do not look like black-body spectra. The atmospheres affect these results by "limb darkening", that is, the radiation that emerges is directed more vertically than simple black-body radiation, so these pulse profiles are not really right. This can be done better, but I just wanted to write a quick hack (using the usual suspects, python, numpy, scipy, and matplotlib) and get a feel for the effect.