»
S
I
D
E
B
A
R
«
Co-routines in Haskell
Jul 21st, 2010 by blackh

It is easy to implement co-routines in Haskell… but only if you know how. No fewer than three people asked me to blog about it, so here’s a quick guide to rolling your own co-routines. To understand this blog, you will need to have a basic understanding of monad transformers.

There are co-routine packages on Hackage, but I have not had much luck with them. The point here, really, is to show you how it all works.

What’s a co-routine?

A co-routine (called a ‘generator’ in Python) is where you create two interleaved flows of control on a single thread. Unlike threads, co-routines switch co-operatively using a ‘yield’ operation. (This is quite a good trick in GUI programming for implementing complex workflow that spans multiple GUI events, since most GUI libraries require everything to be on one thread.)

The example I’m presenting here works in this way: A caller executes a CoroutineT monad transformer, which adds the ‘yield’ operation to the underlying monad (which can be anything). From the caller’s point of view, the ‘yield’ looks like the monad has returned, but with a continuation. In the callee, ‘yield’ appears to block until the caller executes the continuation. In addition, we add the ability to pass a value in both directions.

So we’ve inverted the flow of control in the callee. Continuation passing style (CPS) can also do this, but co-routines are better than CPS because 1. it’s a bit neater, and 2. it allows for recursion.

One application of co-routines is to separate I/O from logic. By way of example I am going to implement an expert system for identifying fruit. I try not to use contrived examples, and as you can see, this time I have completely failed.

In this example, the CoroutineT sits on top of the identity monad, so it’s pure, but the approach works just the same on top of IO or anything else. This example is not deeply nested, but this approach happily supports any level of recursion or nesting.

So here’s our expert system logic. We’ll define CoroutineT shortly. You can read ‘yield’ as ‘askUser’:

type Question = String
data Answer = Y | N deriving Eq
type Expert a = CoroutineT Answer Question Identity a
data Fruit
    = Apple
    | Kiwifruit
    | Banana
    | Orange
    | Lemon
    deriving Show
identifyFruit :: Expert Fruit
identifyFruit = do
    yellow <- yield "Is it yellow?"
    if yellow == Y then do
        long <- yield "Is it long?"
        if long == Y then
            return Banana
          else
            return Lemon
      else do
        orange <- yield "Is it orange?"
        if orange == Y then
           return Orange
         else do
           fuzzy <- yield "Is it fuzzy?"
           if fuzzy == Y then
               return Kiwifruit
             else
               return Apple

Our ‘Expert’ type above…

type Expert a = CoroutineT Answer Question Identity a

…specifies the type we are sending into our co-routine (Answer) and the type we are getting out of it (Question) as viewed from the caller.

Now we just need a main program to drive it. Because the I/O is separated out, we can later replace this with a nice touch-screen GUI for the seriously fruit-impaired.

main :: IO ()
main = do
    putStrLn $ "Expert system for identifying fruit"
    run identifyFruit
  where
    run :: Expert Fruit -> IO ()
    run exp = handle $ runIdentity $ runCoroutineT exp
    handle (Yield q cont) = do
        putStrLn q
        l <- getLine
        case map toLower l of
            "y"   -> run $ cont Y
            "yes" -> run $ cont Y
            "n"   -> run $ cont N
            "no"  -> run $ cont N
            _   -> putStrLn "Please answer 'yes' or 'no'" >> handle (Yield q cont)
    handle (Result fruit) = do
        putStrLn $ "The fruit you have is: "++show fruit

When we run our co-routine, it returns with one of these two events:

  • Yield, which happens when yield is executed. It gives us the output value (Question) and the continuation, which when passed the input value (Answer) gives us the same ‘Expert’ type we started with.
  • Result, which happens when the co-routine has finished executing.

So how does CoroutineT work?

We’ll start with the types:

data Result i o m a = Yield o (i -> CoroutineT i o m a) | Result a
-- | Co-routine monad transformer
--
--   * i = input value returned by yield
--
--   * o = output value, passed to yield
--
--   * m = next monad in stack
--
--   * a = monad return value
data CoroutineT i o m a = CoroutineT {
        runCoroutineT :: m (Result i o m a)
    }

Hopefully that’s pretty straightforward. ‘yield’ is defined like this:

-- | Suspend processing, returning a @o@ value and a continuation to the caller
yield :: Monad m => o -> CoroutineT i o m i
yield o = CoroutineT $ return $ Yield o (\i -> CoroutineT $ return $ Result i)

The key point here is that the continuation does nothing except return the value, which is what we want it to do when we run a monad that contains only a yield.

Most of the magic is in the definition of >>=, thus:

instance Monad m => Monad (CoroutineT i o m) where
    return a = CoroutineT $ return $ Result a
    f >>= g = CoroutineT $ do
        res1 <- runCoroutineT f
        case res1 of
            Yield o c -> return $ Yield o (\i -> c i >>= g)
            Result a  -> runCoroutineT (g a)
    -- Pass fail to next monad in the stack
    fail err = CoroutineT $ fail err

