BCMtools 0.1.0 → 0.1.1
raw patch · 13 files changed
+372/−362 lines, 13 filesdep +cerealdep −binarydep −data-binary-ieee754dep ~basedep ~bytestringdep ~matrices
Dependencies added: cereal
Dependencies removed: binary, data-binary-ieee754
Dependency ranges changed: base, bytestring, matrices
Files
- BCMtools.cabal +13/−13
- exe/BCMtools/Convert.hs +36/−30
- exe/BCMtools/Types.hs +1/−1
- exe/BCMtools/View.hs +22/−40
- exe/Main.hs +7/−3
- src/BCM.hs +8/−8
- src/BCM/Binary.hs +180/−0
- src/BCM/DiskMatrix.hs +63/−104
- src/BCM/IOMatrix.hs +13/−12
- src/BCM/Matrix/Instances.hs +0/−112
- src/BCM/Visualize.hs +4/−4
- src/BCM/Visualize/Internal.hs +4/−5
- src/BCM/Visualize/Internal/Types.hs +21/−30
BCMtools.cabal view
@@ -1,7 +1,7 @@ name: BCMtools-version: 0.1.0+version: 0.1.1 synopsis: Big Contact Map Tools--- description: +description: Tools for working with Big Contact Map generated by HiC experiments license: MIT license-file: LICENSE author: Kai Zhang@@ -16,18 +16,18 @@ exposed-modules: BCM BCM.Visualize- BCM.Visualize.Internal- BCM.Visualize.Internal.Types BCM.DiskMatrix BCM.IOMatrix- BCM.Matrix.Instances+ BCM.Binary - -- other-modules: + other-modules: + BCM.Visualize.Internal+ BCM.Visualize.Internal.Types+ build-depends:- base >=4.0 && <5.0- , bytestring >=0.9- , data-binary-ieee754 >=0.4- , binary >=0.7+ base >=4.7 && <5.0+ , bytestring >=0.10.0+ , cereal , conduit , conduit-extra , colour@@ -51,15 +51,15 @@ BCMtools.View build-depends:- base >=4.0 && <5.0+ base >=4.7 && <5.0 , BCMtools- , binary+ , cereal , optparse-applicative , conduit , conduit-extra , resourcet , bytestring-lexing- , bytestring+ , bytestring >=0.10.0 , unordered-containers , split , data-default-class
exe/BCMtools/Convert.hs view
@@ -20,8 +20,7 @@ import BCMtools.Types import BCM (ContactMap, createContactMap, saveContactMap, closeContactMap)-import BCM.IOMatrix (DMatrix, MCSR, DSMatrix)-import BCM.Matrix.Instances ()+import BCM.IOMatrix (DMatrix, MCSR, DSMatrix, MSMatrix, MMatrix) convertOptions :: Parser Command convertOptions = fmap Convert $ ConvertOptions@@ -52,58 +51,65 @@ readInt' x = let (Just (i, left)) = B.readInt $ B.pack x in case () of _ | B.null left -> i- | left == "K" -> i * 1000- | left == "M" -> i * 1000000+ | left == "k" || left == "K" -> i * 1000 | otherwise -> i -convert :: FilePath -> FilePath -> ConvertOptions -> IO ()-convert input output opt = do+convert :: FilePath -> FilePath -> Bool -> ConvertOptions -> IO ()+convert input output onDisk opt = do genome <- case _genome opt of "hg19" -> return hg19 fl -> readGenome fl inputLength <- runResourceT $ Bin.sourceFile input $= Bin.lines $$ CL.fold (\i _ -> i+1) 0 line1 <- B.split '\t' <$> withFile input ReadMode B.hGetLine- let (fn, n) = case line1 of- [f1,f2,_] -> (field2 f1 f2, inputLength - 1)- [_,_,_,_,_] -> (field5, inputLength)+ let fn = case line1 of+ [_,_,_] -> field2+ [_,_,_,_,_] -> field5 _ -> error "Please check your input format"- runner x = runResourceT $+ runner f = runResourceT $ Bin.sourceFile input $=- Bin.lines $= fn $$- createContactMap output rows cols (_resolution opt) (Just x)+ Bin.lines $= f $$+ createContactMap output rows cols (_resolution opt) (Just inputLength) rows = getChrSize genome $ _rownames opt cols = getChrSize genome $ _colnames opt case () of _ | _sparse opt && _symmetric opt -> do- cm <- runner n :: IO (ContactMap MCSR)- saveContactMap cm- closeContactMap cm- | _sparse opt -> do- cm <- runner n :: IO (ContactMap MCSR)- saveContactMap cm- closeContactMap cm- | _symmetric opt -> do- cm <- runner n :: IO (ContactMap DSMatrix)- saveContactMap cm- closeContactMap cm+ if onDisk+ then createMapWith (runner fn :: IO (ContactMap MCSR))+ else createMapWith (runner fn :: IO (ContactMap MCSR))+ | _sparse opt ->+ if onDisk+ then createMapWith (runner fn :: IO (ContactMap MCSR))+ else createMapWith (runner fn :: IO (ContactMap MCSR))+ | _symmetric opt ->+ if onDisk+ then createMapWith (runner fn :: IO (ContactMap DSMatrix))+ else createMapWith (runner fn :: IO (ContactMap MSMatrix)) | otherwise -> do- cm <- runner n :: IO (ContactMap DMatrix)- saveContactMap cm- closeContactMap cm+ if onDisk+ then createMapWith (runner fn :: IO (ContactMap DMatrix))+ else createMapWith (runner fn :: IO (ContactMap MMatrix)) where+ createMapWith run = do cm <- run+ saveContactMap cm+ closeContactMap cm readGenome x = do c <- B.readFile x return $ M.fromList $ map ((\[a,b] -> (B.unpack a, readInt b)) . B.words) $ B.lines c- getChrSize g = map (B.pack &&& flip (M.lookupDefault errMsg) g)- where errMsg = error "Unknown chromosome"- field2 chr1 chr2 = do+ getChrSize g = map lookup'+ where+ lookup' x = (B.pack x, M.lookupDefault errMsg x g)+ where+ errMsg = error $ "Unknown chromosome: " ++ x+ field2 = do _ <- await CL.map f where f l = let [x1,x2,v] = B.split '\t' l- in (chr1, readInt x1, chr2, readInt x2, readDouble' v)+ in (B.pack chr1, readInt x1, B.pack chr2, readInt x2, readDouble' v)+ [chr1] = _rownames opt+ [chr2] = _colnames opt field5 = CL.map f where f l = let [x1,x2,x3,x4,x5] = B.split '\t' l
exe/BCMtools/Types.hs view
@@ -9,6 +9,7 @@ { _input :: FilePath , _output :: FilePath , _command :: Command+ , _onDisk :: Bool } data Command = Convert ConvertOptions@@ -25,5 +26,4 @@ data ViewOptions = ViewOptions { _valueRange :: (Double, Double)- , _inMem :: Bool }
exe/BCMtools/View.hs view
@@ -4,8 +4,8 @@ , viewOptions ) where -import qualified Data.ByteString.Lazy as L-import Data.Binary.Get (runGet, getWord32le)+import qualified Data.ByteString as B+import Data.Serialize (runGet, getWord32le) import System.IO import Options.Applicative import Data.Conduit@@ -17,8 +17,8 @@ import BCMtools.Types import BCM (ContactMap, _matrix, closeContactMap, openContactMap) import BCM.IOMatrix (DMatrix, MCSR, DSMatrix, MMatrix, MSMatrix)-import BCM.Matrix.Instances import BCM.Visualize+import BCM.Binary (ds_matrix_magic, d_matrix_magic, sp_matrix_magic) viewOptions :: Parser Command viewOptions = fmap View $ ViewOptions@@ -26,59 +26,41 @@ ( long "range" <> short 'r' <> metavar "Heatmap range" ))- <*> switch- ( long "memory"- <> help "store matrix in memory" ) where f x = let a = read $ takeWhile (/='-') x b = read $ tail $ dropWhile (/='-') x in (a, b) -view :: FilePath -> FilePath -> ViewOptions -> IO ()-view input output opt = do- magic <- withFile input ReadMode $ \h -> do+view :: FilePath -> FilePath -> Bool -> ViewOptions -> IO ()+view input output onDisk opt = do+ Right magic <- withFile input ReadMode $ \h -> do hSeek h AbsoluteSeek 4- p <- fromIntegral . runGet getWord32le <$> L.hGet h 4+ Right p <- fmap fromIntegral . runGet getWord32le <$> B.hGet h 4 hSeek h AbsoluteSeek p- runGet getWord32le <$> L.hGet h 4+ runGet getWord32le <$> B.hGet h 4 case () of _ | magic == d_matrix_magic -> do hPutStrLn stderr "Format: Dense matrix"- if _inMem opt- then do- cm <- openContactMap input :: IO (ContactMap MMatrix)- draw cm- closeContactMap cm- else do- cm <- openContactMap input :: IO (ContactMap DMatrix)- draw cm- closeContactMap cm+ if onDisk+ then openWith (openContactMap input :: IO (ContactMap DMatrix))+ else openWith (openContactMap input :: IO (ContactMap MMatrix)) | magic == ds_matrix_magic -> do hPutStrLn stderr "Format: Dense symmetric matrix"- if _inMem opt- then do- cm <- openContactMap input :: IO (ContactMap MSMatrix)- draw cm- closeContactMap cm- else do- cm <- openContactMap input :: IO (ContactMap DSMatrix)- draw cm- closeContactMap cm- | magic == sp_matrix_magic -> do+ if onDisk+ then openWith (openContactMap input :: IO (ContactMap DSMatrix))+ else openWith (openContactMap input :: IO (ContactMap MSMatrix))+ _ | magic == sp_matrix_magic -> do hPutStrLn stderr "Format: Sparse matrix"- if _inMem opt- then do- cm <- openContactMap input :: IO (ContactMap MCSR)- draw cm- closeContactMap cm- else do- cm <- openContactMap input :: IO (ContactMap MCSR)- draw cm- closeContactMap cm+ if onDisk+ then openWith (openContactMap input :: IO (ContactMap MCSR))+ else openWith (openContactMap input :: IO (ContactMap MCSR)) | otherwise -> error "unknown format" where- draw x = runResourceT $ drawMatrix (_matrix x) drawopt $= CL.map L.toStrict $$ Bin.sinkFile output+ openWith fn = do cm <- fn+ draw cm+ closeContactMap cm+ draw x = runResourceT $ drawMatrix (_matrix x) drawopt $$ Bin.sinkFile output drawopt = def { _range = _valueRange opt } {-# INLINE view #-}
exe/Main.hs view
@@ -16,6 +16,10 @@ <> short 'o' <> metavar "OUTPUT" ) <*> cmd+ <*> switch+ ( long "disk"+ <> help "store matrix on disk" )+ bcmtoolsOptions :: Parser BCMtoolsOptions bcmtoolsOptions = subparser $@@ -29,9 +33,9 @@ ) runBCMtools :: BCMtoolsOptions -> IO () -runBCMtools (BCMtoolsOptions input output bcmopt) = case bcmopt of- Convert opt -> convert input output opt- View opt -> view input output opt+runBCMtools (BCMtoolsOptions input output bcmopt onDisk) = case bcmopt of+ Convert opt -> convert input output onDisk opt+ View opt -> view input output onDisk opt main :: IO () main = execParser opts >>= runBCMtools
src/BCM.hs view
@@ -12,11 +12,10 @@ import Control.Applicative ((<$>)) import Control.Monad (when, guard) import Control.Monad.IO.Class (MonadIO(..))-import qualified Data.ByteString.Lazy as L+import qualified Data.ByteString as L import qualified Data.ByteString.Char8 as B import qualified Data.HashMap.Strict as M-import Data.Binary.Put-import Data.Binary.Get+import Data.Serialize import Data.Conduit import qualified Data.Conduit.List as CL import System.IO@@ -93,22 +92,23 @@ rows = encodeLab . M.toList $ rowChr cols = encodeLab . M.toList $ colChr encodeLab xs = L.concat $ concatMap (\(chr, (a,b)) ->- [L.fromStrict chr, "\0", DM.toByteString a, DM.toByteString b]) xs+ [chr, "\0", runPut $ putWord64le $ fromIntegral a, runPut $ putWord64le $ fromIntegral b]) xs offset = fromIntegral $ 4 + 4 + 4 + L.length rowAndcol + 1 openContactMap :: (IOM.IOMatrix m t a, MonadIO io, mat ~ m t a) => FilePath -> io (ContactMap mat) openContactMap fl = liftIO $ do h <- openFile fl ReadWriteMode - magic <- runGet getWord32le <$> L.hGet h 4+ Right magic <- runGet getWord32le <$> L.hGet h 4 guard $ magic == contact_map_magic _ <- runGet getWord32le <$> L.hGet h 4- res <- fromIntegral . runGet getWord32le <$> L.hGet h 4+ Right res <- fmap fromIntegral . runGet getWord32le <$> L.hGet h 4 rows <- M.fromList <$> getChrs [] h cols <- M.fromList <$> getChrs [] h _ <- runGet getWord8 <$> L.hGet h 1 mat <- IOM.hReadMatrix h+ x <- IOM.unsafeIndexM mat (1,1) return $ ContactMap rows cols res mat h where getChrs acc h = do@@ -116,8 +116,8 @@ if B.null chr then return acc else do- a <- fromIntegral . runGet getWord64le <$> L.hGet h 8- b <- fromIntegral . runGet getWord64le <$> L.hGet h 8+ Right a <- fmap fromIntegral . runGet getWord64le <$> L.hGet h 8+ Right b <- fmap fromIntegral . runGet getWord64le <$> L.hGet h 8 getChrs ((chr, (a, b)) : acc) h getByteStringNul h = B.concat <$> go [] where
+ src/BCM/Binary.hs view
@@ -0,0 +1,180 @@+{-# LANGUAGE FlexibleInstances #-}+module BCM.Binary where++import Control.Applicative ((<$>))+import Foreign+import System.IO+import System.IO.Unsafe (unsafePerformIO)++import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Storable as S++import qualified Data.Matrix.Dense.Generic as DM+import qualified Data.Matrix.Symmetric as DS+import qualified Data.Matrix.Sparse.Generic as SM++class Swappable a where+ byteSwap :: a -> a++instance Swappable Double where+ byteSwap x = unsafePerformIO $ with x $ \ptr -> do+ word64 <- byteSwap64 <$> peek (castPtr ptr)+ poke (castPtr ptr) word64+ peek ptr + {-# INLINE byteSwap #-}++instance Swappable Int64 where+ byteSwap = fromIntegral . byteSwap64 . fromIntegral+ {-# INLINE byteSwap #-}++class BinaryData a where+ hGetData :: Handle -> Bool -> IO a+ hPutData :: Handle -> a -> IO ()+ size :: a -> Int++instance BinaryData Int64 where+ hGetData hdl = fmap fromIntegral . (hGetData hdl :: Bool -> IO Word64)+ hPutData hdl = (hPutData hdl :: Word64 -> IO ()). fromIntegral+ size _ = 8+ {-# INLINE hGetData #-}+ {-# INLINE hPutData #-}++instance BinaryData Double where+ hGetData hdl byteSwapped = alloca $ \ptr -> do+ hGetBuf hdl ptr 8+ word64 <- if byteSwapped+ then byteSwap64 <$> peek ptr+ else peek ptr+ poke ptr word64+ peek $ castPtr ptr+ hPutData hdl x = with x $ \ptr -> hPutBuf hdl ptr 8+ size _ = 8+ {-# INLINE hGetData #-}+ {-# INLINE hPutData #-}++instance BinaryData Word64 where+ hGetData hdl byteSwapped = alloca $ \ptr -> do+ hGetBuf hdl ptr 8+ if byteSwapped+ then byteSwap64 <$> peek ptr+ else peek ptr+ hPutData hdl x = with x $ \ptr -> hPutBuf hdl ptr 8+ size _ = 8+ {-# INLINE hGetData #-}+ {-# INLINE hPutData #-}++instance BinaryData Word32 where+ hGetData hdl byteSwapped = alloca $ \ptr -> do+ hGetBuf hdl ptr 4+ if byteSwapped+ then byteSwap32 <$> peek ptr+ else peek ptr+ hPutData hdl x = with x $ \ptr -> hPutBuf hdl ptr 4+ size _ = 4+ {-# INLINE hGetData #-}+ {-# INLINE hPutData #-}++d_matrix_magic :: Word32+d_matrix_magic = 0x22D20B77++instance BinaryData (DM.Matrix U.Vector Double) where+ hPutData hdl (DM.Matrix r c _ _ vec) = do+ hPutData hdl d_matrix_magic+ hPutData hdl (fromIntegral r :: Int64)+ hPutData hdl (fromIntegral c :: Int64)+ G.mapM_ (hPutData hdl) vec++ hGetData hdl _ = do+ magic <- hGetData hdl False+ let byteSwapped | magic == d_matrix_magic = False+ | byteSwap32 magic == d_matrix_magic = True+ | otherwise = error "Read matrix fail: wrong signature"+ r <- fromIntegral <$> (hGetData hdl byteSwapped :: IO Int64)+ c <- fromIntegral <$> (hGetData hdl byteSwapped :: IO Int64)+ vec <- hGetDoubleVector hdl (r*c)+ if byteSwapped+ then return $ DM.Matrix r c c 0 $ G.map byteSwap $ G.convert vec+ else return $ DM.Matrix r c c 0 $ G.convert vec+ size = undefined+ {-# INLINE hGetData #-}+ {-# INLINE hPutData #-}++ds_matrix_magic :: Word32+ds_matrix_magic = 0x33D31A66++instance BinaryData (DS.SymMatrix U.Vector Double) where+ hPutData hdl (DS.SymMatrix n vec) = do+ hPutData hdl ds_matrix_magic+ hPutData hdl (fromIntegral n :: Int64)+ G.mapM_ (hPutData hdl) vec++ hGetData hdl _ = do+ magic <- hGetData hdl False+ let byteSwapped | magic == ds_matrix_magic = False+ | byteSwap32 magic == ds_matrix_magic = True+ | otherwise = error "Read matrix fail: wrong signature"+ n <- fromIntegral <$> (hGetData hdl byteSwapped :: IO Int64)+ let len = ((n+1)*n) `shiftR` 1+ vec <- hGetDoubleVector hdl len+ if byteSwapped+ then return $ DS.SymMatrix n $ G.map byteSwap $ G.convert vec+ else return $ DS.SymMatrix n $ G.convert vec+ size = undefined+ {-# INLINE hGetData #-}+ {-# INLINE hPutData #-}++sp_matrix_magic :: Word32+sp_matrix_magic = 0x177BFFA0++instance BinaryData (SM.CSR U.Vector Double) where+ hPutData hdl (SM.CSR r c vec ci rp) = do+ hPutData hdl sp_matrix_magic+ hPutData hdl (fromIntegral r :: Int64)+ hPutData hdl (fromIntegral c :: Int64)+ hPutData hdl (fromIntegral n :: Int64)+ G.mapM_ (hPutData hdl) vec+ G.mapM_ (hPutData hdl . (fromIntegral :: Int -> Int64)) ci+ G.mapM_ (hPutData hdl . (fromIntegral :: Int -> Int64)) rp+ where+ n = G.length vec++ hGetData hdl _ = do+ magic <- hGetData hdl False+ let byteSwapped | magic == sp_matrix_magic = False+ | byteSwap32 magic == sp_matrix_magic = True+ | otherwise = error "Read matrix fail: wrong signature"+ r <- fromIntegral <$> (hGetData hdl byteSwapped :: IO Int64)+ c <- fromIntegral <$> (hGetData hdl byteSwapped :: IO Int64)+ n <- fromIntegral <$> (hGetData hdl byteSwapped :: IO Int64)+ vec <- G.convert <$> hGetDoubleVector hdl n+ ci <- G.convert <$> hGetIntVector hdl n+ rp <- G.convert <$> hGetIntVector hdl c+ if byteSwapped+ then return $ SM.CSR r c (G.map byteSwap vec)+ (G.map (fromIntegral . (byteSwap :: Int64 -> Int64) . fromIntegral) ci)+ (G.map (fromIntegral . (byteSwap :: Int64 -> Int64) . fromIntegral) rp)+ else return $ SM.CSR r c vec ci rp+ {-# INLINE hGetData #-}+ {-# INLINE hPutData #-}++ size (SM.CSR _ _ vec ci rp) = 28 + 8 * (G.length vec + G.length ci + G.length rp)+++------------------------------------------------------------------------------++hGetIntVector :: Handle -> Int -> IO (S.Vector Int)+hGetIntVector hdl n = do+ ptr <- mallocBytes (n*8)+ hGetBuf hdl ptr (n*8)+ ptr' <- newForeignPtr_ ptr+ return $ S.unsafeFromForeignPtr0 ptr' n+{-# INLINE hGetIntVector #-}++hGetDoubleVector :: Handle -> Int -> IO (S.Vector Double)+hGetDoubleVector hdl n = do+ ptr <- mallocBytes (n*8)+ hGetBuf hdl ptr (n*8)+ ptr' <- newForeignPtr_ ptr+ return $ S.unsafeFromForeignPtr0 ptr' n+{-# INLINE hGetDoubleVector #-}
src/BCM/DiskMatrix.hs view
@@ -5,7 +5,6 @@ {-# LANGUAGE ScopedTypeVariables #-} module BCM.DiskMatrix ( Offset- , DiskData(..) , DiskMatrix(..) , read , write@@ -17,67 +16,20 @@ import Control.Monad (replicateM_) import Control.Monad.IO.Class (MonadIO(..)) import Control.Applicative ((<$>))-import qualified Data.ByteString.Lazy as L+import qualified Data.ByteString as B import Data.Bits (shiftR)-import Data.Binary.IEEE754 (getFloat64le, putFloat64le)-import Data.Binary.Put (putWord32le, putWord64le, runPut)-import Data.Binary.Get (getWord32le, getWord64le, runGet)+import Data.Int+import Data.Word import qualified Data.Vector.Generic as G import System.IO-import Numeric (showHex)--import BCM.Matrix.Instances (d_matrix_magic, ds_matrix_magic)+import BCM.Binary type Offset = Integer -class DiskData a where- -- | The default value- zero :: a-- sizeOf :: a -> Int-- fromByteString :: L.ByteString -> a- - toByteString :: a -> L.ByteString-- hRead1 :: MonadIO m => Handle -> m a- hRead1 h = liftIO $ fmap fromByteString $ L.hGet h $ sizeOf (undefined :: a)- {-# INLINE hRead1 #-}-- hWrite1 :: MonadIO m => Handle -> a -> m ()- hWrite1 h = liftIO . L.hPut h . toByteString- {-# INLINE hWrite1 #-}-- {-# MINIMAL zero, sizeOf, fromByteString, toByteString #-}--instance DiskData Double where- zero = 0.0-- sizeOf _ = 8- {-# INLINE sizeOf #-}-- fromByteString = runGet getFloat64le- {-# INLINE fromByteString #-}-- toByteString = runPut . putFloat64le- {-# INLINE toByteString #-}--instance DiskData Int where- zero = 0-- sizeOf _ = 8- {-# INLINE sizeOf #-}-- fromByteString = fromIntegral . runGet getWord64le- {-# INLINE fromByteString #-}-- toByteString = runPut . putWord64le . fromIntegral- {-# INLINE toByteString #-}- -- | Matrix stored in binary file-class DiskData a => DiskMatrix m a where+class BinaryData a => DiskMatrix m a where elemType :: m a -> a- elemType _ = undefined+ elemType = undefined {-# INLINE elemType #-} hReadMatrixEither :: MonadIO io => Handle -> io (Either String (m a))@@ -107,94 +59,100 @@ (r,c) = dim mat {-# INLINE write #-} -data DMatrix a = DMatrix !Int -- ^ rows- !Int -- ^ cols+data DMatrix a = DMatrix !Int -- rows+ !Int -- cols !Offset -- offset- !Handle -- ^ file handle+ !Bool -- byte swap+ !Handle -- file handle -instance DiskData a => DiskMatrix DMatrix a where+instance (Swappable a, BinaryData a) => DiskMatrix DMatrix a where hReadMatrixEither h = liftIO $ do p <- hTell h- magic <- runGet getWord32le <$> L.hGet h 4- if magic == d_matrix_magic- then do- r <- hRead1 h- c <- hRead1 h- return $ Right $ DMatrix r c (p+20) h- else return $ Left $ "Read fail: wrong signature: 0x" ++ showHex magic ""+ magic <- hGetData h False+ let byteSwapped | magic == d_matrix_magic = False+ | byteSwap32 magic == d_matrix_magic = True+ | otherwise = error "Read matrix fail: wrong signature"+ r <- fromIntegral <$> (hGetData h byteSwapped :: IO Int64)+ c <- fromIntegral <$> (hGetData h byteSwapped :: IO Int64)+ return $ Right $ DMatrix r c (p+20) byteSwapped h {-# INLINE hReadMatrixEither #-} - dim (DMatrix r c _ _) = (r,c)+ dim (DMatrix r c _ _ _) = (r,c) {-# INLINE dim #-} replicate h (r,c) x = liftIO $ do p <- hTell h- L.hPut h $ runPut $ putWord32le d_matrix_magic- hWrite1 h r- hWrite1 h c- replicateM_ (r*c) $ hWrite1 h x- return $ DMatrix r c (p+20) h+ hPutData h d_matrix_magic+ hPutData h (fromIntegral r :: Word64)+ hPutData h (fromIntegral c :: Word64)+ replicateM_ (r*c) $ hPutData h x+ return $ DMatrix r c (p+20) False h {-# INLINE replicate #-} - unsafeRead mat@(DMatrix _ c offset h) (i,j) = liftIO $ do- hSeek h AbsoluteSeek $ offset + fromIntegral (sizeOf (elemType mat) * idx c i j)- hRead1 h+ unsafeRead mat@(DMatrix _ c offset byteSwapped h) (i,j) = liftIO $ do+ hSeek h AbsoluteSeek $ offset + fromIntegral (size (elemType mat) * idx c i j)+ hGetData h byteSwapped {-# INLINE unsafeRead #-} - unsafeWrite mat@(DMatrix _ c offset h) (i,j) x = liftIO $ do- hSeek h AbsoluteSeek $ offset + fromIntegral (sizeOf (elemType mat) * idx c i j)- hWrite1 h x+ unsafeWrite mat@(DMatrix _ c offset byteSwapped h) (i,j) x = liftIO $ do+ hSeek h AbsoluteSeek $ offset + fromIntegral (size (elemType mat) * idx c i j)+ let x' | byteSwapped = byteSwap x+ | otherwise = x+ hPutData h x' {-# INLINE unsafeWrite #-} - unsafeReadRow mat@(DMatrix _ c offset h) i = liftIO $ do- hSeek h AbsoluteSeek $ offset + fromIntegral (size * c * i)- G.replicateM c $ hRead1 h+ unsafeReadRow mat@(DMatrix _ c offset byteSwapped h) i = liftIO $ do+ hSeek h AbsoluteSeek $ offset + fromIntegral (s * c * i)+ G.replicateM c $ hGetData h byteSwapped where- size = sizeOf $ elemType mat+ s = size $ elemType mat {-# INLINE unsafeReadRow #-} - close (DMatrix _ _ _ h) = liftIO $ hClose h+ close (DMatrix _ _ _ _ h) = liftIO $ hClose h -- | Symmetric matrix-data DSMatrix a = DSMatrix !Int -- ^ size- !Offset -- ^ offset- !Handle -- ^ file handle+data DSMatrix a = DSMatrix !Int -- size+ !Offset -- offset+ !Bool -- byteSwapped+ !Handle -- file handle -instance DiskData a => DiskMatrix DSMatrix a where+instance (Swappable a, BinaryData a) => DiskMatrix DSMatrix a where hReadMatrixEither h = liftIO $ do p <- hTell h- magic <- runGet getWord32le <$> L.hGet h 4- if magic == ds_matrix_magic- then do- n <- hRead1 h- return $ Right $ DSMatrix n (p+12) h- else return $ Left $ "Read fail: wrong signature: 0x" ++ showHex magic ""+ magic <- hGetData h False+ let byteSwapped | magic == ds_matrix_magic = False+ | byteSwap32 magic == ds_matrix_magic = True+ | otherwise = error "Read matrix fail: wrong signature"+ n <- fromIntegral <$> (hGetData h byteSwapped :: IO Int64)+ return $ Right $ DSMatrix n (p+12) byteSwapped h {-# INLINE hReadMatrixEither #-} - dim (DSMatrix n _ _) = (n,n)+ dim (DSMatrix n _ _ _) = (n,n) {-# INLINE dim #-} replicate h (r,c) x | r /= c = error "Not a sqaure matrix" | otherwise = liftIO $ do p <- hTell h- L.hPut h $ runPut $ putWord32le ds_matrix_magic- hWrite1 h r- replicateM_ (((r+1)*r) `shiftR` 1) $ hWrite1 h x- return $ DSMatrix r (p+12) h+ hPutData h ds_matrix_magic+ hPutData h (fromIntegral r :: Int64)+ replicateM_ (((r+1)*r) `shiftR` 1) $ hPutData h x+ return $ DSMatrix r (p+12) False h {-# INLINE replicate #-} - unsafeRead mat@(DSMatrix n offset h) (i,j) = liftIO $ do- hSeek h AbsoluteSeek $ offset + fromIntegral (sizeOf (elemType mat) * idx' n i j)- hRead1 h+ unsafeRead mat@(DSMatrix n offset byteSwapped h) (i,j) = liftIO $ do+ hSeek h AbsoluteSeek $ offset + fromIntegral (size (elemType mat) * idx' n i j)+ hGetData h byteSwapped {-# INLINE unsafeRead #-} - unsafeWrite mat@(DSMatrix n offset h) (i,j) x = liftIO $ do- hSeek h AbsoluteSeek $ offset + fromIntegral (sizeOf (elemType mat) * idx' n i j)- hWrite1 h x+ unsafeWrite mat@(DSMatrix n offset byteSwapped h) (i,j) x = liftIO $ do+ hSeek h AbsoluteSeek $ offset + fromIntegral (size (elemType mat) * idx' n i j)+ let x' | byteSwapped = byteSwap x+ | otherwise = x+ hPutData h x' {-# INLINE unsafeWrite #-} - close (DSMatrix _ _ h) = liftIO $ hClose h+ close (DSMatrix _ _ _ h) = liftIO $ hClose h ------------------------------------------------------------------------------ -- helper functions@@ -210,3 +168,4 @@ idx' n i j | i <= j = (i * (2 * n - i - 1)) `shiftR` 1 + j | otherwise = (j * (2 * n - j - 1)) `shiftR` 1 + i {-# INLINE idx' #-}+
src/BCM/IOMatrix.hs view
@@ -14,8 +14,8 @@ import Control.Monad (when, forM_) import Control.Applicative ((<$>)) import Control.Monad.IO.Class (MonadIO(..))+import qualified Data.ByteString as B import qualified Data.ByteString.Lazy as L-import Data.Binary (Binary, encode, decode) import qualified Data.Matrix.Generic as MG import qualified Data.Matrix.Generic.Mutable as MGM@@ -31,6 +31,7 @@ import Text.Printf (printf) import System.IO +import BCM.Binary import qualified BCM.DiskMatrix as DM type MMatrix = IOMat (Matrix U.Vector) Double@@ -64,7 +65,7 @@ -- | Just a wrapper newtype IODMat m a = IODMat { unwrapD :: m a } -instance (DM.DiskData a, DM.DiskMatrix m a) => IOMatrix IODMat m a where+instance (Zero a, BinaryData a, DM.DiskMatrix m a) => IOMatrix IODMat m a where dim = DM.dim . unwrapD {-# INLINE dim #-} @@ -82,7 +83,7 @@ {-# INLINE hReadMatrix #-} hCreateMatrix handle (r,c) _ = do- mat <- DM.replicate handle (r,c) DM.zero+ mat <- DM.replicate handle (r,c) zero CL.mapM_ $ \((i,j), x) -> DM.unsafeWrite mat (i,j) x return $ IODMat mat {-# INLINE hCreateMatrix #-}@@ -91,8 +92,9 @@ newtype IOMat m a = IOMat { unwrap :: m a } instance ( U.Unbox a- , DM.DiskData a- , Binary (m U.Vector a)+ , Zero a+ , BinaryData a+ , BinaryData (m U.Vector a) , MG.Matrix m U.Vector a ) => IOMatrix IOMat (m U.Vector) a where dim = MG.dim . unwrap@@ -105,15 +107,15 @@ {-# INLINE unsafeTakeRowM #-} hReadMatrix handle = liftIO $ do- mat <- decode <$> L.hGetContents handle+ mat <- hGetData handle False return $ IOMat mat {-# INLINE hReadMatrix #-} - hSaveMatrix handle (IOMat mat) = liftIO . L.hPutStr handle . encode $ mat+ hSaveMatrix handle (IOMat mat) = liftIO . hPutData handle $ mat {-# INLINE hSaveMatrix #-} hCreateMatrix _ (r,c) _ = do- mat <- liftIO $ MGM.replicate (r,c) DM.zero+ mat <- liftIO $ MGM.replicate (r,c) zero CL.mapM_ $ \((i,j), x) -> liftIO $ MGM.unsafeWrite mat (i,j) x mat' <- liftIO $ MG.unsafeFreeze mat return $ IOMat mat'@@ -124,8 +126,7 @@ instance ( Zero a , U.Unbox a- , DM.DiskData a- , Binary (CSR U.Vector a)+ , BinaryData (CSR U.Vector a) ) => IOMatrix IOSMat (CSR U.Vector) a where dim = MG.dim . unwrapS {-# INLINE dim #-}@@ -137,11 +138,11 @@ {-# INLINE unsafeTakeRowM #-} hReadMatrix handle = liftIO $ do- mat <- decode <$> L.hGetContents handle+ mat <- hGetData handle undefined return $ IOSMat mat {-# INLINE hReadMatrix #-} - hSaveMatrix handle (IOSMat mat) = liftIO . L.hPutStr handle . encode $ mat+ hSaveMatrix handle (IOSMat mat) = liftIO . hPutData handle $ mat {-# INLINE hSaveMatrix #-} hCreateMatrix _ (r,c) (Just n) = do
− src/BCM/Matrix/Instances.hs
@@ -1,112 +0,0 @@-{-# LANGUAGE FlexibleInstances #-}-{-# LANGUAGE FlexibleContexts #-}-module BCM.Matrix.Instances where--import Data.Bits (shiftR)-import Control.Applicative ((<$>))-import Data.Binary (Binary(..), Get, Put, Word32)-import Data.Binary.IEEE754 (getFloat64le, putFloat64le)-import Data.Binary.Get (getWord64le, getWord32le)-import Data.Binary.Put (putWord64le, putWord32le)-import Control.Monad (guard)-import qualified Data.Vector.Unboxed as U-import qualified Data.Vector.Generic as G--import qualified Data.Matrix.Dense.Generic as DM-import qualified Data.Matrix.Symmetric as DS-import qualified Data.Matrix.Sparse.Generic as SM--d_matrix_magic :: Word32-d_matrix_magic = 0x22D20B77--instance Binary (DM.Matrix U.Vector Double) where- put = putMatrix putFloat64le- get = getMatrix getFloat64le--instance Binary (DM.Matrix U.Vector Int) where- put = putMatrix $ putWord64le . fromIntegral- get = getMatrix $ fromIntegral <$> getWord64le--putMatrix :: G.Vector v a => (a -> Put) -> DM.Matrix v a -> Put-putMatrix putElement (DM.Matrix r c _ _ vec) = do- putWord32le d_matrix_magic- putWord64le . fromIntegral $ r- putWord64le . fromIntegral $ c- G.mapM_ putElement vec-{-# INLINE putMatrix #-}--getMatrix :: G.Vector v a => Get a -> Get (DM.Matrix v a)-getMatrix getElement = do- m <- getWord32le- guard $ m == d_matrix_magic- r <- fromIntegral <$> getWord64le- c <- fromIntegral <$> getWord64le- vec <- G.replicateM (r*c) getElement- return $ DM.Matrix r c c 0 vec-{-# INLINE getMatrix #-}---ds_matrix_magic :: Word32-ds_matrix_magic = 0x33D31A66--instance Binary (DS.SymMatrix U.Vector Double) where- put = putSymMatrix putFloat64le- get = getSymMatrix getFloat64le--instance Binary (DS.SymMatrix U.Vector Int) where- put = putSymMatrix $ putWord64le . fromIntegral- get = getSymMatrix $ fromIntegral <$> getWord64le--putSymMatrix :: G.Vector v a => (a -> Put) -> DS.SymMatrix v a -> Put-putSymMatrix putElement (DS.SymMatrix n vec) = do- putWord32le ds_matrix_magic- putWord64le . fromIntegral $ n- G.mapM_ putElement vec-{-# INLINE putSymMatrix #-}--getSymMatrix :: G.Vector v a => Get a -> Get (DS.SymMatrix v a)-getSymMatrix getElement = do- m <- getWord32le- guard $ m == ds_matrix_magic- n <- fromIntegral <$> getWord64le- vec <- G.replicateM (((n+1)*n) `shiftR` 1) getElement- return $ DS.SymMatrix n vec-{-# INLINE getSymMatrix #-}---sp_matrix_magic :: Word32-sp_matrix_magic = 0x177BFFA0--instance Binary (SM.CSR U.Vector Double) where- put = putCSR putFloat64le- get = getCSR getFloat64le--instance Binary (SM.CSR U.Vector Int) where- put = putCSR $ putWord64le . fromIntegral- get = getCSR $ fromIntegral <$> getWord64le--putCSR :: G.Vector v a => (a -> Put) -> SM.CSR v a -> Put-putCSR putElement (SM.CSR r c vec ci rp) = do- putWord32le sp_matrix_magic- putWord64le . fromIntegral $ r- putWord64le . fromIntegral $ c- putWord64le . fromIntegral $ n- G.mapM_ putElement vec- G.mapM_ (putWord64le . fromIntegral) ci- G.mapM_ (putWord64le . fromIntegral) rp- where- n = G.length vec-{-# INLINE putCSR #-}--getCSR :: G.Vector v a => Get a -> Get (SM.CSR v a)-getCSR getElement = do- m <- getWord32le- guard $ m == sp_matrix_magic- r <- fromIntegral <$> getWord64le- c <- fromIntegral <$> getWord64le- n <- fromIntegral <$> getWord64le- vec <- G.replicateM n getElement- ci <- G.replicateM n $ fromIntegral <$> getWord64le- rp <- G.replicateM c $ fromIntegral <$> getWord64le- return $ SM.CSR r c vec ci rp-{-# INLINE getCSR #-}
src/BCM/Visualize.hs view
@@ -7,8 +7,8 @@ ) where import Control.Monad.IO.Class-import qualified Data.ByteString.Lazy as L-import Data.Binary (encode)+import qualified Data.ByteString as B+import Data.Serialize (encode) import Data.Default.Class import Data.Colour import Data.Colour.Names@@ -37,14 +37,14 @@ blueRed :: [Colour Double] blueRed = interpolate 30 blue white ++ interpolate 30 white red -drawMatrix :: (MonadIO io, IOM.IOMatrix mat t Double) => mat t Double -> DrawOpt -> Source io L.ByteString+drawMatrix :: (MonadIO io, IOM.IOMatrix mat t Double) => mat t Double -> DrawOpt -> Source io B.ByteString drawMatrix mat opt = do yield pngSignature yield $ encode header yield . encode . preparePalette . coloursToPalette . _palette $ opt cs <- liftIO $ loop mat 0 $= toPngData $$ CL.consume- yield $ encode $ prepareIDatChunk $ L.fromChunks cs+ yield $ encode $ prepareIDatChunk $ B.concat cs yield $ encode endChunk where
src/BCM/Visualize/Internal.hs view
@@ -14,7 +14,6 @@ import Foreign.Storable( Storable, sizeOf ) import Data.Word import qualified Data.ByteString as B-import qualified Data.ByteString.Lazy as L import Data.Vector.Storable (Vector, unsafeToForeignPtr) import qualified Data.ByteString.Internal as S import qualified Data.Vector.Generic as G@@ -39,9 +38,9 @@ , interlaceMethod = PngNoInterlace } -prepareIDatChunk :: L.ByteString -> PngRawChunk+prepareIDatChunk :: B.ByteString -> PngRawChunk prepareIDatChunk imgData = PngRawChunk- { chunkLength = fromIntegral $ L.length imgData+ { chunkLength = fromIntegral $ B.length imgData , chunkType = iDATSignature , chunkCRC = pngComputeCrc [iDATSignature, imgData] , chunkData = imgData@@ -51,7 +50,7 @@ endChunk = PngRawChunk { chunkLength = 0 , chunkType = iENDSignature , chunkCRC = pngComputeCrc [iENDSignature]- , chunkData = L.empty+ , chunkData = B.empty } preparePalette :: Palette -> PngRawChunk@@ -61,7 +60,7 @@ , chunkCRC = pngComputeCrc [pLTESignature, binaryData] , chunkData = binaryData }- where binaryData = L.fromChunks [toByteString pal]+ where binaryData = B.concat [toByteString pal] toByteString :: forall a. (Storable a) => Vector a -> B.ByteString toByteString vec = S.PS (castForeignPtr ptr) offset (len * size)
src/BCM/Visualize/Internal/Types.hs view
@@ -8,26 +8,17 @@ #endif import Control.Monad (when)-import Data.Binary (Binary(..), Get)-import Data.Binary.Get( getWord8- , getWord32be- , getLazyByteString- )-import Data.Binary.Put( runPut- , putWord8- , putWord32be- , putLazyByteString- )+import Data.Serialize import Data.Bits( xor, (.&.), unsafeShiftR ) import qualified Data.Vector.Unboxed as U import Data.Word import Data.List (foldl')-import qualified Data.ByteString.Lazy as L-import qualified Data.ByteString.Lazy.Char8 as LS+import qualified Data.ByteString as B+import qualified Data.ByteString.Char8 as BS import Data.Vector.Storable (Vector) -- | Value used to identify a png chunk, must be 4 bytes long.-type ChunkSignature = L.ByteString+type ChunkSignature = B.ByteString type Palette = Vector Word8 @@ -48,7 +39,7 @@ { chunkLength :: Word32 , chunkType :: ChunkSignature , chunkCRC :: Word32- , chunkData :: L.ByteString+ , chunkData :: B.ByteString } -- | What kind of information is encoded in the IDAT section@@ -71,20 +62,20 @@ | PngInterlaceAdam7 deriving (Enum, Show) -instance Binary PngRawChunk where+instance Serialize PngRawChunk where put chunk = do putWord32be $ chunkLength chunk- putLazyByteString $ chunkType chunk+ putByteString $ chunkType chunk when (chunkLength chunk /= 0)- (putLazyByteString $ chunkData chunk)+ (putByteString $ chunkData chunk) putWord32be $ chunkCRC chunk get = do size <- getWord32be- chunkSig <- getLazyByteString (fromIntegral $ L.length iHDRSignature)+ chunkSig <- getByteString (fromIntegral $ B.length iHDRSignature) imgData <- if size == 0- then return L.empty- else getLazyByteString (fromIntegral size)+ then return B.empty+ else getByteString (fromIntegral size) crc <- getWord32be let computedCrc = pngComputeCrc [chunkSig, imgData]@@ -98,7 +89,7 @@ chunkType = chunkSig } -instance Binary PngImageType where+instance Serialize PngImageType where put PngGreyscale = putWord8 0 put PngTrueColour = putWord8 2 put PngIndexedColor = putWord8 3@@ -115,11 +106,11 @@ imageTypeOfCode 6 = return PngTrueColourWithAlpha imageTypeOfCode _ = fail "Invalid png color code" -instance Binary PngIHdr where+instance Serialize PngIHdr where put hdr = do putWord32be 13 let inner = runPut $ do- putLazyByteString iHDRSignature+ putByteString iHDRSignature putWord32be $ width hdr putWord32be $ height hdr putWord8 $ bitDepth hdr@@ -128,12 +119,12 @@ put $ filterMethod hdr put $ interlaceMethod hdr crc = pngComputeCrc [inner]- putLazyByteString inner+ putByteString inner putWord32be crc get = do _size <- getWord32be- ihdrSig <- getLazyByteString (L.length iHDRSignature)+ ihdrSig <- getByteString (B.length iHDRSignature) when (ihdrSig /= iHDRSignature) (fail "Invalid PNG file, wrong ihdr") w <- getWord32be@@ -154,7 +145,7 @@ interlaceMethod = interlace } -instance Binary PngInterlaceMethod where+instance Serialize PngInterlaceMethod where get = getWord8 >>= \w -> case w of 0 -> return PngNoInterlace 1 -> return PngInterlaceAdam7@@ -168,11 +159,11 @@ -- | Signature signalling that the following data will be a png image -- in the png bit stream pngSignature :: ChunkSignature-pngSignature = L.pack [137, 80, 78, 71, 13, 10, 26, 10]+pngSignature = B.pack [137, 80, 78, 71, 13, 10, 26, 10] -- | Helper function to help pack signatures. signature :: String -> ChunkSignature-signature = LS.pack +signature = BS.pack -- | Signature for the header chunk of png (must be the first) iHDRSignature :: ChunkSignature @@ -196,8 +187,8 @@ -- | Compute the CRC of a raw buffer, as described in annex D of the PNG -- specification.-pngComputeCrc :: [L.ByteString] -> Word32-pngComputeCrc = (0xFFFFFFFF `xor`) . L.foldl' updateCrc 0xFFFFFFFF . L.concat+pngComputeCrc :: [B.ByteString] -> Word32+pngComputeCrc = (0xFFFFFFFF `xor`) . B.foldl' updateCrc 0xFFFFFFFF . B.concat where updateCrc crc val = let u32Val = fromIntegral val lutVal = pngCrcTable U.! fromIntegral ((crc `xor` u32Val) .&. 0xFF)