[Haskell-cafe] my Fasta is slow ;(

[Bryan O'Sullivan]

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 <bmaxa at hotmail.com>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 at maxa:~/shootout/fasta$ time ./fastahs 25000000 > /dev/null
>
> real    0m5.262s
> user    0m5.228s
> sys     0m0.020s
>
> bmaxa at 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 at haskell.org
> http://www.haskell.org/mailman/listinfo/haskell-cafe
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.haskell.org/pipermail/haskell-cafe/attachments/20121219/d35adbf8/attachment-0001.htm>