A typical monad would normally execute f then pass its result to g and execute that, and this is in fact exactly what we do in the Result case. Ho hum.

But there’s no law that says you have to execute g. This is Haskell so we can do whatever we like. g is just a plain old closure representing the continuation.

So what we do in the Yield case is take the continuation that executing f gave us, and bind that to the continuation g, then bail out of the monad (in the same way ErrorT does when it gets an error), handing our constructed continuation to the caller. So we end up with a closure that represents the entire execution state of the monad, and it doesn’t matter how deeply nested we are. It just puts the continuation together in the right way as we unravel everything on our way back to the caller.

Here’s the code in downloadable form:

This code is released in the public domain.

Stephen Blackheath, Manawatu, New Zealand

AusHac2010 Day 2 progress
Jul 17th, 2010 by axman6
Day 2 of AusHac2010 is coming to an end, and we’ve made a lot of progress:
Bernie Pope has been making great progress with a new MPI binding for Haskell
Ben Lippmeier, Erik de Castro Lopo and Ben Sinclair have been busily hacking on DDC, with 13 commits today alone
Stephen Blackheath has been working on some code using the Accelerate library that rasterises triangles for use in a commercial computer game.
Hamish Mackenzie, Jens Petersen and Matthew Sellers have been working on better Yi integration for Leksah, working on using Yi’s current configuration file, and improving “launch experience”, focusing on eliminating the requirement of creating an initial workspace file.
Lang Hames has been using his experience with LLVM from working at Apple as an intern to improve various low level problems in LLVM. His work should help resolve some of the problems the LLVM backend to GHC has, but should also be very beneficial to many other LLVM users. While doing this, he’s written a very nice tool that illustrates register liveness, with further work focusing on colouring the HTML output to show register pressure. The LLVM guys seem quite excited about this work, which is great.
Mark Wotton and Sohum Banerjea have been trying to extend Hubris, the Haskell-Ruby bridge, to work with polymorphic functions. Their heads are quite sore from all the head banging. Raphael Speyer has been working on an install script to make installation much easier for users.. but only if you use Ubuntu so far.
Ivan Miljenovic has been prematurely optimising his containers library, before finalising the API. This library is designed to let library writers leave the choice of which container data structure to output to the library consumer as well as making it easier to change which data structure you want to use in your code, with minimal code change. See his blog post for more details. **********
Trevor McDonell has been working on the CUDA backend to Accelerate, adding support for efficient nested tuple types, and other bug fixes. Sean Lee has been helping out with testing of this code, along with Manuel Chakravarty.
With one more full day to go, I think we’ll be getting a lot of awesome work done tomorrowQ!

Day 2 of AusHac2010 is coming to an end, and we’ve made a lot of progress:

Bernie Pope has been making great progress with a new MPI binding for Haskell

Ben Lippmeier, Erik de Castro Lopo and Ben Sinclair have been busily hacking on DDC, with 13 commits today alone

Stephen Blackheath has been working on some code using the Accelerate library that rasterises triangles for use in a commercial computer game.

Hamish Mackenzie, Jens Petersen and Matthew Sellers have been working on better Yi integration for Leksah, working on using Yi’s current configuration file, and improving “launch experience”, focusing on eliminating the requirement of creating an initial workspace file.

Lang Hames has been using his experience with LLVM from working at Apple as an intern to improve various low level problems in LLVM. His work should help resolve some of the problems the LLVM backend to GHC has, but should also be very beneficial to many other LLVM users. While doing this, he’s written a very nice tool that illustrates register liveness, with further work focusing on colouring the HTML output to show register pressure. The LLVM guys seem quite excited about this work, which is great.

Mark Wotton and Sohum Banerjea have been trying to extend Hubris, the Haskell-Ruby bridge, to work with polymorphic functions. Their heads are quite sore from all the head banging. Raphael Speyer has been working on an install script to make installation much easier for users.. but only if you use Ubuntu so far.

Ivan Miljenovic has been prematurely optimising his containers library, before finalising the API. This library is designed to let library writers leave the choice of which container data structure to output to the library consumer as well as making it easier to change which data structure you want to use in your code, with minimal code change. See his blog post for more details.

Trevor McDonell has been working on the CUDA backend to Accelerate, adding support for efficient nested tuple types, and other bug fixes. Sean Lee has been helping out with testing of this code, along with Manuel Chakravarty.

With one more full day to go, I think we’ll be getting a lot of awesome work done tomorrow!

Testing BlogLiterally, and HTML DSL changes
Nov 3rd, 2009 by Axman6

Robert Greayer has just released his BlogLiterally program on hackage, so I thought I’d try it out. I’ve made some small changes to the HTML DSL I’ve been working on since i last posted, which I’ll outline here:

The HTML type is now

data HTML a = Doctype String String String (HTML a)
            | Content a
            | Text String
            | Node Tag Attrs [HTML a]
            | NoContent Tag Attrs
            | Comment String
            deriving Show

