[Haskell-cafe] Squashing space leaks

[Josef Svenningsson]

On 5/6/05, Daniel Fischer <daniel.is.fischer at web.de> wrote:
> Am Freitag, 6. Mai 2005 02:24 schrieb Greg Buchholz:
> >     We're heading in the right direction anyway.  I can now compute 1
> > million iteration in about 2 minutes (with +RTS -H750M -K100M).  Well
> > almost, since it now doesn't compute the right answer, so something must
> > be amiss in the shuffling section.  Now if we can get it to us a little
> > less than 1G of memory, we'll be in good shape.
> >
> > Thanks,
> >
> > Greg Buchholz
> 
> The problem is that a prime crept in in Josef's code,
> so to calculate the positions and velocities, the updated versions of the
> planets are used, it should be
> 
> update f newlist (a:as) = a' : update f (newlist . (a:)) as
>                where a' = f a (newlist as)
> 
> instead of
> update f newlist (a:as) = a' : update f (newlist . (a':)) as
>                where a' = f a (newlist as).
> 
The prime didn't creep in, I put it there on purpose. Which just shows
that I didn't understand the algorithm properly :). But your
correction makes the algorithm a little easier to express. But first
to the space problem.

Greg, you express the sequential computation of updating the set of
planets as an iteration over a list. This is a sweet way of expressing
it but doesn't work very well when an element of the list contain
pointers to the previous element. This will force all the elements to
be in memory at the same time. What causes this is too little
strictness, all computations which depend on the previous element
haven't been forced and thus still points to it. It this specific
program the problem is the list of planets. Even if we may have
computer ourselves a new list of planets the planets themselves
haven't been computed yet and will contain pointers to the previous
list of planets. Devastating! Hence we will need a strict list to do
this or some other strict data structure.

So, here's my code. It's not beautiful but it solves the space
problem. I hope I got it right this time :)
\begin{code}
data StrictList a = Cons !a !(StrictList a) | Nil

slToList (Cons a sl) = a : slToList sl
slToList Nil = []

listToSl (a:as) = Cons a (listToSl as)
listToSl [] = Nil

mapSL f (Cons a sl) = Cons (f a) (mapSL f sl)
mapSL f Nil = Nil

partitionSL :: StrictList a -> StrictList (a,StrictList a)
partitionSL sl = go id sl
  where go prev Nil = Nil
	go prev (Cons a sl) = Cons (a,prev sl) (go (prev . (Cons a)) sl)

advance dt sl = mapSL newplanet (partitionSL sl)
  where newplanet (p,sl) = Planet (pos p + delta_x) new_v (mass p)
	  where ps = slToList sl
		delta_v     = sum (map (\q -> 
		 (pos p - pos q) `scale` ((mass q)*dt/(dist p q)^3)) ps)
		new_v       = (vel p) - delta_v
		delta_x     = new_v `scale` dt

offset_momentum ((Planet p v m):ps) = Cons (Planet p new_v m) (listToSl ps)
  where new_v = (sum (map (\p->(vel p) `scale` (mass p)) ps)) 
                `scale` ((-1.0)/solar_mass)

energy:: Int -> StrictList Planet -> Double
energy n ls = kinetic - potential
  where
    ps = slToList ls
    kinetic   = 0.5 * (sum $ map (\q-> (mass q)*((vel q)`dot`(vel q))) ps)
    potential = sum [(mass (ps!!i))*(mass (ps!!j))/(dist (ps!!i) (ps!!j))
                      | i<-[0..n-1], j<-[i+1..n-1]]
\end{code}

Cheers,

/Josef