Category: data-analysis

OK, so we're done with this epic of climate data analysis. I've prepared an index of the articles in this series, on the off chance that it might be useful for someone. The goal of this exercise was mostly to try doing some "basic" climate data analysis tasks in Haskell, things…

This is going to be the last substantive post of this series (which is probably as much of a relief to you as it is to me...). In this article, we're going to look at phase space partitioning for our dimension-reduced PCA data and we're going to calculate Markov transition matrices for our partitions to try to pick out consistent non-diffusive transitions in atmospheric flow regimes.

(There's no code in this post, just some examples to explain what we're going to do next.)

Suppose we define the state of the system whose evolution we want to study by a probability vector -- at any moment in time, we have a probability distribution over a finite partition of the state space of the system (so that if we partition the state space into components, then ). Evolution of the system as a Markov chain is then defined by the evolution rule

where is a Markov matrix. This approach to modelling the evolution of probability densities has the benefit both of being simple to understand and to implement (in terms of estimating the matrix from data) and, as we'll see, of allowing us to distinguish between random "diffusive" evolution and conservative "non-diffusive" dynamics.

We'll see how this works by examining a very simple example.

The analysis of preferred flow regimes in the previous article is all very well, and in its way quite illuminating, but it was an entirely static analysis -- we didn't make any use of the fact that the original data we used was a time series, so we couldn't gain any information about transitions between different states of atmospheric flow. We'll attempt to remedy that situation now.

What sort of approach can we use to look at the dynamics of changes in patterns of ? Our parameterisation of flow patterns seems like a good start, but we need some way to model transitions between different flow states, i.e. between different points on the sphere. Each of our original
maps corresponds to a point on this sphere, so we might hope that we can some up with a way of looking at trajectories of points in space that will give us some insight into the dynamics of atmospheric flow.

A quick post today to round off the "static" part of our atmospheric flow analysis.

Now that we've satisfied ourselves that the bumps in the spherical PDF in article 8 of this series are significant (in the narrowly defined sense of the word "significant" that we've discussed), we might ask what to sort of atmospheric flow regimes these bumps correspond. Since each point on our unit sphere is really a point in the three-dimensional space spanned by the first three PCA eigenpatterns that we calculated earlier, we can construct composite maps to look at the spatial patterns of flow for each bump just by combining the first three PCA eigenpatterns in proportions given by the " " coordinates of points on the unit sphere.

The spherical PDF we constructed by kernel density estimation in the article before last appeared to have "bumps", i.e. it's not uniform in and . We'd like to interpret these bumps as preferred regimes of atmospheric flow, but before we do that, we need to decide whether these bumps are significant. There is a huge amount of confusion that surrounds this idea of significance, mostly caused by blind use of "standard recipes" in common data analysis cases. Here, we have some data analysis that's anything but standard, and that will rather paradoxically make it much easier to understand what we really mean by significance.

The Haskell kernel density estimation code in the last article does work, but it's distressingly slow. Timing with the Unix time command (not all that accurate, but it gives a good idea of orders of magnitude) reveals that this program takes about 6.3 seconds to run. For a one-off, that's not too bad, but in the next article, we're going to want to run this type of KDE calculation thousands of times, in order to generate empirical distributions of null hypothesis PDF values for significance testing. So we need something faster.

Up to this point, all the analysis that we've done has been what might be called "normal", or "pedestrian" (or even "boring"). In climate data analysis, you almost always need to do some sort of spatial and temporal subsetting and you very often do some sort of anomaly processing. And everyone does PCA! So there's not really been anything to get excited about yet.

Now that we have our PCA-transformed anomalies though, we can start to do some more interesting things. In this article, we're going to look at how we can use the new representation of atmospheric flow patterns offered by the PCA eigenpatterns to reduce the dimensionality of our data, making it much easier to handle. We'll then look at our data in an interesting geometrical way that allows us to focus on the patterns of flow while ignoring the strengths of different flows, i.e. we'll be treating strong and weak blocking events as being the same, and strong and weak "normal" flow patterns as being the same. This simplification of things will allow us to do some statistics with our data to get an idea of whether there are statistically significant (in a sense we'll define) flow patterns visible in our data.

