Re: [Haskell-cafe] my Fasta is slow ;(
I've already submitted it, thanks. The Fortran program commits the same sin as the C++ one, of doing floating point arithmetic in the inner loop; that's why it's slow. On Dec 27, 2012, at 18:05, Branimir Maksimovic bm...@hotmail.com wrote: Thank you. Your entry is great. Faster than fortran entry! Dou you want to contribute at the site, or you want me to do it for you? -- Date: Thu, 27 Dec 2012 15:58:40 -0800 Subject: Re: [Haskell-cafe] my Fasta is slow ;( From: b...@serpentine.com To: bm...@hotmail.com CC: haskell-cafe@haskell.org On Tue, Dec 18, 2012 at 12:42 PM, Branimir Maksimovic bm...@hotmail.comwrote: Seems to me that culprit is in function random as I have tested rest of code and didn't found speed related problems. The problem with your original program was that it was not pure enough. Because you stored your PRNG state in an IORef, you forced the program to allocate and case-inspect boxed Ints in its inner loop. I refactored it slightly to make genRand and genRandom pure, and combined with using the LLVM back end, this doubled the program's performance, so that the Haskell program ran at the same speed as your C++ version. The next bottleneck was that your program was performing floating point arithmetic in the inner loop. I changed it to precompute a small lookup table, followed by only using integer arithmetic in the inner loop (the same technique used by the fastest C fasta program). This further improved performance: the new Haskell code is 40% faster than the C++ program, and only ~20% slower than the C program that currently tops the shootout chart. The Haskell source is a little over half the size of the C source. You can follow the work I did here: https://github.com/bos/shootout-fasta ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] my Fasta is slow ;(
On Tue, Dec 18, 2012 at 12:42 PM, Branimir Maksimovic bm...@hotmail.comwrote: Seems to me that culprit is in function random as I have tested rest of code and didn't found speed related problems. The problem with your original program was that it was not pure enough. Because you stored your PRNG state in an IORef, you forced the program to allocate and case-inspect boxed Ints in its inner loop. I refactored it slightly to make genRand and genRandom pure, and combined with using the LLVM back end, this doubled the program's performance, so that the Haskell program ran at the same speed as your C++ version. The next bottleneck was that your program was performing floating point arithmetic in the inner loop. I changed it to precompute a small lookup table, followed by only using integer arithmetic in the inner loop (the same technique used by the fastest C fasta program). This further improved performance: the new Haskell code is 40% faster than the C++ program, and only ~20% slower than the C program that currently tops the shootout chart. The Haskell source is a little over half the size of the C source. You can follow the work I did here: https://github.com/bos/shootout-fasta ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] my Fasta is slow ;(
Thank you. Your entry is great. Faster than fortran entry!Dou you want to contribute at the site, or you want me to do it for you? Date: Thu, 27 Dec 2012 15:58:40 -0800 Subject: Re: [Haskell-cafe] my Fasta is slow ;( From: b...@serpentine.com To: bm...@hotmail.com CC: haskell-cafe@haskell.org On Tue, Dec 18, 2012 at 12:42 PM, Branimir Maksimovic bm...@hotmail.com wrote: Seems to me that culprit is in function random as I have tested rest of codeand didn't found speed related problems. The problem with your original program was that it was not pure enough. Because you stored your PRNG state in an IORef, you forced the program to allocate and case-inspect boxed Ints in its inner loop. I refactored it slightly to make genRand and genRandom pure, and combined with using the LLVM back end, this doubled the program's performance, so that the Haskell program ran at the same speed as your C++ version. The next bottleneck was that your program was performing floating point arithmetic in the inner loop. I changed it to precompute a small lookup table, followed by only using integer arithmetic in the inner loop (the same technique used by the fastest C fasta program). This further improved performance: the new Haskell code is 40% faster than the C++ program, and only ~20% slower than the C program that currently tops the shootout chart. The Haskell source is a little over half the size of the C source. You can follow the work I did here: https://github.com/bos/shootout-fasta ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] my Fasta is slow ;(
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 bm...@hotmail.comwrote:
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=fastalang=gppid=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 2500 /dev/null
real0m5.262s
user0m5.228s
sys 0m0.020s
bmaxa@maxa:~/shootout/fasta$ time ./fastacpp 2500 /dev/null
real0m2.075s
user0m2.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\
\CCAGCCTGGCCAACATGGTGAAAGTCTCTACTAT\
\ACATTAGCCGGGCGTGGTGGCGCGCGCCTGTAATCCCA\
\GCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGG\
\AGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTGCACTCC\
\AGCCTGGGCGACAGAGCGAGACTCCGTCTCA
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)
[Haskell-cafe] my Fasta is slow ;(
This time I have tried fasta benchmark since current entries does notdisplay
correct output.Program is copy of mine
http://benchmarksgame.alioth.debian.org/u64q/program.php?test=fastalang=gppid=1c++
benchmark, but unfortunately executes more than twice time.
Seems to me that culprit is in function random as I have tested rest of
codeand didn't found speed related problems.
bmaxa@maxa:~/shootout/fasta$ time ./fastahs 2500 /dev/null
real0m5.262suser0m5.228ssys 0m0.020s
bmaxa@maxa:~/shootout/fasta$ time ./fastacpp 2500 /dev/null
real0m2.075suser0m2.056ssys 0m0.012s
Since I am planning to contribute program, perhaps someone cansee 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.Environmentimport System.IO.Unsafe
import Data.IORefimport Data.Array.Unboxedimport Data.Array.Storableimport
Data.Array.Baseimport Data.Word
import Foreign.Ptrimport Foreign.C.Types
type A = UArray Int Word8type B = StorableArray Int Word8type 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 = don - 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 iubmake
THREE Homo sapiens frequency (n*5) $ random homosapiens
make :: String - String - Int - IO Word8 - IO (){-# INLINE make #-}make id
desc n f = dolet lst = ++ id ++ ++ desca - (newListArray
(0,length lst) $ map (fromIntegral. fromEnum) lst:: IO B)
unsafeWrite a (length lst) 0pr amake' n 0where make' :: Int
- Int - IO ()make' !n !i = dolet line = (unsafePerformIO
$ newArray (0,60) 0 :: B)if n 0
then do!c - funsafeWrite line i c
if i+1 = 60 then do
pr linemake' (n-1) 0
else make' (n-1) (i+1)else do
unsafeWrite line i 0l - len line
if l /= 0then pr line
else return ()
pr :: B - IO ()pr line = withStorableArray line (\ptr - puts ptr)len :: B -
IO CIntlen line = withStorableArray line (\ptr - strlen ptr)
repeat :: A - Int - IO Word8repeat xs !n = dolet v = unsafePerformIO
$ newIORef 0!i - readIORef vif i+1 = nthen
writeIORef v 0else writeIORef v (i+1)return $ xs `unsafeAt`
i
random :: C - IO Word8random (a,b) = do !rnd - randlet
find :: Int - IO Word8find !i = let
!c = a `unsafeAt` i!p = b `unsafeAt` i
in if p = rndthen return celse
find (i+1)find 0
rand :: IO Double{-# INLINE rand #-}rand = do!seed - readIORef lastlet
newseed = (seed * ia + ic) `rem` imnewran = fromIntegral
newseed * rimdrimd = 1.0 / (fromIntegral im)im, ia, ic ::
Intim = 139968ia = 3877ic = 29573writeIORef last
newseedreturn newranwhere last = unsafePerformIO $ newIORef 42
alu:: [Char]alu = GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGG\
\GAGGCCGAGGCGGGCGGATCACCTGAGGTCAGGAGTTCGAGA\
\CCAGCCTGGCCAACATGGTGAAAGTCTCTACTAT\
\ACATTAGCCGGGCGTGGTGGCGCGCGCCTGTAATCCCA\
\GCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGG\
\AGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTGCACTCC\
\AGCCTGGGCGACAGAGCGAGACTCCGTCTCA
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')
