Sharing Computation of Force in Haskell

Question

I'm implementing a N-Body simulation in Haskell. https://github.com/thorlucas/N-Body-Simulation

Right now, each particle calculates its Force, then acceleration against each other particle. In other words, O(n²) computations of force. I could reduce this down to O(n choose 2) if I were to calculate each combination once.

let combs = [(a, b) | (a:bs) <- tails ps, b <- bs ]
    force = map (\comb -> gravitate (fst comb) (snd comb)) combs

But I cant figure out how to apply these to the particles without using state. In the above example, ps is [Particle] where

data Particle = Particle Mass Pos Vel Acc deriving (Eq, Show)

Theoretically, in an stateful language I would simply be able to loop through the combinations, calculate the relevant acceleration from the force from each a and b, then update each Particle in ps acceleration as I do that.

I've thought about doing something like foldr f ps combs. The starting accumulator would be the current ps and f would be some function that takes each comb and updates the relevant Particle in ps, and returns that accumulator. That seems really memory intensive and quite complicated for such a simple process.

Any ideas?

Answer 1

Taking the code from https://github.com/thorlucas/N-Body-Simulation

updateParticle :: Model -> Particle -> Particle
updateParticle ps p@(Particle m pos vel acc) =
    let accs = map (gravitate p) ps
        acc' = foldr (\(accx, accy) (x, y) -> (accx + x, accy + y)) (0, 0) accs
        vel' = (fst vel + fst acc, snd vel + snd acc)
        pos' = (fst pos + fst vel, snd pos + snd vel)
    in  Particle m pos' vel' acc'

step :: ViewPort -> Float -> Model -> Model
step _ _ ps = map (updateParticle ps) ps

and modifying it so the accelerations are worked out in a matrix (well, list of lists...) separately from updating each particle, we get...

updateParticle :: Model -> (Particle, [Acc]) -> Particle
updateParticle ps (p@(Particle m pos vel acc), accs) =
    let acc' = foldr (\(accx, accy) (x, y) -> (accx + x, accy + y)) (0, 0) accs
        vel' = (fst vel + fst acc, snd vel + snd acc)
        pos' = (fst pos + fst vel, snd pos + snd vel)
    in  Particle m pos' vel' acc'

step :: ViewPort -> Float -> Model -> Model
step _ _ ps = map (updateParticle ps) $ zip ps accsMatrix where
    accsMatrix = [map (gravitate p) ps | p <- ps]

... so the problem is essentially how to reduce the number of calls to gravitate when working out accsMatrix, using the fact that gravitate a b = -1 * gravitate b a.

If we were to print out accsMatrix, it would look like...

[[( 0.0,  0.0), ( 1.0,  2.3), (-1.0,  0.0), ...
[[(-1.0, -2.3), ( 0.0,  0.0), (-1.2,  5.3), ...
[[( 1.0,  0.0), ( 1.2, -5.3), ( 0.0,  0.0), ...
...

... so we see accsMatrix !! i !! j == -1 * accsMatrix !! j !! i.

So to use the above fact we need access to some indexes. Firstly, we index the outer list...

accsMatrix = [map (gravitate p) ps | (i,p) <- zip [0..] ps]

... and replace the inner list with a list comprehension...

accsMatrix = [[ gravitate p p' | p' <- ps] | (i,p) <- zip [0..] ps]

... get some more indexes available via zip...

accsMatrix = [[ gravitate p p' | (j, p') <- zip [0..] ps] | (i,p) <- zip [0..] ps]

... and then, the key, is make accsMatrix depend on itself for half of the matrix...

accsMatrix = [[ if i == j then 0 else if i < j then gravitate p p' else -1 * accsMatrix !! j !! i | (j, p') <- zip [0..] ps] | (i, p) <- zip [0..] ps]

We can split it up a bit as well, as below...

accsMatrix = [[ accs (j, p') (i, p) | (j, p') <- zip [0..] ps] | (i, p) <- zip [0..] ps]
accs (j, p') (i, p) 
    | i == j    = 0
    | i < j     = gravitate p p'
    | otherwise = -1 * accsMatrix !! j !! i

... or avoid list comprehensions by using map

accsMatrix = map (flip map indexedPs) $ map accs indexedPs  
indexedPs = zip [0..] ps
accs (i, p) (j, p')
    | i == j    = 0
    | i < j     = gravitate p p'
    | otherwise = -1 * accsMatrix !! j !! i

... or by using the list monad...

accsMatrix = map accs indexedPs >>= (:[]) . flip map indexedPs

... although it's harder (for me) to see what's going on in these.

There are probably some horrible performance issues with this list of lists approach, especially using !!, and the fact that you're still running O(n²) operations due to the traversals, and the fact that O (n · (n ­– 1)) ≡ O (n²) as @leftaroundabout mentioned, but each iteration should call gravitate n * (n-1) / 2 times.