Although the basics of the "project onto eigenvectors of the covariance matrix" prescription do hold just the same in the case of spatio-temporal data as in the simple two-dimensional example we looked at in the earlier article, there are a number of things we need to think about when we come to look at PCA for spatio-temporal data. Specifically, we need to think bout data organisation, the interpretation of the output of the PCA calculation, and the interpretation of PCA as a change of basis in a spatio-temporal setting. Let's start by looking at data organisation.

The pre-processing that we've done hasn't really got us anywhere in terms of the main analysis we want to do -- it's just organised the data a little and removed the main source of variability (the seasonal cycle) that we're not interested in. Although we've subsetted the original geopotential height data both spatially and temporally, there is still a lot of data: 66 years of 181-day winters, each day of which has values. This is a very common situation to find yourself in if you're dealing with climate, meteorological, oceanographic or remote sensing data. One approach to this glut of data is something called dimensionality reduction, a term that refers to a range of techniques for extracting "interesting" or "important" patterns from data so that we can then talk about the data in terms of how strong these patterns are instead of what data values we have at each point in space and time.

I've put the words "interesting" and "important" in quotes here because what's interesting or important is up to us to define, and determines the dimensionality reduction method we use. Here, we're going to side-step the question of determining what's interesting or important by using the de facto default dimensionality reduction method, principal components analysis (PCA). We'll take a look in detail at what kind of "interesting" and "important" PCA give us a little later.

Note: there are a couple of earlier articles that I didn't tag as "haskell" so they didn't appear in Planet Haskell. They don't contain any Haskell code, but they cover some background material that's useful to know (#3 talks about reanalysis data and what is, and #4 displays some of the characteristics of the data we're going to be using). If you find terms here that are unfamiliar, they might be explained in one of these earlier articles.

The code for this post is available in a Gist.

Update: I missed a bit out of the pre-processing calculation here first time round. I've updated this post to reflect this now. Specifically, I forgot to do the running mean smoothing of the mean annual cycle in the anomaly calculation -- doesn't make much difference to the final results, but it's worth doing just for the data manipulation practice...

Before we can get into the "main analysis", we need to do some pre-processing of the data. In particular, we are interested in large-scale spatial structures, so we want to subsample the data spatially. We are also going to look only at the Northern Hemisphere winter, so we need to extract temporal subsets for each winter season. (The reason for this is that winter is the season where we see the most interesting changes between persistent flow regimes. And we look at the Northern Hemisphere because it's where more people live, so it's more familiar to more people.) Finally, we want to look at variability about the seasonal cycle, so we are going to calculate "anomalies" around the seasonal cycle.

We'll do the spatial and temporal subsetting as one pre-processing step and then do the anomaly calculation seperately, just for simplicity.

In the last article, I talked a little about geopotential height and the data we're going to use for this analysis. Earlier, I talked about how to read data from the NetCDF files that the NCEP reanalysis data comes in. Now we're going to take a look at some of the features in the data set to get some idea of what we might see in our analysis. In order to do this, we're going to have to produce some plots. As I've said before, I tend not to be very dogmatic about what software to use for plotting -- for simple things (scatter plots, line plots, and so on) there are lots of tools that will do the job (including some Haskell tools, like the Chart library), but for more complex things, it tends to be much more efficient to use specialised tools. For example, for 3-D plotting, something like Paraview or Mayavi is a good choice. Here, we're mostly going to be looking at geospatial data, i.e. maps, and for this there aren't really any good Haskell tools. Instead, we're going to use something called NCL (NCAR Command Language). This isn't by any stretch of the imagination a pretty language from a computer science point of view, but it has a lot of specialised features for plotting climate and meteorological data and is pretty perfect for the needs of this task (the sea level pressure and plots in the last post were made using NCL). I'm not going to talk about the NCL scripts used to produce the plots here, but I might write about NCL a bit more later since it's a very good tool for this sort of thing.

In this article, we're going to look at some of the details of the data that we're going to be using in our study of non-diffusive flow in the atmosphere. This is still all background material, so there's no Haskell code here!

As I said in the last article, the next bit of this data analysis series is going to attempt to use Haskell to reproduce the analysis in the paper: D. T. Crommelin (2004). Observed nondiffusive dynamics in large-scale atmospheric flow. J. Atmos. Sci. 61(19), 2384--2396. Before we can do this, we need to cover some background, which I'm going to do in this and the next couple of articles. There won't be any Haskell code in any of these three articles, so I'm not tagging them as "Haskell" so that they don't end up on Planet Haskell, annoying category theorists who have no interest in atmospheric dynamics. I'll refer to these background articles from the later "codey" articles as needed.

I never really intended the FFT stuff to go on for as long as it did, since that sort of thing wasn't really what I was planning as the focus for this Data Analysis in Haskell series. The FFT was intended primarily as a "warm-up" exercise. After fourteen blog articles and about 10,000 words, everyone ought to be sufficiently warmed up now...

Instead of trying to lay out any kind of fundamental principles for data analysis before we get going, I'm just going to dive into a real example. I'll talk about generalities as we go along when we have some context in which to place them.

All of the analysis described in this next series of articles closely follows that in the paper: D. T. Crommelin (2004). Observed nondiffusive dynamics in large-scale atmospheric flow. J. Atmos. Sci. 61(19), 2384--2396. We're going to replicate most of the data analysis and visualisation from this paper, maybe adding a few interesting extras towards the end.

It's going to take a couple of articles to lay out some of the background to this problem, but I want to start here with something very practical and not specific to this particular problem. We're going to look at how to gain access to meteorological and climate data stored in the NetCDF file format from Haskell. This will be useful not only for the low-frequency atmospheric variability problem we're going to look at, but for other things in the future too.

After releasing the arb-fft package last week, a couple of issues were pointed out by Daniel Díaz, Pedro Magalhães and Carter Schonwald. I've made another couple of releases to fix these things. They were mostly little Haskell things that I didn't know about, so I'm going to quickly describe them here. (I'm sure these are things that most people know about already...)

