Category: haskell

One of the things that C2HS is lacking is a good tutorial. So I'm going to write one (or try to, anyway).

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.

I took over the day-to-day support for C2HS about 18 months ago and have now finally cleaned up all the issues on the GitHub issue tracker. It took a lot longer than I was expecting, mostly due to pesky "real work" getting in the way. Now seems like a good time to announce the 0.25.1 "Snowmelt" release of C2HS and to summarise some of the more interesting new C2HS features.

This is going to be the oldest of old hat for the cool Haskell kids who invent existential higher-kinded polymorphic whatsits before breakfast, but it amused me, and it's the first time I've used some of these more interesting language extensions for something "real".

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.

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.

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.

A useful but little used feature of Haddock is the ability to include inline images in Haddock pages. Here are a few examples. You can use images for diagrams or for inserting mathematics into your documentation in a readable way. In order to use an image, you just put <<path-to-image>> in a Haddock comment. The image can be of any format that browsers support: PNG, JPEG, SVG, whatever.

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 .

For my blog, I use Jasper Van der Jeugt's Hakyll system. When I was first looking for a blogging platform, I rejected all the usual choices (Wordpress, etc.), mostly because I'm a borderline obsessive control freak and they didn't give enough configurability. Plus, they weren't Haskell. I started writing my own blogging system, but then I found Hakyll. "Great," I thought, "that's perfect, except I want something almost completely different!"

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...

