Re: [Haskell-cafe] my Fasta is slow ;(

19 Dec 2012


      I took your Haskell program as a base and have refactored it into a version
that is about the same speed as your original C++ program. Will follow up
with details when I have a little more time.

On Tue, Dec 18, 2012 at 12:42 PM, Branimir Maksimovic wrote:
...
This time I have tried fasta benchmark since current entries does not
display correct output.
Program is copy of mine
http://benchmarksgame.alioth.debian.org/u64q/program.php?test=fasta&lang=gpp&id=1
c++ benchmark, but unfortunately executes more than twice time.
Seems to me that culprit  is in function random as I have tested rest of
code
and didn't found speed related  problems.
bmaxa@maxa:~/shootout/fasta$ time ./fastahs 25000000 > /dev/null
real    0m5.262s
user    0m5.228s
sys     0m0.020s
bmaxa@maxa:~/shootout/fasta$ time ./fastacpp 25000000 > /dev/null
real    0m2.075s
user    0m2.056s
sys     0m0.012s
Since I am planning to contribute program, perhaps someone can
see a problem to speed it up at least around 3.5 secs which is
speed of bench that display incorrect result  (in 7.6.1).
Program follows:
{-# LANGUAGE BangPatterns #-}
{-  The Computer Language Benchmarks Game
http://shootout.alioth.debian.org/
contributed by Branimir Maksimovic
-}
import System.Environment
import System.IO.Unsafe
import Data.IORef
import Data.Array.Unboxed
import Data.Array.Storable
import Data.Array.Base
import Data.Word
import Foreign.Ptr
import Foreign.C.Types
type A = UArray Int Word8
type B = StorableArray Int Word8
type C = (UArray Int Word8,UArray Int Double)
foreign import ccall unsafe "stdio.h"
     puts  :: Ptr a -> IO ()
foreign import ccall unsafe "string.h"
     strlen :: Ptr a -> IO CInt
main :: IO ()
main = do
    n <- getArgs >>= readIO.head
let !a = (listArray (0,(length alu)-1)
             $ map (fromIntegral. fromEnum) alu:: A)
    make "ONE" "Homo sapiens alu" (n*2) $ Main.repeat a (length alu)
    make "TWO"  "IUB ambiguity codes" (n*3) $ random iub
    make "THREE" "Homo sapiens frequency" (n*5) $ random homosapiens
make :: String -> String -> Int -> IO Word8 -> IO ()
{-# INLINE make #-}
make id desc n f = do
    let lst = ">" ++ id ++ " " ++ desc
    a <- (newListArray (0,length lst)
        $ map (fromIntegral. fromEnum) lst:: IO B)
    unsafeWrite a (length lst) 0
    pr a
    make' n 0
    where
        make' :: Int -> Int -> IO ()
        make' !n !i = do
            let line = (unsafePerformIO $
                        newArray (0,60) 0 :: B)
            if n > 0
                then do
                    !c <- f
                    unsafeWrite line i c
                    if i+1 >= 60
                        then do
                            pr line
                            make' (n-1) 0
                        else
                            make' (n-1) (i+1)
                else do
                    unsafeWrite line i 0
                    l <- len line
                    if l /= 0
                        then pr line
                        else return ()
pr :: B -> IO ()
pr line = withStorableArray line (\ptr -> puts ptr)
len :: B -> IO CInt
len line  = withStorableArray line (\ptr -> strlen ptr)
repeat :: A -> Int -> IO Word8
repeat xs !n = do
        let v = unsafePerformIO $ newIORef 0
        !i <- readIORef v
        if i+1 >= n
            then writeIORef v 0
            else writeIORef v (i+1)
        return $ xs `unsafeAt` i
random :: C -> IO Word8
random (a,b) = do
        !rnd <- rand
        let
            find :: Int -> IO Word8
            find !i =
                let
                    !c = a `unsafeAt` i
                    !p = b `unsafeAt` i
                in if p >= rnd
                    then return c
                    else find (i+1)
        find 0
rand :: IO Double
{-# INLINE rand #-}
rand = do
    !seed <- readIORef last
    let
        newseed = (seed * ia + ic) `rem` im
        newran  =  fromIntegral newseed * rimd
        rimd      = 1.0 / (fromIntegral im)
        im, ia, ic :: Int
        im  = 139968
        ia  = 3877
        ic  = 29573
    writeIORef last newseed
    return newran
    where
        last = unsafePerformIO $ newIORef 42
alu    :: [Char]
alu =
    "GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGG\
    \GAGGCCGAGGCGGGCGGATCACCTGAGGTCAGGAGTTCGAGA\
    \CCAGCCTGGCCAACATGGTGAAACCCCGTCTCTACTAAAAAT\
    \ACAAAAATTAGCCGGGCGTGGTGGCGCGCGCCTGTAATCCCA\
    \GCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGG\
    \AGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTGCACTCC\
    \AGCCTGGGCGACAGAGCGAGACTCCGTCTCAAAAA"
mkCum :: [(Char,Double)] -> [(Word8,Double)]
mkCum lst = map (\(c,p) -> ((fromIntegral.fromEnum) c,p)) $
              scanl1 (\(_,p) (c',p') -> (c', p+p')) lst
homosapiens, iub :: C
iub' = mkCum [('a',0.27),('c',0.12),('g',0.12),('t',0.27),('B',0.02)
        ,('D',0.02),('H',0.02),('K',0.02),('M',0.02),('N',0.02)
        ,('R',0.02),('S',0.02),('V',0.02),('W',0.02),('Y',0.02)]
homosapiens' = mkCum [('a',0.3029549426680),('c',0.1979883004921)
                ,('g',0.1975473066391),('t',0.3015094502008)]
iub = (listArray (0, (length iub')-1) $ map fst iub',
        listArray (0, (length iub')-1) $ map snd iub')
homosapiens = (listArray (0, (length homosapiens')-1) $ map fst
homosapiens',
                listArray (0, (length homosapiens')-1) $ map snd
homosapiens')
_______________________________________________
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe

Re: [Haskell-cafe] my Fasta is slow ;(

Bryan O'Sullivan