SuffixStructures (empty) → 0.0.1.0
raw patch · 9 files changed
+482/−0 lines, 9 filesdep +ListLikedep +SuffixStructuresdep +aesonsetup-changed
Dependencies added: ListLike, SuffixStructures, aeson, base, binary, bytestring, cereal, cmdargs, containers, criterion, data-default-class, deepseq, mwc-random, primitive, vector, vector-algorithms, vector-binary-instances
Files
- Data/SuffixStructure/ESA.hs +64/−0
- Data/SuffixStructure/NaiveArray.hs +80/−0
- LICENSE +30/−0
- README.md +17/−0
- Setup.hs +2/−0
- SuffixStructures.cabal +118/−0
- changelog.md +8/−0
- src/SAperformance.hs +56/−0
- src/mkesa.hs +107/−0
+ Data/SuffixStructure/ESA.hs view
@@ -0,0 +1,64 @@++-- | The suffix array data structure. Supports (de-) serialization via+-- aeson,cereal,binary.+--+-- Reading and writing to and from specialized "bio" formats is currently+-- open.+--+-- TODO compression during serialization?+-- TODO versioning?+-- TODO read sam/bam format?+-- TODO what about mmap for really large indices?++module Data.SuffixStructure.ESA where++import Data.Aeson+import Data.Binary+import Data.Default.Class+import Data.Int (Int8)+import Data.IntMap.Strict (IntMap)+import Data.Serialize+import Data.Vector.Binary+import Data.Vector.Cereal+import Data.Vector.Unboxed (Vector)+import GHC.Generics (Generic)+import qualified Data.IntMap.Strict as IM+import qualified Data.Vector.Unboxed as VU++++-- | The Suffix Array data type, together with the longest common prefix+-- table.+--+-- TODO skip table?+-- TODO inverse suffix array?+--+-- TODO maybe parametrize on the Int type (Int,Int64,Int32,Word's) This+-- will require better specialization of operations in @NaiveArray@ and+-- elsewhere. Otherwise performance drops quite noticable by @x5@ to @x10@.++data SA = SA+ { sa :: !(Vector Int) -- ^ the actual suffix array using 8byte Ints+ , lcp :: !(Vector Int8) -- ^ 1byte longest common prefix vector, negative number indicates to look at lcpLong+ , lcpLong :: !(IntMap Int) -- ^ lcp's that are unusual long, but this is sparse+ }+ deriving (Eq,Ord,Show,Generic)++instance Binary SA+instance FromJSON SA+instance Serialize SA+instance ToJSON SA++instance Default SA where+ def = SA VU.empty VU.empty IM.empty++-- | Automatically check 'lcp' and 'lcpLong' to return the real prefix+-- length in 'Int' (as opposed to 'Int8' storage of 'lcp').++lcpAt :: SA -> Int -> Int+lcpAt SA{..} k+ | p >= 0 = fromIntegral p+ | Just p' <- IM.lookup k lcpLong = p' -- by construction!+ where p = VU.unsafeIndex lcp k+{-# INLINE lcpAt #-}+
+ Data/SuffixStructure/NaiveArray.hs view
@@ -0,0 +1,80 @@++-- |+--+-- TODO need to check performance of 'drop' vs 'unsafeDrop'+--+-- TODO use 'Word' instead of 'Int' -- but check performance due to all+-- those conversions++module Data.SuffixStructure.NaiveArray where++import Data.ByteString (ByteString(..))+import Data.Int+import Data.IntMap.Strict (IntMap(..))+import Data.ListLike (ListLike)+import Data.Vector.Unboxed (Vector(..))+import Data.Word+import qualified Data.ByteString as B+import qualified Data.IntMap.Strict as IM+import qualified Data.ListLike as LL+import qualified Data.Vector.Algorithms.AmericanFlag as AA+import qualified Data.Vector.Algorithms.Intro as AI+import qualified Data.Vector.Unboxed as VU++import Data.SuffixStructure.ESA++++-- | Create Suffix Array via Introsort++genSA :: (ListLike ll a, Ord ll, Eq a) => ll -> SA+genSA ll = SA sa lcp lcpLong where+ sa = VU.modify (AI.sortBy srt) $ VU.enumFromN 0 (LL.length ll)+ (lcp,lcpLong) = buildLCP ll sa+ srt i j = LL.drop i ll `compare` LL.drop j ll+{-# INLINE genSA #-}++-- | Create Suffix Array via American Flag sort++genSAaf :: (ListLike ll a, Ord ll, AA.Lexicographic ll, Eq a) => ll -> SA+genSAaf ll = SA sa lcp lcpLong where+ sa = VU.modify (AA.sortBy srt strp bckt rdx) $ VU.enumFromN 0 (LL.length ll)+ (lcp,lcpLong) = buildLCP ll sa+ srt i j = LL.drop i ll `compare` LL.drop j ll+ strp _ i = i >= LL.length ll+ bckt = AA.size ll+ rdx i _ = AA.index i ll+{-# INLINE genSAaf #-}++-- | Build LCP array++buildLCP :: (ListLike ll a, Eq a) => ll -> VU.Vector Int -> (VU.Vector Int8, IM.IntMap Int)+buildLCP inp sa = (lcp,lcpLong) where+ lcp = (-1) `VU.cons` VU.zipWith golcp sa (VU.tail sa)+ lcpLong = VU.foldl' golcpLong IM.empty $ VU.zip4 (VU.enumFromN 1 $ VU.length sa) sa (VU.tail sa) (VU.tail lcp)+ golcp p k = let cpl = commonPrefixLength (LL.drop p inp) (LL.drop k inp)+ in if cpl <= 127+ then fromIntegral cpl+ else (-2)+ golcpLong im (k,s,t,l)+ | l== -2 = IM.insert k (commonPrefixLength (LL.drop s inp) (LL.drop t inp)) im+ | otherwise = im+{-# INLINE buildLCP #-}++-- | Return the shared prefix of two strings.++commonPrefix :: (ListLike ll a, Eq a) => ll -> ll -> ll+commonPrefix xs ys = LL.take k xs where+ k = commonPrefixLength xs ys+{-# INLINE commonPrefixLength #-}++-- | Return the length of the common prefix.++commonPrefixLength :: (ListLike ll a, Eq a) => ll -> ll -> Int+commonPrefixLength = go 0 where+ go !k !x !y+ | LL.null x || LL.null y = k+ | LL.head x == LL.head y = go (k+1) (LL.tail x) (LL.tail y)+ | otherwise = k+{-# INLINE commonPrefix #-}+
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright Christian Hoener zu Siederdissen 2011-2012++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++ * Neither the name of Christian Hoener zu Siederdissen nor the names of other+ contributors may be used to endorse or promote products derived+ from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ README.md view
@@ -0,0 +1,17 @@+# OrderedBits++[](https://travis-ci.org/choener/SuffixStructures)++A collection of suffix structures. While we currently only provide a naive+construction method, this method should still be acceptably fast in many cases+up to around 100M characters -- assuming the input sequence is well behaved in+some sense.++++#### Contact++Christian Hoener zu Siederdissen+choener@bioinf.uni-leipzig.de+Leipzig University, Leipzig, Germany+
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ SuffixStructures.cabal view
@@ -0,0 +1,118 @@+name: SuffixStructures+version: 0.0.1.0+author: Christian Hoener zu Siederdissen+copyright: Christian Hoener zu Siederdissen, 2014 - 2015+homepage: http://www.bioinf.uni-leipzig.de/~choener/+maintainer: choener@bioinf.uni-leipzig.de+category: Data, Data Structures+license: BSD3+license-file: LICENSE+build-type: Simple+stability: experimental+cabal-version: >= 1.10+tested-with: GHC == 7.8.4, GHC == 7.10.1+synopsis: Suffix array construction+description:+ Suffix array construction in Haskell. Currently, only a naive+ method is provided. More advanced construction methods might+ follow.++++Extra-Source-Files:+ README.md+ changelog.md++++flag llvm+ description: use llvm backend+ default: False+ manual: True++++library+ build-depends: base >= 4.7 && < 4.9+ , aeson == 0.8.*+ , binary == 0.7.*+ , bytestring == 0.10.*+ , cereal == 0.4.*+ , containers == 0.5.*+ , data-default-class == 0.0.1+ , ListLike >= 4.1.0.0 && < 4.3+ , primitive >= 0.5 && < 0.7+ , vector == 0.10.*+ , vector-algorithms == 0.6.*+ , vector-binary-instances == 0.2.*+ default-language:+ Haskell2010+ default-extensions: BangPatterns+ , DeriveGeneric+ , FlexibleContexts+ , PatternGuards+ , RankNTypes+ , RecordWildCards+ , TypeFamilies++ exposed-modules:+ Data.SuffixStructure.ESA+ Data.SuffixStructure.NaiveArray+ ghc-options:+ -O2 -funbox-strict-fields++++executable mkesa+ build-depends: base+ , cmdargs == 0.10.*+ , SuffixStructures+ , aeson+ , binary+ , containers+ , bytestring+ , cereal+ , vector+ default-language:+ Haskell2010+ hs-source-dirs:+ src+ main-is:+ mkesa.hs+ ghc-options:+ -O2 -funbox-strict-fields+ if flag(llvm)+ ghc-options:+ -fllvm -optlo-O3 -optlo-std-compile-opts++++benchmark BenchmarkSuffixStructures+ build-depends: base+ , bytestring+ , cmdargs == 0.10.*+ , criterion == 1.1.*+ , deepseq+ , mwc-random == 0.13.*+ , SuffixStructures+ , vector+ default-language:+ Haskell2010+ hs-source-dirs:+ src+ main-is:+ SAperformance.hs+ type:+ exitcode-stdio-1.0+ ghc-options:+ -O2 -funbox-strict-fields+ if flag(llvm)+ ghc-options:+ -fllvm -optlo-O3 -optlo-std-compile-opts++++source-repository head+ type: git+ location: git://github.com/choener/SuffixStructures+
+ changelog.md view
@@ -0,0 +1,8 @@+0.0.0.1+-------++- naive suffix arrays+- Criterion performance measurements+- QuickCheck+- travis-ci integration+
+ src/SAperformance.hs view
@@ -0,0 +1,56 @@++{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE RecordWildCards #-}+{-# LANGUAGE ScopedTypeVariables #-}++{-# OPTIONS_GHC -fno-cse #-}++module Main where++import Control.DeepSeq+import Criterion.Main+import Data.Char (ord,chr)+import Data.Vector.Unboxed (Vector(..))+import qualified Data.ByteString.Char8 as BS+import qualified Data.Vector.Unboxed as VU+import System.Console.CmdArgs+import System.Environment (withArgs)+import System.Random.MWC (withSystemRandom, asGenST, uniformVector)++import Data.SuffixStructure.ESA+import Data.SuffixStructure.NaiveArray++++data Options = Options+ { from :: Int+ , step :: Int+ , count :: Int+ , remaining :: [String]+ }+ deriving (Show,Data,Typeable)++options = Options+ { from = 1 &= help ""+ , step = 10 &= help ""+ , count = 5 &= help ""+ , remaining = def &= args+ }++main = do+ o@Options{..} <- cmdArgs options+ is :: [Vector Int] <- sequence . map (\k -> withSystemRandom . asGenST $ \gen -> uniformVector gen k) . take count $ (iterate (*step) from)+ let cs :: [Vector Char] = map (VU.map (chr . flip mod 256)) is+ let bs :: [BS.ByteString] = map (BS.pack . VU.toList) cs+ deepseq (is,cs,bs) $ putStrLn "vectors generated"+ withArgs remaining $ defaultMain+ [ bgroup "naive/introsort"+ $ zipWith (\k r -> bench (show k) $ whnf (VU.length . sa . genSA) r)+ (iterate (*step) from)+ cs+ , bgroup "naive/americanflag"+ $ zipWith (\k r -> bench (show k) $ whnf (VU.length . sa . genSAaf) r)+ (iterate (*step) from)+ bs+ ]+
+ src/mkesa.hs view
@@ -0,0 +1,107 @@++{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE RecordWildCards #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE TupleSections #-}++-- | Create ESA data structures from multiple sources and serialize to+-- disk.++module Main where++import Control.Applicative ((<$>))+import Control.Monad (forM_)+import qualified Data.Aeson as A+import Data.Function (on)+import Data.List (groupBy,sort)+import Data.Tuple (swap)+import qualified Data.Binary as DB+import qualified Data.ByteString as B+import qualified Data.ByteString.Lazy as BL+import qualified Data.IntMap.Strict as IM+import qualified Data.Serialize as DS+import qualified Data.Vector.Unboxed as VU+import System.Console.CmdArgs+import Text.Printf++import Data.SuffixStructure.ESA+import Data.SuffixStructure.NaiveArray++++data Type = Binary | Cereal | JSON+ deriving (Show,Data,Typeable)++data Options+ = RawBS+ { infile :: String+ , outfile :: String+ , outtype :: Type+ }+ | ReadTest+ { infile :: String+ , outfile :: String+ , intype :: Type+ }+ deriving (Show,Data,Typeable)++-- | Read a file as a simple raw bytestring and create the enhanced suffix+-- array. Note: the input is a strict bytestring. We need to load the+-- complete bytestring anyway, so we make it strict. We don't mmap because+-- the whole string is loaded.++oRawBS = RawBS+ { infile = def &= help ""+ , outfile = def &= help ""+ , outtype = Cereal &= help ""+ }++-- | Will read in the enhanced suffix array, print some statistics, and+-- quit.+--+-- TODO should actually become @Statistics@ I think.++oReadTest = ReadTest+ { intype = Cereal &= help ""+ }++main = do+ o <- cmdArgs $ modes [oRawBS, oReadTest]+ case o of+ RawBS{..} -> do i <- if null infile then B.getContents else B.readFile infile+ let !ar = genSA i+ let writer = if null outfile then BL.putStr else BL.writeFile outfile+ writer $ case outtype of+ Binary -> DB.encode ar+ Cereal -> DS.encodeLazy ar+ JSON -> A.encode ar+ ReadTest{..} -> do i <- if null infile then BL.getContents else BL.readFile infile+ let ar :: SA = case intype of+ Binary -> DB.decode i+ Cereal -> case DS.decodeLazy i of+ Right d -> d+ Left err -> error err+ JSON -> case A.decode i of+ Just ar -> ar+ Nothing -> error "not JSON input"+ printf "Suffix array size: %d\n" . VU.length $ sa ar+ let lcpdist = IM.fromListWith (+) $+ (VU.toList . VU.map (,1::Int) . VU.filter (>=0) . VU.map fromIntegral $ lcp ar) +++ (map (,1) . IM.elems $ lcpLong ar)+ printf "LCP array distribution\n"+ {-+ let tmax = maximum $ IM.keys lcpdist+ let ts = groupBy ((==) `on` (`div` 10)) [0 .. tmax]+ forM_ ts $ \tt -> do+ mapM_ (printf "%8d") $ tt+ printf "\n"+ mapM_ (printf "%8d") $ map (maybe 0 id . flip IM.lookup lcpdist) $ tt+ printf "\n\n"+ -}+ let hs = take 14 . reverse . sort . map swap $ IM.assocs lcpdist+ mapM_ (printf "%8d") $ map snd hs+ printf "\n"+ mapM_ (printf "%8d") $ map fst hs+ printf "\n\n"+