I got hold of a BeagleBone Black a few weeks ago (courtesy of Tom Nielsen, who'd had it for a while, but had no time to play with it). This is a small (credit card sized) Linux machine with an ARM processor, intended for, well, pretty much anything you can use a little computer for. It has 512 Mb of RAM, plus 2 Gb of flash storage, so it's not a completely trivial machine.

Obviously, to do anything really significant with it requires some hardware work (there are lots of general purpose I/O pins to play with, plus UARTs, Ethernet, USB, SPI, a couple of analogue-to-digital converters, PWM drivers for motor control, and even HDMI video output!), but there are a few LEDs on the board itself that lend themselves to a blinkenlights demo...

The BeagleBone comes with some built-in software to let people write code quickly, using JavaScript of all things. There's also a bit of documentation about programming the thing in C. But (of course) I wanted Haskell blinkenlights!

So, the goal was to write a little bit of code to count in binary on the four user-addressable LEDs on the board, first in Javascript, then in C, then (somehow) in Haskell.

I've recently started helping with the maintenance of C2HS, a tool for generating Haskell foreign function interface (FFI) bindings from C header files. I started using C2HS because of a Haskell library I was writing to read Unidata NetCDF files. I didn't fancy writing all of the bindings to the hundreds of functions in the C NetCDF library by hand and C2HS seemed like a good way to get started with the Haskell FFI.

For BayesHive, one of the things we need to do a lot of is program transformation. The main reason for this is to generate code to represent probability distribution functions derived from descriptions of probability models -- the models are described using code written in the probability monad, and this needs to be manipulated quite extensively to determine the PDF, which is needed for the Markov chain Monte Carlo sampling we use for estimating parameters in probability models.

Program transformation is something that Haskell is really great at, since it basically just involves pattern matching against lots of different cases of syntax trees in the language you're transforming. Unfortunately, this leads to lots of code that looks like this (picking a simple case at random):

...
case last rvs of
   Observed (EApp (EApp (EVar "map") (EVar f)) (EVar nm)) -> rewriteFst f nm
   _     -> rvs
...

All that stuff with EApp and EVar is the explicit AST representation of the expression map f nm in the language (Baysig) that we're manipulating.

At first, doing things like this is pretty unavoidable, but wouldn't it be nice if you could write something like this instead?

...
case last rvs of
   Observed [baysig|map  nm|] -> rewriteFst f nm
   _     -> rvs
...

For BayesHive, we have a fairly big Yesod web application plus a lot of client-side JavaScript code. We've gone through a couple of iterations of how we organise all this. We started with a quick and nasty manual approach to chaining state between different pages of the app. That was pretty horrible. Then, as described in an earlier article, we switched to using AngularJS in the browser, which made the JavaScript side of things much nicer, along with Michael Snoyman's yesod-angular module for managing the interaction between the client and server. Eventually, we found that Angular's default routing provider wasn't flexible enough for our needs in the browser (we have quite a few "wizard"-type stateful interactions, plus a fairly complex dashboard with lots of more independent states).

So we decided to switch over to using the Angular ui-router state-based routing system. This is pretty good and is definitely flexible enough for our current needs. It allows you to define a hierarchical structure of UI states, cleanly controlling transitions between states and assigning URLs to states to allow good interaction with the browsers back and forward buttons and to allow bookmarking of application states. The only problem? Not something that the yesod-angular code supports...

I've started using c2hs recently, for a Haskell NetCDF library I'm writing (which might see the light of day in a couple of months). I like c2hs, so I volunteered to help out with maintenance. The first thing we did was transfer the source code repository from Darcs to GitHub, which was easy, using Steve Purcell's darcs-to-git utility. Once that was done, the next task was to transfer all the tickets on the c2hs Trac site to GitHub's issue tracking. That's what I want to talk about here.

For the past few months, I've been spending most of my working hours on BayesHive, the brainchild of Tom Nielsen, CEO of OpenBrain. This is a web platform for doing Bayesian data analysis, intended to make these powerful methods more accessible for scientists and other data analysts.

We opened the system to some alpha users a while ago, but this week we've decided we're ready for wider beta testing, for people to have a look, kick the tyres and tell us what they think.

I've mostly been working on the web app part of the system, but in this post I want to talk a little bit about the general idea behind the system and the Bayesian approach to statistics. (As far as BayesHive goes, the system is made up of a reasonably complicated single page web app using the AngularJS framework, a Yesod server process to deal with the management of documents, data and models, and a model compiler and inference engine that does all the clever stuff.)

I've been working on BayesHive with Tom Nielsen for something like six months now, and it's been a lot of fun. We've been ironing bugs out of the back-end of the system, which implements a probabilistic functional language (called Baysig). We found an entertaining thing the other day -- I raised the bug report, Tom found out what was going on, and I laughed quite a lot when he told me.

We use the Scrap Your Boilerplate generics framework for doing various transformations of the AST of the Baysig language. It's a great way to define transformations of complex hierarchical data structures in a nice compact way. It does have a trap for the unwary though...

TL; DR: use package version number bounds in the build-depends clauses of your Cabal files. You may save yourself some time and frustration.

Most Haskell programmers seem to have Cabal war stories to tell. Library version management in Haskell is complicated by the aggressive inter-module optimisation that GHC tries to do, which makes Cabal's life tough and sometimes necessitates blowing away both Cabal and GHC's knowledge of installed packages to start over. It's also frustratingly easy to inadvertently get into a state where you have versions of packages installed that prevent Cabal from finding a consistent set of packages to install to satisfy the constraints for a new package you're trying to install.

I thought I'd cleverly avoided most of these problems by using the hsenv sandboxing tool, which allows you to isolate the set of packages needed for a particular project, and even to maintain projects that use different versions of GHC, and cabal-meta, which allows you to build multiple packages from source together so that Cabal can find a set of package versions that are acceptable for all of the packages at once. I often found myself succumbing to a smug sense of self-satisfaction as I read about other peoples' Cabal problems on haskell-cafe.

Of course, pride goeth before a fall. My hubristic little bubble was burst the other day when I decided to reset my sandbox for a project I'm currently working on. For some reason, typing cabal-meta install just wasn't working, and I was being told that there was no set of package versions that could consistently satisfy the requirements of the packages I was trying to build from source.

In the last three days, I've averaged 10 hours of billable time per day. That's ten hours per day actually standing at my desk coding, perhaps 10-20% more than that "on the job". I've been working on social and search aspects of the BayesHive web app: we have a pile of different…

So, I lied a little in my last post. I didn't quite get to go off right away and update my Beeminder totals, because I had a bit a problem with my blogging software. I use Jasper Van der Jeugt's Hakyll, a Haskell static site generator. Jasper recently (well, not that recently...) released a new major version of Hakyll, Hakyll 4. I had been putting off upgrading because I have so much custom code that it's kind of a big job -- my blog is a little different to most Hakyll sites, so I had to write quite a bit of extra stuff to work around that.

So, I've not done any blogging since the ill-fated Midwinter Polyphasic Sleep Experiment. I have a backlog of articles I want to write, but most of those involve some significant investment of time and effort. So you get this instead...

Context

I've been doing a lot of JavaScript programming recently. One of the two contracts I'm working on at the moment is for a UK university spin-off called OpenBrain, who are trying to bring Bayesian statistics to the masses by making tools that allow scientists to build complex (or simple!) Bayesian models and have them analysed using Markov chain Monte Carlo methods. All the really clever stuff is in the server-side MCMC inference code, but the user-facing side of things is a web app with a Yesod backend. My mission is to prettify this and make it ready for prime-time.

It's a fairly complicated application that needs to manage literate Markdown documents describing statistical models, data sets for analysis and a bunch of other stuff. The first thing I looked at was a tool for building statistical models based on dynamical systems represented as systems of ordinary or stochastic differential equations. This tool allows you to give algebraic representations of ODEs or SDEs, which are then rendered as MathML (for prettiness and familiarity), you can do simulations of the systems and visualise the results, you can specify data sets that correspond to observations of system variables, you can set up prior distributions for system parameters, and you can then get some code (in a proprietary Haskell-like language specialised for statistical modelling) that can be used to drive MCMC inference for the parameter values.

All that stuff happens on the client side, except for the listing of the available data sets and the rendering of the model into code. This means that you need client side code for expression parsing, analysis of coupled systems of ODEs and SDEs (to turn them into a canonical form that you can simulate from), numerical integration of ODEs and SDEs, graphing, and a framework to tie it all together. That's a lot of client side code!

This won't sound like fun (or Fay!) to start with, but we'll get there...

In the previous article, we looked at some of the details of a scheme for calculating the Lyapunov exponents of a system of differential equations. In this post, we'll see how to implement this scheme in Haskell. The code we describe is available as a Gist here.

The archetypal chaotic dynamical system is the Lorentz system, defined by the differential equation system

where , and are real constants. On the face of it, this looks like a fairly innocuous set of differential equations. However, as is well known, for a wide range of choices for the parameters , and , this system of equations has chaotic orbits, with sensitive dependence on initial conditions and an attractor of fractal dimension.

Determining whether a given dynamical system will exhibit chaotic behaviour is, in general, difficult . The orbit structure of chaotic systems can be very complex, and it can be difficult to distinguish true chaos from a number of related phenomena.

One approach to characterising the orbit structure of dynamical systems is via Lyapunov exponents, a spectrum of values that measure the stretching and squashing of orbits -- "normal" chaotic systems normally have at least one positive Lyapunov exponent and are dissipative, which means that the sum of all of their Lyapunov exponents is negative. While it is relatively straightforward to define Lyapunov exponents, in practice calculating them is tricky. We'll see how to do this in some practical cases in this and the next couple of articles, following mostly the approach of Rangarajan et al. (1998).

The motivation for this series of articles was basically that I have recently been wondering about how practical it would be to use Haskell for more "traditional" mathematics: there is a big body of work and code relating to category theory and type theory in Haskell, but I've not seen so much about matrix algebra, differential equations, dynamical systems and so on. I implemented the Rangarajan method in Mathematica some years ago and thought it might be a good test case. I'm going to gloss over a lot of technical details about the computation of Lyapunov exponents since my main interest is in the Haskell implementation issues.

So, GHC 7.6.1 was released last month. And there was a new release of Cabal too, bumping the version number to 1.16 from 0.14. The Arch Haskell people did a great job of quickly producing Arch packages for the new version, which meant that everyone pulled the new versions when they did an update via pacman -Su. Unfortunately, there seems to be a problem. The new cabal-install version can't read its own default configuration file and just fails more or less immediately. There seems to be a fix, but it's not yet percolated into the Arch packages.

Michael Snoyman recently had a blog post about composability in the Yesod web framework for Haskell, where he responded to comments about it being difficult to build reusable components for Yesod by, well, building one. I'm pretty much a newbie with Yesod, but the composability aspect of things had never looked too difficult to me, so I thought I'd also have a go at the same kind of exercise.

I decided to implement a slider form field, using one of the nicest jQuery-based slider widgets I've found. All the code for this example is available as a Gist. I'll follow the usual convention here of putting more or less all of the code in the blog post too.

I was originally planning to write a quick-and-dirty implementation of the simplex algorithm myself to demonstrate some of the gritty details, but I decided to leave that kind of thing for later, since there are going to be other, less well-known, constraint solving algorithms that I want to look at . There is a Haskell implementation of a simple version of simplex algorithm in the Matrix.Simplex module in the dsp package, although that comes with a comment that says "I only guarantee that this module wastes inodes", so let's view this more as a way to play with some interface issues than a serious attempt to implement any kind of constraint solver! The code that goes with this article is available here.

I've recently started making some tiny contributions to Brent Yorgey's Haskell diagrams library, mostly bug fixes. The bug fixes are primarily just a way of familiarising myself with the diagrams codebase though. I have my eye on a rather chunkier task in the open issues list. That's writing a constraints solver for diagram layout. I've been interested in this sort of problem for some time, and have never had a good reason to get down to it, so I'm going to see what I can come up with.

And I'm going to try something I've not done before, which is to document experiments, design and development as blog articles. It may all turn out to be hideously embarrassing, it may end up that none of what I do actually makes it into the mainstream of the diagrams library, but having to write about what I'm doing will keep me honest, force me to be clear about what I'm doing, and will provide an incentive to work on this stuff.

The second week of the online Stanford AI class was about probability, reasoning when there is uncertainty and, more specifically Bayes networks. I don't much want to talk about the theory behind these things here since that's covered in a lot of detail on the linked Wikipedia page, in the lectures and in AIMA.

Instead, I want to present a mildly amusing, although not very efficient, Haskell implementation of Bayes networks. It's often said that if you know the full joint probability distribution function (PDF) for a system, you can calculate any marginal or conditional probabilities you want. Well, Haskell is a functional programming language, a PDF is a function, so can we represent the PDF explicitly and use it for calculating? The answer, of course, is yes.

As blogging software, I use the really rather nice Hakyll, by Jasper Van der Jeugt. This is a static website generator written in Haskell, and it's a great solution for smaller blogs or personal websites.

Hakyll works just fine out of the box, but the blog setup I wanted was a bit different from what I've seen other people do with it, so a bit of hacking was required. I wanted to share some of the things I've done, since they might be of interest to other Hakyllers (Hakyllites?). All the code is available on Github.

The first week of the online Stanford AI class has gone by. The format of the presentation is pretty good, with lots of short lecture videos chained together with little quizzes embedded. It's really pretty neat. There are some minor errors in some of the lectures, which is normal, though in this kind of online setting, it's more difficult for the keener students to correct the lecturer as the lecture goes along, so some people are likely to spend a bit of time confused because of that. Some people have complained on Reddit about the relatively superficial level of the coverage of the lectures, but as far as I'm concerned, that's what the book (and papers) are for. The course is called "Introduction to Artificial Intelligence", after all.

The two things we covered this week were a very quick introduction to the field, then some stuff about search. Although there is a big pile of code available for download associated with Russell & Norvig's AI book which would have made playing with some of this stuff easier, I decided to write some code of my own in Haskell, mostly for my own amusement. I wanted to implement A* search, with a clean interface for setting up problems.

Last time, we looked at the basics for a webcomic scraper. In this post, we'll look at extending our code a little to incorporate some date handling and some simple XML processing.

I'm a relative beginner with Haskell, and like many people, to start with I was a little perplexed by the Haskell approach to I/O. A small worked example helped a lot. I was curious to see how easy it would be to do something like the webcomic scraper application implemented in Clojure here and here. This is a simple application, but it does do realistic I/O, downloading files from the web, writing them to disk, and also doing some computations on the file contents. Over the course of two articles, I'm going to build something comparable in Haskell. It turns out to be pretty easy!