instead of

data HTML a = Doctype String String String [HTML a]
            | Content a
            | Text String
            | Node Tag Attrs [HTML a]
            | NoContent Tag Attrs
            | Comment String
            deriving Show

I’ve also added some more helper functions:

addAttrs :: Attrs -> HTML a -> HTML a
addAttrs attrs (Node tag attrs' rest) = Node tag (attrs' ++ attrs) rest
addAttrs attrs (NoContent tag attrs') = NoContent tag (attrs'++attrs)
addAttrs _ x = x

toTable :: [[String]] -> HTML a toTable xss = table . map (tr . map (\x -> td [text x])) $ xss

toTable' :: Show a => [[a]] -> HTML b toTable' = toTable . map (map show)

And updated my example:

example :: HTML ()
example =
    doctype "html PUBLIC"
            "-//W3C//DTD XHTML 1.0 Transitional//EN"
            "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd" $
        xhtml [
            header [
                title "My random page",
                meta "name" "Alex"
            ],
            body [
                comment "This is a comment",
                para [
                    text "The hr tag defines a horizontal rule:"
                ],
                hr,
                para [
                    text "This is a paragraph"
                ],
                hr,
                para [
                    text "This is a paragraph"
                ],
                addAttrs [("border","1px")] . toTable' $ [[0..n] | n <- [0..3]]
            ]
        ]

Which produces this html:

    <!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
       "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
    <html xmlns="http://www.w3.org/1999/xhtml" xml:lang="en" lang="en">
        <head >
            <title >
                My random page
            </title>
            <meta name="name" content="Alex"/>
        </head>
        <body >
            <!-- This is a comment -->
            <p >
                The hr tag defines a horizontal rule:
            </p>
            <hr />
            <p >
                This is a paragraph
            </p>
            <hr />
            <p >
                This is a paragraph
            </p>
            <table border="1px">
                <tr >
                    <td >
                        0
                    </td>
                </tr>
                <tr >
                    <td >
                        0
                    </td>
                    <td >
                        1
                    </td>
                </tr>
                <tr >
                    <td >
                        0
                    </td>
                    <td >
                        1
                    </td>
                    <td >
                        2
                    </td>
                </tr>
                <tr >
                    <td >
                        0
                    </td>
                    <td >
                        1
                    </td>
                    <td >
                        2
                    </td>
                    <td >
                        3
                    </td>
                </tr>
            </table>
        </body>
    </html>

An HTML DSL in Haskell in one day
Nov 2nd, 2009 by Axman6

So, after a discussion on #haskell, I decided to try writing a simple HTML DSL in Haskell. In about 15 minutes, I had a rather rudimentary implementation, and after a day or so of working, I’ve got a somewhat useful DSL for writing HTML, and pretty printing it.

The main type behind the package is:

data HTML a = Doctype String String String [HTML a]
            | Content a
            | Text String
            | Node Tag Attrs [HTML a]
            | NoContent Tag Attrs
            | Comment String
            deriving Show

where

type Attr = (String, String)
type Attrs = [Attr]
data Tag = Html
         | ARef
         | BlockQuote
         | Body
         ...
         | Title
         | TR
         | UL

I’ve also provided helper functions to make this an actual DSL (I guess?):

example :: HTML ()
example = 
    doctype "html PUBLIC" "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd" [
        xhtml [
            header [
                title "My random page",
                meta "name" "Alex"
            ],
            body [
                comment "This is a comment",
                para [
                    text "The hr tag defines a horizontal rule:"
                ],
                hr,
                para [
                    text "This is a paragraph"
                ],
                hr,
                para [
                    text "This is a paragraph"
                ],
                table [
                    tr [
                        td [
                            text "data"
                        ],
                        td [
                            text "more data"
                        ]
                    ],
                    tr [
                        td [
                            text "more data"
                        ]
                    ]
                ]
            ]
        ]
    ]

produces:

<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml" xml:lang="en" lang="en">
    <head >
        <title >
            My random page
        </title>
        <meta name="name" content="Alex"/>
    </head>
    <body >
        <!-- This is a comment -->
        <p >
            The hr tag defines a horizontal rule:
        </p>
        <hr />
        <p >
            This is a paragraph
        </p>
        <hr />
        <p >
            This is a paragraph
        </p>
        <table >
            <tr >
                <td >
                    data
                </td>
                <td >
                    more data
                </td>
            </tr>
            <tr >
                <td >
                    more data
                </td>
            </tr>
        </table>
    </body>
</html>

I don’t claim that anyone should use this (there’s a good chance I won’t be releasing it, unless someone really wants me to), but I’m still quite proud of it.

Reddit fame for Axman!
Aug 2nd, 2009 by Axman6

So, I posted a link to reddit yesterday, thought it was funny, and thought it might get a few votes… I was wrong. It got a lot of votes, and is now the 11th top link submitted to reddit of all time. If you’re curious what all the fuss is about, check it out here.

