Last night in #haskell, there was some interesting discussion going on, in which we discovered a class of data structures that seems to fit in between Applicative and Monad. A foreword though: Most of the credit for this article should go to Brent Yorgey, Cale Gibbard, Peaker and Tristan Seligmann – the article is the product of a discussion with them last night in which they did most of the work.

### Why doesn’t Applicative cut the mustard

One of the major drawbacks of Applicative is that you can’t make choices without evaluating too much. If we attempt to use liftA3 if', we find that we evaluate both branches, no matter what we get back from our predicate.

### So just use Monad, right?

Well no, bind/join give us too much power, they let us make choices between an infinite number of functions, we only want the ability to make a choice between a finite number of branches.

We can define this class then:

class Applicative b => Branchy b where
ifB :: b Bool -> b a -> b a -> b a

We can add a pair of laws to this that force this to really be an if like operation, and to make sure that we get lazyness properly:

ifB (pure True ) x y = x
ifB (pure False) x y = y

These laws though may not be complete in getting us to a useful class. Cale has suggested this additional law to get us closer, but holes in the laws are a definite place to do some work hashing this out.

f <*> ifB b x y = ifB b (f <*> x) (f <*> y)

So then, do we actually have something useful here? For that, we would have to find two things – an Applicative which is not a Branchy (to show that we’ve not just reinvented applicatives), and a Branchy that’s not a Monad. Well, clearly Const a is not a Branchy, a Const a Bool does not actually contain a Bool with which to decide which branch to take. And, here’s the useful Branchy that’s not a monad:

instance Branchy ZipList where
ifB (ZipList ps) (ZipList xs) (ZipList ys)
= ZipList \$ zipWith3 if' ps xs ys

Andrea Vezzosi commented this morning that the interface would be rather nicer if rather than using ifB, we used eitherB :: b (Either c d) -> b (c -> e) -> b (d -> e) -> b e. While each can be implemented in terms of the other, eitherB allows for a little more efficiency. The following implementation of eitherB evaluates the either twice:

eitherB fe fl fr = ifB (isLeft       <\$> fe)
((.fromLeft)  <\$> fl)
((.fromRight) <\$> fr) <*> fe

This version of the class also makes the connection with ArrowChoice rather clearer.

I’d really appreciate everyone’s comments on this topic.

# Simulating n-bodies and functional programming

### The n-bodies problem

There’s been some activity in the #haskell IRC channel recently towards coming up with a better solution to the n-bodies problem.  In this problem, you’re given a set of bodies, and are asked to simulate Newtonian gravity on them.  Most of the solutions I’ve seen go past have involved some sort of crazy state monad based stuff, reading from and writing to STUArrays and other imperative ideas.  So, I thought I’d try and solve the problem using a pure functional style instead.

### Reactive

I’ve recently been playing with Reactive a lot, initially for my work, but also just because it’s so damned cool.  Reactive is a library for Haskell that lets you describe time varying values in a purely functional way.  That is, out goes the concept of state, and in comes the concept of describing exactly what is going on.

### Back to n-bodies

In order to simulate newtonian physics, we need the bodies’ mass,the positions of the bodies, and the velocities of the bodies:

data Body = Body { mass     :: Double
, position :: Point3  Double
, velocity :: Vector3 Double
}

That was pretty straight forward, lets get down to the core of the problem – simulating gravity:

accelBetween :: Body -> Body -> Vector3 Double
accelBetween p1 p2 =
(mass p2 *^ vecBetween) / ((mag vecBetween) ^ 3)
where
vecBetween = position p2 ^-^ position p1

computeAccels :: Body -> [Body] -> Vector3 Double
computeAccels p = sum . map (accelBetween p)

Well, that was surprisingly easy!  The accelBetween function describes the newtonian gravity equation just as the maths does, and summing up the acceleration due to all the other planets in the system is fairly straightforward.

What we’ve not seen yet though, is any of that lovely functional reactive programming I was talking about.  We know now how to compute the acceleration affecting any one body at any one time, but what we want to know is the acceleration affecting a body at *any* time.  To do this, we’re going to need to move from doing computations on Bodys, to instead doing computation on Behaviors of Bodys.

bodyAccels :: [Behaviour Body] -> [Behaviour (Vector3 Double)]
bodyAccels ps = withOthers (liftA2 computeAccels) ps

withOthers :: (a -> [a] -> b) -> [a] -> [b]
withOthers f xs =
withOthers' f [] xs
where
withOthers' _ _  []       = []
withOthers' f os (x : xs) = f x (os ++ xs) : withOthers' f (x : os) xs

The withOthers function here is just like map, but it passes in all the other values in the list in a second argument to the function.

So then, bodyAccels computes a continuous acceleration function for all the bodies in the system.  For each body, it runs the computeAccels function, giving it all other bodies in the system as its second argument.  Crucially, liftA2 allows us to do this in the Behavior Applicative, so we are no longer computing it on rigid, static bodies, but instead, on all the positions and velocities the bodies in the system will ever have (isn’t lazyness great!).

Finally, we can get from these accelerations down to the velocities, and then positions of the bodies using integration on the acceleration:

bodyVel :: Body -> Behaviour (Vector3 Double) -> Behaviour (Vector3 Double)
bodyVel p acc = (velocity p ^+^) <\$> integral dt acc

I’m sure you can imagine what the bodyPos function looks like.  You may wonder what the dt in here is talking about.  This is an unfortunate effect of not being able to mathematically integrate arbitrary functions.  Instead, we must use euler integration, and that requires us to provide times at which to take samples.  The dt argument is an event which ticks reasonably fast, and progresses our simulation:

dt :: Event ()
dt = atTimes [0,0.01..]

So, now we are able to combine all our efforts together, and solve the whole n-bodies problem:

nbodies :: [Body] -> [Behaviour Body]
nbodies ps = pbs
where
pbs = Body <\$> (mass <\$> ps) <*> pps <*> pvs
pps = zipWith bodyPos ps pvs
pvs = zipWith bodyVel ps pas
pas = bodyAccels pbs

### Conclusions

A lot of assumptions are made about how we must write programs.  Often, even beautiful mathematical problems end up described as horrible state-full masses of code that obscure what it is we’re trying to compute.  I’ve presented a solution to the  n-bodies problem using the Reactive library to get a handle on time in a purely functional setting, it turned out that this was rather beautiful!

Full code for my solution can be found below.