After thirteen blog articles and about two-and-a-half weeks' equivalent work time spread out over three months, I think it's time to stop with the FFT stuff for a while. I told myself when I started all this that if I could get within a factor of 20 of the performance of FFTW, I would be "quite obscenely pleased with myself". For most , the performance of our FFT code is within a factor of 5 or so of the performance of the vector-fftw package. For a pure Haskell implementation with some but not lots of optimisation work, that's not too shabby.

But now, my notes are up to nearly 90 pages and, while there are definitely things left to do, I'd like to move on for a bit and do some other data analysis tasks. This post is just to give a quick round-up of what we've done, and to lay out the things remaining to be done (some of which I plan to do at some point and some others that I'd be very happy for other people to do!). I've released the latest version of the code to Hackage. The package is called arb-fft. There is also an index of all the posts in this series.

In this article, we'll finish up with optimisation, tidying up a couple of little things and look at some "final" benchmarking results. The code described here is the pre-release-4 version in the GitHub repository. I'm not going to show any code snippets in this section, for two reasons: first, if you've been following along, you should be familiar enough with how things work to find the sections of interest; and second, some of the optimisation means that the code is getting a little more unwieldy, and it's getting harder to display what's going on in small snippets. Look in the GitHub repository if you want to see the details.

In the last article, we did some basic optimisation of our FFT code. In this article, we're going to look at ways of reordering the recursive Cooley-Tukey FFT decomposition to make things more efficient. In order to make this work well, we're going to need more straight line transform "codelets". We'll start by looking at our example in detail, then we'll develop a general approach.

Based on what we saw concerning the performance of our FFT code in the last article, we have a number of avenues of optimisation to explore. Now that we've got reasonable looking scaling for all input sizes, we're going to try to make the basic Danielson-Lanczos step part of the algorithm faster, since this will provide benefits for all input sizes. We can do this by looking at the performance for a single input length (we'll use ). We'll follow the usual approach of profiling to find parts of the code to concentrate on, modifying the code, then looking at benchmarks to see if we've made a positive difference.

Once we've got some way with this "normal" kind of optimisation, there are some algorithm-specific things we can do: we can include more hard-coded base transforms for one thing, but we can also try to determine empirically what the best decomposition of our input vector length is -- for example, for , we could decompose as , using length-2 base transforms and seven Danielson-Lanczos steps to form the final transform, or as , using a length-16 base transform and a single Danielson-Lanczos step to form the final result.

So far in this series, we've been looking at code piecemeal, writing little test modules to explore the algorithms we've been developing. Before we start trying to optimise this code, it makes sense to put it into a Cabal package with a more organised module structure, to provide a sensible API and to make a few other small changes.

From now on, the code we're going to be talking about will be the arb-fft package, which is hosted on GitHub. In this article, we'll be talking about the version of the code tagged pre-release-1. (We'll be uploading a version of the code to Hackage once it's ready, but there will be a few pre-release versions that we'll look at before we get to that point.)

The code from this article is available as a Gist

In the last article in this series, benchmarking results gave a pretty good indication that we needed to do something about prime-length FFTs if we want to get good scaling for all input sizes. Falling back on the DFT algorithm just isn't good enough.

In this article, we're going to look at an implementation of Rader's FFT algorithm for prime-length inputs. Some aspects of this algorithm are not completely straightforward to understand since they rely on a couple of results from number theory and group theory. We'll try to use small examples to motivate what's going on.

We're going to do a number of different things to optimise the code developed in the previous article, but before we do that, we need to do some baseline performance measurements, so that we know that our "improvements" really are improving things. The performance results we see will also help to guide us a little in the optimisations that we need to do.

Benchmarking is a fundamental activity when working with numerical code, but it can be quite difficult to do it correctly in a lazy language like Haskell. It's very easy to get into a situation where a benchmarking loop (perhaps to run a computation 1000 times to collect accurate timing information) performs the computation on the first time round the loop, then reuses the computed result for all the other 999 times round.

We can bypass these worries quite effectively using Bryan O'Sullivan's Criterion library. As well as providing a clean framework for ensuring that pure calculations really do get rerun for timing purposes, Criterion gives us a nice API for setting up benchmarks, running them and collecting results.

Let's summarise what we've done so far, as well as pointing out some properties of our general mixed-radix FFT algorithm that will be important going forward. We'll also think a little bit about just what we are going to do next. Here's what we've done:

1. Written down the basic expression for the DFT and implemented it directly in Haskell. POST
2. Built some "toy algebra" tools to help us understand the Fourier matrix decomposition at the heart of the Cooley-Tukey FFT algorithm. POST
3. Implemented a basic powers-of-two FFT. POST
4. Extended the "toy algebra" code to the case of general input vector lengths. POST
5. Implemented a full mixed-radix FFT in Haskell for general input vector lengths. POST

The code from this article is available as a Gist

In this article, we're going to implement the full "sorted prime factorisation mixed-radix decimation-in-time Cooley-Tukey FFT". That's a bit of a mouthful, and you might want to review some of the algebra in the last article to remind yourself of some of the details of how it works.

The approach we're going to take to exploring the implementation of this mixed-radix FFT is to start with a couple of the smaller and simpler parts (prime factor decomposition, digit reversal ordering), then to look at the top-level driver function, and finally to look at the generalised Danielson-Lanczos step, which is the most complicated part of things. This code has an unavoidably large number of moving parts (which is probably why this algorithm is rarely covered in detail in textbooks!), so we'll finish up by writing some QuickCheck properties to make sure that everything works .

The code from this article is available as a Gist

So far, we've been dealing only with input vectors whose lengths are powers of two. This means that we can use the Danielson-Lanczos result repeatedly to divide our input data into smaller and smaller vectors, eventually reaching vectors of length one, for which the discrete Fourier transform is just the identity. Then, as long as we take account of the bit reversal ordering due to the repeated even/odd splitting we've done, we can use the Danielson-Lanczos lemma or its matrix equivalent to efficiently reassemble the entries in our data vector to form the DFT (at this point, you may want to review the second article to remind yourself how this matrix decomposition works).

This is all standard stuff that you'll find in a lot of textbooks. Starting in this article though, we're going to veer "off piste" and start doing something a bit more interesting. It used to be that, if you wanted to do FFTs of vectors whose lengths weren't powers of two, the professional advice was "don't do that" (that's what my copy of Numerical Recipes says!). However, if you're prepared to brave a little bit of algebra, extending the Cooley-Tukey approach to vectors of arbitrary lengths is not conceptually difficult. It is a bit fiddly to get everything just right, so in this article, we'll start exploring the problem with some more "toy algebra".

The code from this article is available as a Gist.

Before we go on to look at how to deal with arbitrary input vector lengths (i.e. not just powers of two), let's try a simple application of the powers-of-two FFT from the last article.

We're going to do some simple frquency-domain filtering of audio data. We'll start with reading WAV files. We can use the Haskell WAVE package to do this -- it provides functions to read and write WAV files and has data structures for representing the WAV file header information and samples.

We'll use a sort of sliding window FFT to process the samples from each channel in the WAV file. Here's how this works:

We split the audio samples up into fixed-length "windows" (each
samples long) with adjacent windows overlapping by samples. We select to be a power of 2 so that we can use our power-of-two FFT algorithm from the last section. We calculate the FFT of each window, and apply a filter to the spectral components (just by multiplying each spectral component by a frequency-dependent number). Then we do an inverse FFT, cut off the overlap regions and reassemble the modified windows to give a filtered signal.

I'm not going to talk much about the details of audio filtering here. Suffice it to say that what we're doing isn't all that clever and is mostly just for a little demonstration -- we'll definitely be able to change the sound of our audio, but this isn't a particularly good approach to denoising speech. There are far better approaches to that, and I may talk about them at some point in the future.

The code from this article is available as a Gist. I'm going to go on putting the code for these articles up as individual Gists as long as we're doing experimentation -- eventually, there will be a proper repo and a Hackage package for the full "bells and whistles" code.

The decomposition of the Fourier matrix described in the previous article allows us to represent the DFT of a vector of length by two DFTs of length . If is a power of two, we can repeat this decomposition until we reach DFTs of length 1, which are just an identity transform. We can then use the decomposition

to build the final result up from these decompositions.

At each step in this approach, we treat the even and odd indexed elements of our input vectors separately. If we think of the indexes as binary numbers, at the first step we decompose based on the lowest order bit in the index, at the next step we decompose based on the next lowest bit and so on.

The code from this article is available as a Gist.

In the first article in this series, we considered a set of data values from a function sampled at a regular interval :

and defined the discrete Fourier transform of this set of data values as

The original Cooley and Tukey fast Fourier transform algorithm is based on the observation that, for even , the DFT calculation can be broken down into two subsidiary DFTs, one for elements of the input vector with even indexes and one for elements with odd indexes:

where we write , and where is the th component of the DFT of the evenly indexed elements of and is the th component of the DFT of the oddly indexed elements of . This decomposition is sometimes called the Danielson-Lanczos lemma.

In this article, we're going to look in detail at this decomposition to understand just how it works and how it will help us to calculate the DFT in an efficient way. In the next article, we'll see how to put a series of these even/odd decomposition steps together to form a full DFT.

This is the first of a series of posts about implementing the Fast Fourier Transform in Haskell (and trying to make it go fast). This part is going to be pretty pedestrian, but we need to lay some groundwork before the interesting stuff.

Fourier transforms turn up everywhere. They're important in pure mathematics, but they're also used in pretty much any application that deals with time series: filtering, signal analysis, numerical solution of differential equations and many others. While the classical theory of Fourier analysis is based on functions, in applications we most often deal with discrete series of values sampled from a function. The discrete Fourier transform is how we do Fourier analysis in this setting.

In this series of posts, we're not going to look very much at applications (well, maybe we'll have one little exception a bit later) and we're not going to talk a lot about the theory behind the Fourier transform. Instead, after briefly introducing the discrete Fourier transform, we're going to use the fast Fourier transform algorithm and some variations to explore some aspects of Haskell programming. In particular, we'll look at how Haskell deals with complex numbers, vectors, profiling and benchmarking, meta-programming and finally, as our pièce de résistance, we'll use all of this to build a Haskell fast Fourier transform package that does compile-time empirical optimisation for arbitrary-sized transforms.

In practice, if you want to do Fourier transforms of discrete data, you use something called FFTW ("the Fastest Fourier Transform in the West"). Our code has no chance of competing with FFTW, which is incredibly clever , but we can demonstrate some of the methods that they use and learn some techniques along the way that are applicable to more complicated problems. (There are Haskell bindings for FFTW -- the vector-fftw package is a good choice.)

I'm in the process of starting up a blogging series, to write about Data Analysis in Haskell. This is a sort of introduction/manifesto...