Running out of memory with recusion in Haskell

Question

Yesterday I finally decided to start learning Haskell. I started scanning through tutorials, but quickly decided practical exercises would be more beneficial. Therefore I proceeded to port a python script of mine which supposedly simulates gravity into Haskell. To my surprise it actually worked and the generated values match those of python.

I realize the implementation is probably absolutely terrible. The horrible lack of performance does not bother me so much, but what bothers me is that I keep running out of memory when attempting to run the simulation for a longer period of time. Is this because the implementation is inherently flawed or can it be made to work?

I have tried to construct the main loop with three different approaches: "iterate", a recursive function (I read about tail recursion, but didn't manage to make it work) and with a more experimental recursive do function. The functions in question are named simulation, test and test2. I compiled the program with option "-O2".

Why does the program run out of memory, and what can I do to prevent that?

Not so relevant parts of the code:

import System.Environment
import Data.List (tails)
import System.CPUTime
import Text.Printf
import Control.Exception

gConst = 6.674e-11

data Vector = Vector2D Double Double | Vector3D Double Double Double deriving (Show)
deltaVector :: Vector -> Vector -> Vector
deltaVector (Vector2D x1 y1) (Vector2D x2 y2) = Vector2D (x2 - x1) (y2 - y1)
deltaVector (Vector3D x1 y1 z1) (Vector3D x2 y2 z2) = Vector3D (x2 - x1) (y2 - y1) (z2 - z1)

data Position = Position Vector deriving (Show)
data Velocity = Velocity Vector deriving (Show)

distance2DSquared (Vector2D deltaX deltaY) = deltaX ** 2 + deltaY ** 2
distance3DSquared (Vector3D deltaX deltaY deltaZ) = (distance2DSquared $ Vector2D deltaX deltaY) + deltaZ ** 2
distance vector = sqrt (distance3DSquared $ vector)

force vector mass1 mass2 = gConst * (mass1 * mass2) / (distance3DSquared vector)
acceleration force mass = force / mass

vectorComponentDivide (Vector2D x y) c = Vector2D (x/c) (y/c)
vectorComponentDivide (Vector3D x y z) c = Vector3D (x/c) (y/c) (z/c)

vectorComponentMultiply (Vector2D x y) c = Vector2D (x*c) (y*c)
vectorComponentMultiply (Vector3D x y z) c = Vector3D (x*c) (y*c) (z*c)

vectorComponentAdd (Vector2D x1 y1) (Vector2D x2 y2) = Vector2D (x1+x2) (y1+y2)
vectorComponentAdd (Vector3D x1 y1 z1) (Vector3D x2 y2 z2) = Vector3D (x1+x2) (y1+y2) (z1+z2)

invertedVector (Vector2D x1 y1) = Vector2D (-x1) (-y1)
invertedVector (Vector3D x1 y1 z1) = Vector3D (-x1) (-y1) (-z1)

normalizedVector :: Vector -> Vector
normalizedVector vector = vectorComponentDivide vector $ distance vector

velocity vel0 mass1 mass2 vector deltaT = 
    vectorComponentMultiply (vectorComponentAdd vel0 (vectorComponentMultiply (normalizedVector vector) (acceleration (force vector mass1 mass2) mass1))) deltaT

data Actor = Actor String Vector Vector Double deriving (Show)

earth = Actor "Object1" (Vector3D 0 0 0) (Vector3D 0 0 0) 10
moon = Actor "Object2" (Vector3D 10 0 0) (Vector3D 0 0 0) 10

actors = [earth, moon]

combinations :: Int -> [a] -> [[a]]
combinations 0 _  = [ [] ]
combinations n xs = [ y:ys | y:xs' <- tails xs
                           , ys <- combinations (n-1) xs']

updateVelocity [(Actor name1 position1 velocity1 mass1),(Actor name2 position2 velocity2 mass2)] = 
    [(Actor name1 position1 a mass1),(Actor name2 position2 b mass2)]
    where a = velocity velocity1 mass1 mass2 vector deltaT
          b = velocity velocity2 mass2 mass1 (invertedVector vector) deltaT
          vector = deltaVector position1 position2
          deltaT = 1

updatePosition [(Actor name1 position1 velocity1 mass1),(Actor name2 position2 velocity2 mass2)] =
        [Actor name1 (vectorComponentAdd position1 velocity1) velocity1 mass1, Actor name2 (vectorComponentAdd position2 velocity2) velocity2 mass2]

Relevant parts:

update list = map updatePosition (map updateVelocity list)


simulation state n = do
    if n == 0
        then do
            print state
            return ()
        else do
            let newState = update state
            simulation newState $! (n-1)

test list n = iterate update list !! n

test2 list 0 = list
test2 list n = (test2 (update list) (n-1))

time :: IO t -> IO t
time a = do
    start <- getCPUTime
    v <- a
    end   <- getCPUTime
    let diff = (fromIntegral (end - start)) / (10^12)
    printf "Computation time: %0.3f sec\n" (diff :: Double)
    return v

main :: IO ()
main = do
    let combo = combinations 2 actors
    putStrLn "Hello World!"
    let n = 1000000
    print n
    --time $ print (test combo n)
    time $ simulation combo n
    _ <- getLine
    putStrLn "BAI"

Answer 1

I believe laziness harms your code: your code builds large thunks (unevaluated expressions) which lead to OOM.

For instance, iterate is (in)famous for leading to large thunks, when you access the resulting list in the middle without forcing the previous list elements. More precisely

iterate f x !! n

is bad, since it will build the expression f (f (f ...(f x))) before really performing any work. We want to evaluate every list element before accessing the next one. This could be dome by a custom !! function:

(!!!) :: [a] -> Int -> a
[]     !!! _ = error "!!!: out of range"
(x:_ ) !!! 0 = x
(x:xs) !!! n = seq x (xs !!! pred n)

Now we can use iterate f a !!! n without large thunks building up.

This has the same problem:

simulation state n = do
    if n == 0
        then do
            print state
            return ()
        else do
            let newState = update state
            simulation newState $! (n-1)

It will build large update (update (update ...)) thunks without evaluating them. A possible fix could be

            ...
            (simulation $! newState) $! (n-1)

However, keep in mind that in your case newState is a list (of lists!). In such case, seq or $! will only demand the list to be evaluated as far as its first cell constructor -- just enough to check whether the list is empty or not. This "forcing" might be enough or not for your purposes.

There is a library function named deepSeq which will force the full list, if really needed (use Hoogle to find the docs).

Summing up: lazy evaluation has its benefits and its downsides. It usually allows for more efficiency, e.g. sometimes providing constant space list processing without the need of writing carefully crafted functions. It also allows infinite lists tricks which are handy. However, it can also cause unwanted thunks to stick around for too long, wasting memory. So, in those cases, programmers have some burden put on them. Especially when one is used to strict semantics, these issues can be scary at first (we've been there!).