Commit 806920d4 authored by tibbe's avatar tibbe

Add the k-nucleotide shootout benchmark

parent 7c56e6fd
......@@ -58,6 +58,11 @@ shootout/fasta/fasta-c
shootout/fasta/fasta.faststdout
shootout/fasta/fasta.slowstdout
shootout/fasta/fasta.stdout
shootout/k-nucleotide/fasta-c
shootout/k-nucleotide/k-nucleotide
shootout/k-nucleotide/knucleotide-input250000.txt
shootout/k-nucleotide/knucleotide-input2500000.txt
shootout/k-nucleotide/knucleotide-input25000000.txt
shootout/n-body/n-body
shootout/pidigits/pidigits
shootout/reverse-complement/fasta-c
......
......@@ -15,7 +15,4 @@ Notes:
which requires C libraries that aren't installed by default on many
platforms.
* The k-nucleotide benchmark wasn't included as it segfaults on my
64-bit OS X 10.8.2 machine.
1. http://benchmarksgame.alioth.debian.org/
--
-- The Computer Language Benchmarks Game
-- http://benchmarksgame.alioth.debian.org/
--
-- contributed by Stephen Blackheath (with some bits taken from Don Stewart's
-- version), v1.2
import Text.Printf
import Data.ByteString.Internal
import qualified Data.ByteString.Char8 as S
import Control.Applicative
import Control.Monad
import Control.Concurrent
import Foreign.Storable
import Foreign.Marshal.Alloc
import Foreign.Marshal.Array
import Foreign.Ptr
import Foreign.ForeignPtr
import Data.Word
import Data.Bits
import Data.Char
import Data.List
import Data.Maybe
import Data.IORef
import GHC.Exts
main = do
genome <- extract (S.pack ">TH")
let actions = [
do
a <- printFreqsBySize genome 1
b <- printFreqsBySize genome 2
return $ a ++ b
] ++ map (printFreqsSpecific genome) specificSeqs
output <- concat <$> parallel actions
forM_ output putStrLn
-- Drop in replacement for sequence
parallel :: [IO a] -> IO [a]
parallel actions = do
vars <- forM actions $ \action -> do
var <- newEmptyMVar
forkIO $ do
answer <- action
putMVar var answer
return var
forM vars takeMVar
specificSeqs = map S.pack [
"GGT","GGTA","GGTATT","GGTATTTTAATT","GGTATTTTAATTTATAGT"]
extract p = do
s <- S.getContents
let (_, rem) = S.breakSubstring p s
return $! S.map toUpper -- array fusion!
. S.filter ((/=) '\n')
. S.dropWhile ((/=) '\n') $ rem
printFreqsBySize :: S.ByteString -> Int -> IO [String]
printFreqsBySize genome keySize = do
ht0 <- htNew keySize
ht <- hashGenome genome keySize ht0
l <- htToList ht
htFree ht
return $ map draw (sortBy sortRule l) ++ [""]
where
genomeLen = S.length genome
draw :: (S.ByteString, Int) -> String
draw (key, count) = printf "%s %.3f" (S.unpack key) pct
where pct = (100 * (fromIntegral count) / total) :: Double
total = fromIntegral (genomeLen - keySize + 1)
printFreqsSpecific :: S.ByteString -> S.ByteString -> IO [String]
printFreqsSpecific genome seq = do
let keySize = S.length seq
genomeLen = S.length genome
ht0 <- htNew keySize
ht <- hashGenome genome keySize ht0
let (fp, offset, len) = toForeignPtr seq
count <- withForeignPtr fp $ \p_ -> do
htGet ht (p_ `plusPtr` offset)
htFree ht
return [show count ++ ('\t' : S.unpack seq)]
hashGenome :: S.ByteString
-> Int
-> Hashtable
-> IO Hashtable
{-# INLINE hashGenome #-}
hashGenome genome keySize ht = do
let (fp, offset, len) = toForeignPtr genome
withForeignPtr fp $ \p_ -> do
let p = p_ `plusPtr` offset
loop ht idx = do
let key = p `plusPtr` idx
htInc ht key
endIdx = len - keySize + 1
foldMe i ht | i == endIdx = return ht
foldMe i ht = loop ht i >>= foldMe (i+1)
foldMe 0 ht
sortRule :: (S.ByteString, Int) -> (S.ByteString, Int) -> Ordering
sortRule (a1, b1) (a2, b2) =
case compare b2 b1 of
EQ -> compare a1 a2
x -> x
------ Hash table implementation ----------------------------------------------
-- Note: Hash tables are not generally used in functional languages, so there
-- are no available library implementations for Haskell. This benchmark
-- requires a hash table. This is why I have implemented the hash table here.
htNew :: Int -> IO Hashtable
htNew = htNew' (head primes)
htNew' :: Int -> Int -> IO Hashtable
htNew' slots ksz = do
let ssz = spineSize ksz slots
sp <- mallocBytes ssz
memset sp 0 (fromIntegral ssz)
return $ Hashtable {
keySize = ksz,
noOfSlots = slots,
spine = sp
}
primes = [ 1543, 3079, 6151, 12289, 24593,
49157, 98317, 196613, 93241, 786433,
1572869, 3145739, 6291469, 12582917, 25165843,
50331653, 100663319, 201326611, 402653189, 805306457 ]
htFree :: Hashtable -> IO ()
htFree ht = do
htTraverse ht $ \isSpine (Entry ePtr) -> when (not isSpine) $ free ePtr
free (spine ht)
htGet :: Hashtable -> Ptr Word8 -> IO Int
htGet ht key = do
hash <- calcHash ht key
htPayload ht hash key >>= peek
htInc :: Hashtable -> Ptr Word8 -> IO Hashtable
{-# INLINE htInc #-}
htInc !ht !key = do
hash <- calcHash ht key
pPayload <- htPayload ht hash key
value <- peek pPayload
poke pPayload (value+1)
if value == 0 && (hash .&. 0x7f) == 0
then checkGrow ht
else return ht
htPut_ :: Hashtable -> Ptr Word8 -> Int -> IO ()
{-# INLINE htPut_ #-}
htPut_ !ht !key !value = do
hash <- calcHash ht key
pPayload <- htPayload ht hash key
poke pPayload value
checkGrow :: Hashtable -> IO Hashtable
checkGrow ht = do
let pTotal = totalEntriesOf ht
slots = noOfSlots ht
total <- (0x200+) <$> peek pTotal
poke pTotal total
if total >= slots
then do
let newSlots = head $ dropWhile (<= slots*2) primes
ht' <- htNew' newSlots (keySize ht)
htTraverse ht $ \_ -> reinsert ht' -- re-insert all the elts
htFree ht
poke (totalEntriesOf ht') total -- copy the total entry count
return ht'
else return ht
where
reinsert :: Hashtable -> Entry -> IO ()
reinsert ht entry = do
value <- peek (payloadOf entry)
htPut_ ht (keyOf entry) value
htToList :: Hashtable -> IO [(S.ByteString, Int)]
htToList ht =
htMap (\entry -> do
keyStr <- keyString ht (keyOf entry)
payload <- peek (payloadOf entry)
return (keyStr, payload)) ht
htMap :: (Entry -> IO a) -> Hashtable -> IO [a]
htMap f ht = mapM f =<< htEntries ht
keyString :: Hashtable -> Ptr Word8 -> IO S.ByteString
keyString ht key = S.pack . map w2c <$> peekArray (keySize ht) key
isEmptySlot :: Entry -> IO Bool
{-# INLINE isEmptySlot #-}
isEmptySlot entry = do
ch <- peek $ keyOf entry
return $ ch == 0
htEntries :: Hashtable -> IO [Entry]
htEntries ht = do
es <- newIORef []
htTraverse ht $ \_ entry -> modifyIORef es $ \l -> entry:l
readIORef es
htTraverse :: Hashtable -> (Bool -> Entry -> IO ()) -> IO ()
htTraverse ht f = he 0
where
slots = noOfSlots ht
he i | i == slots = return ()
he i = do
let entry = indexEntry ht i
empty <- isEmptySlot entry
if empty
then he (i+1)
else links True i entry
links isSpine i entry = do
next <- peek $ nextPtrOf entry
f isSpine entry
if next == nullPtr
then he (i+1)
else links False i (Entry next)
data Hashtable = Hashtable {
keySize :: Int,
noOfSlots :: Int,
spine :: Ptr Word8
}
wordSize :: Int
wordSize = max (sizeOf (nullPtr :: Ptr Word8)) (sizeOf (0 :: Int))
-- Round up to word size
roundUp :: Int -> Int
{-# INLINE roundUp #-}
roundUp !i = (i + wordSize - 1) .&. complement (wordSize - 1)
slotSize :: Int -> Int
{-# INLINE slotSize #-}
slotSize !ksz = roundUp ksz + wordSize * 2
spineSize :: Int -> Int -> Int
spineSize ksz slots = slots * slotSize ksz + wordSize
calcHash :: Hashtable -> Ptr Word8 -> IO Int
{-# INLINE calcHash #-}
calcHash !ht !key = (`mod` noOfSlots ht) <$> ch 0 0
where
ksz = keySize ht
ch :: Int -> Int -> IO Int
ch !i !acc | i == ksz = return acc
ch !i !acc = do
c <- peek (key `plusPtr` i)
ch (i+1) (acc * 131 + fromIntegral (c::Word8))
newtype Entry = Entry (Ptr Word8)
-- Count of the total number of hash table entries
totalEntriesOf :: Hashtable -> Ptr Int
{-# INLINE totalEntriesOf #-}
totalEntriesOf ht = castPtr $ spine ht
indexEntry :: Hashtable -> Int -> Entry
{-# INLINE indexEntry #-}
indexEntry !ht !hash =
let hOffset = wordSize + hash * slotSize (keySize ht)
in Entry $ spine ht `plusPtr` hOffset
entryMatches :: Hashtable -> Entry -> Ptr Word8 -> IO Bool
{-# INLINE entryMatches #-}
entryMatches !ht !entry !key = do
let eKey = keyOf entry
c <- memcmp key eKey (fromIntegral $ keySize ht)
if c == 0
then return True
else do
empty <- isEmptySlot entry
if empty
then do
memcpy eKey key (fromIntegral $ keySize ht) -- ick
return True
else
return False
nextPtrOf :: Entry -> Ptr (Ptr Word8)
{-# INLINE nextPtrOf #-}
nextPtrOf !(Entry ePtr) = castPtr $ ePtr
payloadOf :: Entry -> Ptr Int
{-# INLINE payloadOf #-}
payloadOf !(Entry ePtr) = castPtr $ ePtr `plusPtr` wordSize
keyOf :: Entry -> Ptr Word8
{-# INLINE keyOf #-}
keyOf !(Entry ePtr) = ePtr `plusPtr` (wordSize*2)
allocEntry :: Hashtable -> Ptr Word8 -> IO Entry
allocEntry !ht !key = do
let esz = slotSize $ keySize ht
ePtr <- mallocBytes esz
memset ePtr 0 (fromIntegral esz)
let entry = Entry ePtr
memcpy (keyOf entry) key (fromIntegral $ keySize ht)
return entry
htPayload :: Hashtable -> Int -> Ptr Word8 -> IO (Ptr Int)
{-# INLINE htPayload #-}
htPayload !ht !hash !key = do
entry <- findEntry (indexEntry ht hash)
return $ payloadOf entry
where
findEntry :: Entry -> IO Entry
findEntry !entry = do
match <- entryMatches ht entry key
if match
then
return entry
else do
let pNext = nextPtrOf entry
next <- peek pNext
if next == nullPtr
then do
newEntry@(Entry ePtr) <- allocEntry ht key
poke pNext ePtr
return newEntry
else
findEntry (Entry next)
TOP = ../..
include $(TOP)/mk/boilerplate.mk
# These values are only used in this file. They are ignored by the
# executable itself.
FAST_OPTS = 250000
NORM_OPTS = 2500000
SLOW_OPTS = 25000000 # official shootout setting
# The benchmark game also uses -fllvm, which we can't since it might
# not be available on the developer's machine.
HC_OPTS += -O2 -XBangPatterns -package bytestring
#------------------------------------------------------------------
# Create input
fasta-c : ../fasta/fasta-c.o
gcc $< -o $@
knucleotide-input250000.txt : fasta-c
./fasta-c $(FAST_OPTS) > $@
knucleotide-input2500000.txt : fasta-c
./fasta-c $(NORM_OPTS) > $@
knucleotide-input25000000.txt : fasta-c
./fasta-c $(SLOW_OPTS) > $@
ifeq "$(mode)" "slow"
INPUT_FILE = knucleotide-input25000000.txt
else
ifeq "$(mode)" "fast"
INPUT_FILE = knucleotide-input250000.txt
else
INPUT_FILE = knucleotide-input2500000.txt
endif
endif
SRC_RUNTEST_OPTS += -i $(INPUT_FILE)
all boot :: $(INPUT_FILE)
include $(TOP)/mk/target.mk
A 30.298
T 30.157
C 19.793
G 19.752
AA 9.177
TA 9.137
AT 9.136
TT 9.094
AC 6.000
CA 5.999
GA 5.986
AG 5.985
TC 5.970
CT 5.970
GT 5.957
TG 5.956
CC 3.915
CG 3.910
GC 3.908
GG 3.902
14717 GGT
4463 GGTA
472 GGTATT
9 GGTATTTTAATT
9 GGTATTTTAATTTATAGT
A 30.295
T 30.151
C 19.800
G 19.754
AA 9.177
TA 9.132
AT 9.131
TT 9.091
CA 6.002
AC 6.001
AG 5.987
GA 5.984
CT 5.971
TC 5.971
GT 5.957
TG 5.956
CC 3.917
GC 3.911
CG 3.909
GG 3.902
1471758 GGT
446535 GGTA
47336 GGTATT
893 GGTATTTTAATT
893 GGTATTTTAATTTATAGT
A 30.297
T 30.151
C 19.798
G 19.755
AA 9.177
TA 9.133
AT 9.131
TT 9.091
CA 6.002
AC 6.001
AG 5.987
GA 5.984
CT 5.971
TC 5.971
GT 5.957
TG 5.956
CC 3.917
GC 3.910
CG 3.909
GG 3.903
147166 GGT
44658 GGTA
4736 GGTATT
89 GGTATTTTAATT
89 GGTATTTTAATTTATAGT
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment