Saturday, 12 June 2010

Orbit Simulator in Haskell

Uncle Bob recently wrote an article showing a simple orbital simulator in Clojure. Aside from the strange formatting, the code is pretty easy to follow, though there is large amounts of duplicate code (see vector.clj vs. position.clj). I guess the reasoning behind this is to try to stop mixing up units?

A phantom type is one way to avoid this duplication in Haskell. A phantom type is only ever used to construct types. Its sole purpose is for validation by the type checker. The example below uses phantom types to encode position, velocity and force. This means we have no code repetition for each of the vector operations.



When we do the physical calculations to move the objects around the phantom types do their job. As an example, if I try to add two unrelated vectors together I get a compile time error showing that I can't. Here's an example from a ghci session.

  Prelude Orbit> let a = Vec 0 0 :: Vec Force
Prelude Orbit> let b = Vec 1 1 :: Vec Position
Prelude Orbit> add a b

:1:6:
Couldn't match expected type `Force'
against inferred type `Position'
Expected type: Vec Force
Inferred type: Vec Position
In the second argument of `add', namely `b'
In the expression: add a b


The code to do the physical calculations is very simple. The record update syntax is new to me and means you can just change the selected fields in an object, rather than creating a new one and copying across all the fields (see the section on named fields).

I think collideAll looks substantially nicer than the Clojure implementation. From the comments on the blog, I suspect it still suffers from a period of inaccuracy if three objects collide simultaneously (they won't get merged in the first round, but almost certainly will in the second). nubBy is an incredibly useful function I hadn't seen before and it allows you to remove duplicates according to your own definition of equality.



The Clojure version has a bunch of unit tests (as you might expect given the author!). The tests are in the 1 + 1 = 2 style by which I mean give some inputs and validate an output. These are the sort of tests you might capture as part of a REPL transcript, but they give you very little assurance that the code is actually correct.

QuickCheck gives an incredibly powerful way of testing applications by making assertions about the code and challenging QuickCheck to falsify them with generated data. In the case of this simple physics model, we can make assertions like the total amount of energy in the system must be conserved and let QuickCheck worry about generating the cases where this doesn't hold.

In order to do this we first have to give QuickCheck a way of generator random objects to simulate. The Arbitrary type class defines generate, so we just need to make Object an instance of this type class and implement the method.

A few simple quick check properties are shown below. This asserts that for any given vector (apart from one with zero length), the magnitude is 1.0 (subject to a rounding error or two) and that for any given set of objects, the kinetic energy is conserved. This is a very simple example of how powerful QuickCheck can be.



Running Quick Check against this property shows that it holds for a large number of automatically generated test cases. This gives a strong hint that this code actually works.

The final task is to visualize this data. I used OpenGL again and knocked up a quick UI, the source code for which is here. The (terrible) video below gives an example of it in action. If anyone knows of a better way of capturing OpenGL I'd love to know. I'm currently using Istanbul which works, but brings my system to its knees.

video

The complete source for this is available on my git hub page. Any comments or suggested improvements greatly appreciated!

Finally, the company I work for is currently recruiting in Cambridge, UK. If you're the sort of person who is looking for challenges such as analysing petabytes of data, writing high performance servers , or developing the richest pure JavaScript interfaces imaginable (and a whole lot more in between) then drop me a note at jeff.foster AT acm.org and I'll give you more information. Oh, and if you're accepted you'll get your choice of an iPad or the latest Android phone.