So, thanks to all of you that voted for it, you’ve made me feel quite loved (in a you liked something I posted but didn’t have anything to do with kind of way).

– Axman

A new web site and a new hackage package.
Jun 29th, 2009 by Axman6

After years of wanting to, I finally bought axman6.com. I’m using the quite awesome nearlyfreespeech.net for my hosting, because they’re damn cheap, and seem to runs things very nicely.

Now for the interesting stuff. I’ve been working over the last few days on implementing a new Set and Map library for haskell, which I released yesterday on Hackage, called imaginatively TernaryTrees. I got the idea from a PC Plus article [Via: Reddit.com], and couldn’t quite get the whole idea into my head at once, so decided to try implementing it in the best language I knew, haskell.

It turns out this was a good idea; without much work at all, I’ve ended up with a rather efficient implementation, which can be used for either efficient Sets of Strings, or Sets or Maps of any Ord a => [a]. The thing I’m most happy about is how nice I managed to make the Data.Binary implementation. By using a little bit of binary number trickery, I’ve managed to save an awful lot of space (The {T,S}End notes only requite 1 bit of information saying whether they’re there or not and these bits are stored in the same byte that marked whether we has a {T,S}Node or {T,S}End). When using this to write the StringSet out to a binary file, I’ve been able to get a the output to be compressed by over 40% in one case (though this doesn’t seem to be typical sadly :( )

There’s not much more I want to say about this for now, maybe later, but I suggest you go and take a look at the source and see how easy this all was.

– Axman

AVar changes
Feb 22nd, 2009 by axman6

Yesterday, i stuck AVar 0.0.4 on hackage. Again, not thinking things through as well ad i should have, i realised later it should have been the 0.1.0 release, since it was the first major change to the package since i first wrote it.

So what’s new? Well, i split it up a bit, into three modules: Data.AVar, the classic interface, Data.AVar.Unsafe, the same as Data.AVar, but it throws any exceptions encountered by the variable, instead of passing them back to the user, and Data.AVar.Internal, which contains the code for the actual for AVars, and the datatypes used by them.

I think this is a nice way of doing things, and hopefully will make using the system a little easier if you don’t expect to cause exceptions. So as always, if you have any suggestions, please let me know.

– Axman

N-bodies speedup (50%!)
Feb 3rd, 2009 by axman6

So I’ve been playing with this n-bodies thing again, and while I haven’t managed to get all that much closer to the likes of C, C++, or even the current Haskell entry. BUT! I did manage to make my program almost exactly twice as fast as it was.

I used a few simple tricks from the performance part of the haskell wiki, mainly the use of unsafeRead and unsafeWrite from the Arrays page, and the use of -fexcess-precision to speed up the Double computation.

For those interested in my progress (no one), here’s my current code:

{-# OPTIONS -O2
-funbox-strict-fields
-fexcess-precision
-fvia-C
-optc-mfpmath=sse
-optc-msse2#-}
{-# LANGUAGE BangPatterns          #-}
{-# LANGUAGE UnboxedTuples         #-}
{-# LANGUAGE MagicHash             #-}
{-# LANGUAGE MultiParamTypeClasses #-}
module Main where
import System.Environment
import Data.Array.Base
import GHC.Exts
import GHC.ST
import GHC.Arr
import Control.Monad
import Control.Monad.ST
import Data.List
import qualified Data.Array.Unboxed as U
import Data.Array.Base
import Foreign.Storable
import Text.Printf
import Data.STRef
data Vec3 = V {-# UNPACK #-} !Double {-# UNPACK #-} !Double {-# UNPACK #-} !Double  -- deriving (Show)
data Body = B {
pos ::  {-# UNPACK #-} !Vec3,
vel ::  {-# UNPACK #-} !Vec3,
mass :: {-# UNPACK #-} !Double
} -- deriving (Show)
main = do
n <- getArgs >>= \xs -> if null xs then return 1000 else (readIO.head) xs
let (e,e') = runST $ do
ps <- planets
offset ps
e1 <- energy ps
advn n ps
e2 <- energy ps
return (e1,e2)
prtnum e
prtnum e'
advn n ps = replicateM_ n (advance deltat ps)
------------------
offset ps = do
v  <- newSTRef (V 0 0 0)
p  <- readPos ps 0
ve <- readVel ps 0
forM_ index1 $ \i -> do
vel <- readVel ps i
m <- readMass ps i
modifySTRef v (.+. (m . vel))
modifySTRef v (./ solar_mass)
v' <- readSTRef v
writeVel ps 0 (ve .-. v')
advance dt ps = do
forM_ index1 $ \i -> do
-- (B p v m) <- readBody ps i
p <- readPos ps i
v <- readVel ps i
m <- readMass ps i
vvar <- newSTRef v
forM_ (index2 i) $ \j -> do
-- (B p' v' m') <- readBody ps j
p' <- readPos ps j
v' <- readVel ps j
m' <- readMass ps j
let dp = p .-. p'
dst = mag dp
magnit = dt / (dst  dst  dst)
change = magnit . dp
(sca1,sca2) = (m' . change, m . change)
modifySTRef vvar (.-. sca1)
writeVel ps j (v' .+. sca2)
vvar'@(V nx ny nz) <- readSTRef vvar
writeVel' ps i nx ny nz
writePos ps i $ move' dt (B p vvar' m)
index1 = [0..numPs]
-- index1 = map (7) [0..numPs]
index2 i = [(i+1)..numPs]
-- index2 = map (\x -> (x7,x)) [0..numPs]
-- index3 = tail . tails $ index1
energy ps = do
e <- newSTRef (0 :: Double)
forM_ index1 $ \i -> do
b1 <- readBody ps i
modifySTRef e ((+) (kinetic b1))
forM_ (index2 i) $ \j -> do
b2 <- readBody ps j
let x = potential b1 b2
modifySTRef e (subtract x)
readSTRef e
move dt (B p v m) = B (p .+. dt . v) v m
move' dt (B p v m) = p .+. dt . v
-----------------------
writePos :: STUArray s Int Double -> Int -> Vec3  -> ST s ()
writePos !arr !i !(V x y z) =
case i7 of
i'@(I# n) ->
do
unsafeWrite arr i'              x
unsafeWrite arr (I# (n +# 1#))  y
unsafeWrite arr (I# (n +# 2#))  z
writePos' :: STUArray s Int Double -> Int -> Double -> Double -> Double -> ST s ()
writePos' !arr !i !x !y !z =
case i7 of
i'@(I# n) ->
do
unsafeWrite arr i'              x
unsafeWrite arr (I# (n +# 1#))  y
unsafeWrite arr (I# (n +# 2#))  z
writeVel :: STUArray s Int Double -> Int -> Vec3  -> ST s ()
writeVel !arr !i !(V vx vy vz) =
case i7 of
i'@(I# n) ->
do
unsafeWrite arr (I# (n +# 3#)) vx
unsafeWrite arr (I# (n +# 4#)) vy
unsafeWrite arr (I# (n +# 5#)) vz
writeVel' :: STUArray s Int Double -> Int -> Double -> Double -> Double -> ST s ()
writeVel' !arr !i !vx !vy !vz =
case i7 of
i'@(I# n) ->
do
unsafeWrite arr (I# (n +# 3#)) vx
unsafeWrite arr (I# (n +# 4#)) vy
unsafeWrite arr (I# (n +# 5#)) vz
-- writeMass :: STUArray s Int Double -> Int -> Double  -> ST s ()
-- writeMass arr i m =
--     do
--         unsafeWrite arr (i+6) m
writeBody :: STUArray s Int Double -> Int -> Body  -> ST s ()
writeBody !arr !i !(B !(V x y z) !(V vx vy vz) m) =
case i7 of
i'@(I# n) ->
do
unsafeWrite arr (i')            x
unsafeWrite arr (I# (n +# 1#))  y
unsafeWrite arr (I# (n +# 2#))  z
unsafeWrite arr (I# (n +# 3#))  vx
unsafeWrite arr (I# (n +# 4#))  vy
unsafeWrite arr (I# (n +# 5#))  vz
unsafeWrite arr (I# (n +# 6#))  m
readPos :: STUArray s Int Double -> Int -> ST s Vec3
readPos !arr !i =
case i7 of
i'@(I# n) ->
do
x <- unsafeRead arr i'
y <- unsafeRead arr (I# (n +# 1#))
z <- unsafeRead arr (I# (n +# 2#))
return (V x y z)
readVel :: STUArray s Int Double -> Int -> ST s Vec3
readVel !arr !i =
case i7 of
i'@(I# n) ->
do
vx <- unsafeRead arr (I# (n +# 3#))
vy <- unsafeRead arr (I# (n +# 4#))
vz <- unsafeRead arr (I# (n +# 5#))
return (V vx vy vz)
readMass :: STUArray s Int Double -> Int -> ST s Double
readMass !arr !i =
case i7 of
i'@(I# n) -> unsafeRead arr (i7+6)
readBody :: STUArray s Int Double -> Int -> ST s Body
readBody !arr !i =
case i7 of
i'@(I# n) ->
do
x   <- unsafeRead arr i'
y   <- unsafeRead arr (I# (n +# 1#))
z   <- unsafeRead arr (I# (n +# 2#))
vx  <- unsafeRead arr (I# (n +# 3#))
vy  <- unsafeRead arr (I# (n +# 4#))
vz  <- unsafeRead arr (I# (n +# 5#))
m   <- unsafeRead arr (I# (n +# 6#))
return (B (V x y z) (V vx vy vz) m)
-----------------------
solar_mass, deltat, days_per_year :: Double
days_per_year = 365.24
dp = days_per_year
solar_mass    = 4  pi  pi
deltat       = 0.01
numPs       = length planetsList - 1
-- planets'' = newArray (0,(numPs+1)7) 0 :: ST s (STUArray s Int Double)
-- (B (V 0 0 0) (V 0 0 0) 1)
-- ms = U.listArray (0,numPs+1) (map mass planetsList)
planets = do
arr <- newArray (0,(numPs+1)7) 0 :: ST s (STUArray s Int Double)
-- unsafeReplaceUArray arr (zip index1 planetsList)
forM_ (zip index1 planetsList) $ \(i,p) -> do
-- unsafeIOToST (print (i, pos p))
writeBody arr i p
return arr
planetsList =
[B  (V  0.0
0.0
0.0)
(V  0.0
0.0
0.0)
solar_mass,
B  (V  (4.84143144246472090e+00)
(-1.16032004402742839e+00)
(-1.03622044471123109e-01))
(V  (1.66007664274403694e-03dp)
(7.69901118419740425e-03dp)
((-6.90460016972063023e-05)dp))
(9.54791938424326609e-04  solar_mass),
B  (V  (8.34336671824457987e+00)
(4.12479856412430479e+00)
(-4.03523417114321381e-01))
(V  ((-2.76742510726862411e-03)dp)
(4.99852801234917238e-03dp)
(2.30417297573763929e-05dp))
(2.85885980666130812e-04  solar_mass),
B  (V  (1.28943695621391310e+01)
(-1.51111514016986312e+01)
(-2.23307578892655734e-01))
(V  (2.96460137564761618e-03dp)
(2.37847173959480950e-03dp)
((-2.96589568540237556e-05)dp))
(4.36624404335156298e-05  solar_mass),
B  (V  (1.53796971148509165e+01)
(-2.59193146099879641e+01)
(1.79258772950371181e-01))
(V  (2.68067772490389322e-03dp)
(1.62824170038242295e-03dp)
((-9.51592254519715870e-05)dp))
(5.15138902046611451e-05  solar_mass)
]
-----------------------
infix 4 .-.
infix 4 .+.
-- infix 5 ..
infixr 5 .
infixl 5 ./
{-# INLINE (.+.) #-}
{-# INLINE (.-.) #-}
-- {-# INLINE (..) #-}
V x y z .+. V x' y' z' = V (x+x') (y+y') (z+z')
V (D# x) (D# y) (D# z) .-. V (D# x') (D# y') (D# z') =
V (D# (x -## x'))
(D# (y -## y'))
(D# (z -## z'))
-- V x y z .. V x' y' z' = V (xx') (yy') (zz')
{-# INLINE (.) #-}
n . (V x y z) = V (nx) (ny) (nz)
{-# INLINE (./) #-}
V x y z ./ n = V (x/n) (y/n) (z/n)
-- {-# INLINE momentum #-}
momentum :: Body -> Vec3
momentum (B _ v m) = m . v
{-# INLINE kinetic #-}
kinetic x =  0.5  mass x  (mag2.vel) x
-- kinetic'
potential a b =
let dst = dist' a b
in --if dst <= 0
--then 0
--else 
(mass a  mass b) / dist' a b
{-# INLINE mag #-}
mag :: Vec3 -> Double
mag = sqrt . mag2
{-# INLINE mag2 #-}
mag2 (V x y z) = xx + yy + zz
{-# INLINE dist #-}
dist v1 v2 = sqrt $ dist2 v1 v2
{-# INLINE dist2 #-}
dist2 v1 v2 = mag2 $ v1 .-. v2
{-# INLINE dist' #-}
dist' b1 b2 = dist (pos b1) (pos b2)
prtnum :: Double -> IO ()
prtnum x = putStrLn $ printf "%.9f" x
-------------------- tests
thing2 pla = do
arr <- pla
unsafeIOToST (putStrLn "thing2")
forM_ [0..numPs7] $ \i -> do
p <- readArray arr i
unsafeIOToST (print (i,p))
thing3 = do
arr <- planets
offset arr
thing2 planets
unsafeIOToST (putStrLn "\n")
-- thing2 arr
e <- energy arr
return e

STUArray woes
Jan 30th, 2009 by axman6

So a little while ago, I was trying to write a nice, clear version of the language shootout n-body problem, using the ST monad, and STUArrays. Though I ran into a snag; I was getting segmentation faults, for very small values of n (< 23), and not knowing enough about such things, I have no idea where to begin in diagnosing the problem. So, I shall post my code here, and hope that someone might have some ideas. I’ve talked to Don Stewart about it, but he hasn’t got back to me yet. (All the MUArray instance stuff is from Don’s work, which I nicked) So, here we have it:

> {-# OPTIONS -O2 -funbox-strict-fields #-}
> {-# LANGUAGE BangPatterns          #-}
> {-# LANGUAGE UnboxedTuples         #-}
> {-# LANGUAGE MagicHash             #-}
> {-# LANGUAGE MultiParamTypeClasses #-}
>
> module Main where
>
> import System.Environment
> import Data.Array.Base
> import GHC.Exts
> import GHC.ST
> import GHC.Arr
> import Control.Monad
> import Control.Monad.ST
> import Data.List
> import Foreign.Storable
> import Text.Printf
> import Data.STRef
>
>
> data Vec3 = V !Double !Double !Double deriving (Show)
> data Body = B {
>                  pos ::  {-# UNPACK #-} !Vec3,
>                  vel ::  {-# UNPACK #-} !Vec3,
>                  mass :: !Double
>                } deriving (Show)
>
>
> main = do
>     n <- getArgs >>= \xs -> if null xs then return 1000 else (readIO.head) xs
>     let (e,e') = runST $ do
>             ps <- planets
>             offset ps
>             e1 <- energy ps
>             advn n ps
>             e2 <- energy ps
>             return (e1,e2)
>     prtnum e
>     prtnum e'
>
>
> offset ps = do
>     v <- newSTRef (V 0 0 0)
>     (B p ve m) <- readArray ps 0
>     forM_ [0..numPs] $ \i -> do
>         p <- readArray ps i
>         modifySTRef v (.+. momentum p)
>     modifySTRef v (./ solar_mass)
>     v' <- readSTRef v
>     writeArray ps 0 (B p (ve.-. v') m)
>
>
>
> advn 0 ps = do
>     return ps
> advn n ps = do
>     advance deltat ps
>     -- e <- energy ps
>     -- unsafeIOToST (prtnum e)
>     advn (n-1) ps
>
> advance dt ps = do
>     forM_ [0..numPs] $ \i -> do
>         (B p v m) <- readArray ps i
>         vvar <- newSTRef v
>         forM_ [i+1..numPs] $ \j -> do
>             (B p' v' m') <- readArray ps j
>             let dp = p .-. p'
>                 dst = mag dp
>                 magnit = dt / (dst  dst  dst)
>                 change = magnit . dp
>                 sca1 = m' . change
>                 sca2 = m . change
>             -- writeArray ps I (B p (v .-. sca1) m)
>             modifySTRef vvar (.-. sca1)
>             writeArray ps j (B p' (v' .+. sca2) m')
>         vvar' <- readSTRef vvar
>         writeArray ps I (move dt (B p vvar' m))
>
>
> energy ps = do
>     e <- newSTRef (0 :: Double)
>
>     forM_ [0..numPs] $ \i -> do
>         b1 <- readArray ps i
>         modifySTRef e ((+) (kinetic b)1)
>
>         forM_ [i+1..numPs] $ \j -> do
>             b2 <- readArray ps j
>             let x = potential b1 b2
>             modifySTRef e (subtract x)
>     e' <- readSTRef e
>     return e'
>
>
>
> move dt (B p v m) = B (p .+. dt . v) v m
>
>
>
>
> -----------------------
> solar_mass, deltat, days_per_year :: Double
> days_per_year = 365.24
> dp = days_per_year
> solar_mass    = 4  pi  pi
> deltat       = 0.01
> numPs       = length planets' - 1
>
> planets'' = newArray_ (0,numPs) :: ST s (STUArray s Int Body)
>     -- (B (V 0 0 0) (V 0 0 0) 1)
>
> planets = do
>     arr <- planets''
>     forM_ (zip [0..] planets') $ \(i,p) -> do
>         writeArray arr I p
>     return arr
>
> planets' =
>     [B  (V  0.0
>             0.0
>             0.0)
>         (V  0.0
>             0.0
>             0.0)
>         solar_mass,
>
>      B  (V  4.84143144246472090e+00
>             (-1.16032004402742839e+00)
>             (-1.03622044471123109e-01))
>         (V  (1.66007664274403694e-03dp)
>             (7.69901118419740425e-03dp)
>             ((-6.90460016972063023e-05)dp))
>         (9.54791938424326609e-04  solar_mass),
>
>      B  (V  8.34336671824457987e+00
>             4.12479856412430479e+00
>             (-4.03523417114321381e-01))
>         (V  ((-2.76742510726862411e-03)dp)
>             (4.99852801234917238e-03dp)
>             (2.30417297573763929e-05dp))
>         (2.85885980666130812e-04  solar_mass),
>
>      B  (V  1.28943695621391310e+01
>             (-1.51111514016986312e+01)
>             (-2.23307578892655734e-01))
>         (V  (2.96460137564761618e-03dp)
>             (2.37847173959480950e-03dp)
>             ((-2.96589568540237556e-05)dp))
>         (4.36624404335156298e-05  solar_mass),
>
>      B  (V  1.53796971148509165e+01
>             (-2.59193146099879641e+01)
>             1.79258772950371181e-01)
>         (V  (2.68067772490389322e-03dp)
>             (1.62824170038242295e-03dp)
>             ((-9.51592254519715870e-05)dp))
>         (5.15138902046611451e-05  solar_mass)
>     ]
>
>
> -----------------------
>
> infix 4 .-.
> infix 4 .+.
> -- infix 5 ..
> infixr 5 .
> infixl 5 ./
>
> {-# INLINE (.+.) #-}
> {-# INLINE (.-.) #-}
> -- {-# INLINE (..) #-}
> V x y z .+. V x' y' z' = V (x+x') (y+y') (z+z')
> V x y z .-. V x' y' z' = V (x-x') (y-y') (z-z')
> -- V x y z .. V x' y' z' = V (xx') (yy') (zz')
>
> {-# INLINE (.) #-}
> n . (V x y z) = V (nx) (ny) (nz)
> {-# INLINE (./) #-}
> V x y z ./ n = V (x/n) (y/n) (z/n)
>
> {-# INLINE momentum #-}
> momentum :: Body -> Vec3
> momentum (B _ v m) = m . v
>
> {-# INLINE kinetic #-}
> kinetic x =  0.5  mass x  (mag2.vel) x
>
>
> potential a b =
>     let dst = dist' a b
>     in if dst <= 0
>         then 0
>         else (mass a  mass b) / dist' a b
>
> {-# INLINE mag #-}
> mag :: Vec3 -> Double
> mag = sqrt . mag2
>
> {-# INLINE mag2 #-}
> mag2 (V x y z) = xx + yy + zz
>
> {-# INLINE dist #-}
> dist v1 v2 = sqrt $ dist2 v1 v2
>
> {-# INLINE dist2 #-}
> dist2 v1 v2 = mag2 $ v1 .-. v2
>
> {-# INLINE dist' #-}
> dist' b1 b2 = dist (pos b1) (pos b2)
>
> prtnum x = putStrLn $ printf "%.9f" x
>
> --------------------
> bodyScale :: Int# -> Int#
> bodyScale n = (scale # 7#) *# n
>      where I# scale = sizeOf (undefined :: Double)
>
> instance MArray (STUArray s) Body (ST s) where
>
>     {-# INLINE getBounds #-}
>     getBounds (STUArray l h _ _) = return (l,h)
>
>     {-# INLINE getNumElements #-}
>     getNumElements (STUArray _ _ n _) = return n
>
>     {-# INLINE unsafeNewArray_ #-}
>     unsafeNewArray_ (l,u) = unsafeNewArraySTUArray_ (l,u) bodyScale -- bytes 
>
>     {-# INLINE newArray_ #-}
>     newArray_ bounds = newArray bounds (B (V 0 0 0) (V 0 0 0) 1)
>
>
>     {-# INLINE unsafeRead #-}
>     unsafeRead (STUArray _ _ _ arr) (I# i) = ST $ \s1 ->
>         let j = bodyScale i in
>         case readDoubleArray# arr j         s1 of { (# s2, xd    #) ->
>         case readDoubleArray# arr (j +# 1#) s2 of { (# s3, yd    #) ->
>         case readDoubleArray# arr (j +# 2#) s3 of { (# s4, zd    #) ->
>         case readDoubleArray# arr (j +# 3#) s4 of { (# s5, vxd   #) ->
>         case readDoubleArray# arr (j +# 4#) s5 of { (# s6, vyd   #) ->
>         case readDoubleArray# arr (j +# 5#) s6 of { (# s7, vzd   #) ->
>         case readDoubleArray# arr (j +# 6#) s7 of { (# s8, massd #) ->
>         (# s8, B { pos    = V (D# xd) (D# yd) (D# zd)
>                     , vel    = V (D# vxd) (D# vyd) (D# vzd)
>                     , mass = D# massd } #) }}}}}}}
>
>     {-# INLINE unsafeWrite #-}
>     unsafeWrite (STUArray _ _ _ arr) (I# i)
>                 (B (V (D# xd) (D# yd) (D# zd)) (V (D# vxd) (D# vyd) (D# vzd)) (D# massd))
>         = ST $ \s1 ->
>         let j = bodyScale i in
>         case writeDoubleArray# arr j          xd   s1  of {s2 ->
>         case writeDoubleArray# arr (j +# 1#)  yd   s2  of {s3 ->
>         case writeDoubleArray# arr (j +# 2#)  zd   s3  of {s4 ->
>         case writeDoubleArray# arr (j +# 3#) vxd   s4  of {s5 ->
>         case writeDoubleArray# arr (j +# 4#) vyd   s5  of {s6 ->
>         case writeDoubleArray# arr (j +# 5#) vzd   s6  of {s7 ->
>         case writeDoubleArray# arr (j +# 6#) massd s7  of {s8 ->
>         (# s8, () #) }}}}}}}
> 

I just finished installing Movable Type 4!
Jan 24th, 2009 by axman6

Welcome to my new blog powered by Movable Type. This is the first post on my blog and was created for me automatically when I finished the installation process. But that is ok, because I will soon be creating posts of my own!

»  Substance: WordPress   »  Style: Ahren Ahimsa
© Alex Mason (Axman6) 2009