packages feed

BCMtools (empty) → 0.1.0

raw patch · 14 files changed

+1399/−0 lines, 14 filesdep +BCMtoolsdep +basedep +binarysetup-changed

Dependencies added: BCMtools, base, binary, bytestring, bytestring-lexing, colour, conduit, conduit-extra, data-binary-ieee754, data-default-class, matrices, optparse-applicative, resourcet, split, transformers, unordered-containers, vector, zlib

Files

+ BCMtools.cabal view
@@ -0,0 +1,71 @@+name:                BCMtools+version:             0.1.0+synopsis:            Big Contact Map Tools+-- description:         +license:             MIT+license-file:        LICENSE+author:              Kai Zhang+maintainer:          kai@kzhang.org+copyright:           (c) 2015 Kai Zhang+category:            Bioinformatics+build-type:          Simple+-- extra-source-files:  +cabal-version:       >=1.10++library+  exposed-modules:     +    BCM+    BCM.Visualize+    BCM.Visualize.Internal+    BCM.Visualize.Internal.Types+    BCM.DiskMatrix+    BCM.IOMatrix+    BCM.Matrix.Instances++  -- other-modules:       +  build-depends:+      base >=4.0 && <5.0+    , bytestring >=0.9+    , data-binary-ieee754 >=0.4+    , binary >=0.7+    , conduit+    , conduit-extra+    , colour+    , data-default-class+    , unordered-containers+    , transformers+    , vector+    , zlib+    , matrices >=0.4.1++  hs-source-dirs:      src+  default-language:    Haskell2010++executable bcmtools+  hs-source-dirs:   exe+  main-is: Main.hs++  other-modules:+    BCMtools.Types+    BCMtools.Convert+    BCMtools.View++  build-depends:+      base >=4.0 && <5.0+    , BCMtools+    , binary+    , optparse-applicative+    , conduit+    , conduit-extra+    , resourcet+    , bytestring-lexing+    , bytestring+    , unordered-containers+    , split+    , data-default-class++  default-language:    Haskell2010++source-repository  head+  type: git+  location: https://github.com/kaizhang/BCMtools.git
+ LICENSE view
@@ -0,0 +1,20 @@+Copyright (c) 2015 Kai Zhang++Permission is hereby granted, free of charge, to any person obtaining+a copy of this software and associated documentation files (the+"Software"), to deal in the Software without restriction, including+without limitation the rights to use, copy, modify, merge, publish,+distribute, sublicense, and/or sell copies of the Software, and to+permit persons to whom the Software is furnished to do so, subject to+the following conditions:++The above copyright notice and this permission notice shall be included+in all copies or substantial portions of the Software.++THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,+EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF+MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.+IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY+CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,+TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE+SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ exe/BCMtools/Convert.hs view
@@ -0,0 +1,146 @@+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE FlexibleContexts #-}+module BCMtools.Convert+    ( convert+    , convertOptions+    ) where++import Control.Arrow ((&&&))+import qualified Data.ByteString.Char8 as B+import Data.List.Split (splitOn)+import qualified Data.HashMap.Strict as M+import Control.Monad.Trans.Resource (runResourceT)+import Data.Conduit+import qualified Data.Conduit.List as CL+import qualified Data.Conduit.Binary as Bin+import Data.Maybe (fromJust)+import Data.ByteString.Lex.Double (readDouble)+import Options.Applicative+import System.IO++import BCMtools.Types+import BCM (ContactMap, createContactMap, saveContactMap, closeContactMap)+import BCM.IOMatrix (DMatrix, MCSR, DSMatrix)+import BCM.Matrix.Instances ()++convertOptions :: Parser Command+convertOptions = fmap Convert $ ConvertOptions+             <$> strOption+                   ( long "genome"+                  <> short 'g'+                  <> metavar "ASSEMBLY"+                  <> help "e.g., hg19, or a file" )+             <*> fmap (splitOn ",") (strOption+                   ( long "rownames"+                  <> short 'r'+                  <> metavar "ROW LABELS" ))+             <*> fmap (splitOn ",") (strOption+                   ( long "colnames"+                  <> short 'c'+                  <> metavar "COLOUMN LABELS" ))+             <*> fmap readInt' (strOption+                   ( long "resolution"+                  <> short 's'+                  <> metavar "RESOLUTION" ))+             <*> switch+                   ( long "sparse"+                  <> help "whether to use sparse encoding" )+             <*> switch+                   ( long "symmetric"+                  <> help "whether to use symmetric encoding" )+  where+    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+                       | otherwise -> i++convert :: FilePath -> FilePath -> ConvertOptions -> IO ()+convert input output 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)+            _ -> error "Please check your input format"+        runner x = runResourceT $+                   Bin.sourceFile input $=+                   Bin.lines $= fn $$+                   createContactMap output rows cols (_resolution opt) (Just x)+        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+          | otherwise -> do+              cm <- runner n :: IO (ContactMap DMatrix)+              saveContactMap cm+              closeContactMap cm+  where+    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+        _ <- await+        CL.map f+      where+        f l = let [x1,x2,v] = B.split '\t' l+              in (chr1, readInt x1, chr2, readInt x2, readDouble' v)+    field5 = CL.map f+      where+        f l = let [x1,x2,x3,x4,x5] = B.split '\t' l+              in (x1, readInt x2, x3, readInt x4, readDouble' x5)+{-# INLINE convert #-}++readInt :: B.ByteString -> Int+readInt = fst . fromJust . B.readInt+{-# INLINE readInt #-}++readDouble' :: B.ByteString -> Double+readDouble' = fst . fromJust . readDouble+{-# INLINE readDouble' #-}++hg19 :: M.HashMap String Int+hg19 = M.fromList [ ("chr1", 249250621)+                  , ("chr2", 243199373)+                  , ("chr3", 198022430)+                  , ("chr4", 191154276)+                  , ("chr5", 180915260)+                  , ("chr6", 171115067)+                  , ("chr7", 159138663)+                  , ("chrX", 155270560)+                  , ("chr8", 146364022)+                  , ("chr9", 141213431)+                  , ("chr10", 135534747)+                  , ("chr11", 135006516)+                  , ("chr12", 133851895)+                  , ("chr13", 115169878)+                  , ("chr14", 107349540)+                  , ("chr15", 102531392)+                  , ("chr16", 90354753)+                  , ("chr17", 81195210)+                  , ("chr18", 78077248)+                  , ("chr20", 63025520)+                  , ("chrY", 59373566)+                  , ("chr19", 59128983)+                  , ("chr22", 51304566)+                  , ("chr21", 48129895)+                  ]
+ exe/BCMtools/Types.hs view
@@ -0,0 +1,29 @@+module BCMtools.Types+    ( BCMtoolsOptions(..)+    , Command(..)+    , ConvertOptions(..)+    , ViewOptions(..)+    ) where++data BCMtoolsOptions = BCMtoolsOptions+    { _input :: FilePath+    , _output :: FilePath+    , _command :: Command+    }++data Command = Convert ConvertOptions+             | View ViewOptions++data ConvertOptions = ConvertOptions+    { _genome :: String    -- ^ genome+    , _rownames :: [String]  -- ^ row chrs+    , _colnames :: [String]  -- ^ column chrs+    , _resolution :: Int   -- ^ resolution+    , _sparse :: Bool+    , _symmetric :: Bool+    } deriving (Show)++data ViewOptions = ViewOptions+    { _valueRange :: (Double, Double)+    , _inMem :: Bool+    }
+ exe/BCMtools/View.hs view
@@ -0,0 +1,84 @@+{-# LANGUAGE FlexibleContexts #-}+module BCMtools.View+    ( view+    , viewOptions+    ) where++import qualified Data.ByteString.Lazy as L+import Data.Binary.Get (runGet, getWord32le)+import System.IO+import Options.Applicative+import Data.Conduit+import qualified Data.Conduit.List as CL+import qualified Data.Conduit.Binary as Bin+import Control.Monad.Trans.Resource (runResourceT)+import Data.Default.Class (def)++import BCMtools.Types+import BCM (ContactMap, _matrix, closeContactMap, openContactMap)+import BCM.IOMatrix (DMatrix, MCSR, DSMatrix, MMatrix, MSMatrix)+import BCM.Matrix.Instances+import BCM.Visualize++viewOptions :: Parser Command +viewOptions = fmap View $ ViewOptions+          <$> fmap f (strOption+                ( 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+        hSeek h AbsoluteSeek 4+        p <- fromIntegral . runGet getWord32le <$> L.hGet h 4+        hSeek h AbsoluteSeek p+        runGet getWord32le <$> L.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+          | 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+                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+          | otherwise -> error "unknown format"+  where+    draw x = runResourceT $ drawMatrix (_matrix x) drawopt $= CL.map L.toStrict $$ Bin.sinkFile output+    drawopt = def { _range = _valueRange opt }+{-# INLINE view #-}+
+ exe/Main.hs view
@@ -0,0 +1,41 @@+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE FlexibleContexts #-}+module Main where++import Options.Applicative++import BCMtools.Types+import BCMtools.Convert (convert, convertOptions)+import BCMtools.View (view, viewOptions)++globalOptions :: Parser Command -> Parser BCMtoolsOptions+globalOptions cmd = BCMtoolsOptions +             <$> argument str (metavar "INPUT") +             <*> strOption+                   ( long "output"+                  <> short 'o'+                  <> metavar "OUTPUT" )+             <*> cmd++bcmtoolsOptions :: Parser BCMtoolsOptions+bcmtoolsOptions = subparser $+      command "convert" ( info (helper <*> globalOptions convertOptions) $+                               fullDesc+                            <> progDesc "file conversion"+                        )+   <> command "view"    ( info (helper <*> globalOptions viewOptions) $+                               fullDesc+                            <> progDesc "view file"+                        )++runBCMtools :: BCMtoolsOptions -> IO () +runBCMtools (BCMtoolsOptions input output bcmopt) = case bcmopt of+    Convert opt -> convert input output opt+    View opt -> view input output opt++main :: IO ()+main = execParser opts >>= runBCMtools+  where+    opts = info (helper <*> bcmtoolsOptions)+         ( fullDesc+        <> header "Big Contact Map (BCM) tools" )
+ src/BCM.hs view
@@ -0,0 +1,143 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE TypeFamilies #-}+module BCM+    ( ContactMap(..)+    , createContactMap+    , saveContactMap+    , openContactMap+    , closeContactMap+    ) where++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.Char8 as B+import qualified Data.HashMap.Strict as M+import Data.Binary.Put+import Data.Binary.Get+import Data.Conduit+import qualified Data.Conduit.List as CL+import System.IO+import Data.List (foldl')+import Data.Word (Word32)+import Text.Printf (printf)++import qualified BCM.DiskMatrix as DM+import qualified BCM.IOMatrix as IOM++-- contact map binary format+-- 4 bytes magic + 4 byte Int (matrix start) + 4 bytes Int (step) + chroms + 1 bytes (reserve) + matrix+data ContactMap m = ContactMap+    { _rowLabels :: M.HashMap B.ByteString (Int, Int)+    , _colLabels :: M.HashMap B.ByteString (Int, Int)+    , _resolution :: Int+    , _matrix :: m+    , _handle :: Handle+    }+++contact_map_magic :: Word32+contact_map_magic = 0x9921ABF0++createContactMap :: (IOM.IOMatrix m t Double, MonadIO io, mat ~ m t Double)+                 => FilePath+                 -> [(B.ByteString, Int)]+                 -> [(B.ByteString, Int)]+                 -> Int+                 -> Maybe Int+                 -> Sink (B.ByteString, Int, B.ByteString, Int, Double) io (ContactMap mat)+createContactMap fl rowChr colChr res len = do+    h <- liftIO $ openFile fl ReadWriteMode++    let source  = CL.mapM $ \(chr1, i, chr2, j, v) -> do+            let (p1, s1) = M.lookupDefault errMsg chr1 rLab+                (p2, s2) = M.lookupDefault errMsg chr2 cLab+                i' = i `div` res + p1+                j' = j `div` res + p2+            when (i > s1 || j > s2) $ error "createContactMap: Index out of bounds"+            when (i `mod` res /= 0 || j `mod` res /=0) $+                liftIO $ hPutStrLn stderr $ printf "(%d,%d) is not divisible by %d" i j res+            return ((i',j'), v)+++    liftIO $ L.hPut h $ L.replicate (fromIntegral offset) 0++    m <- source $= IOM.hCreateMatrix h (r,c) len++    return $ ContactMap rLab cLab res m h+  where+    r = foldl' (\acc (_,x) -> acc + (x - 1) `div` res + 1)  0 rowChr+    c = foldl' (\acc (_,x) -> acc + (x - 1) `div` res + 1)  0 colChr+    rLab = mkLabels rowChr res+    cLab = mkLabels colChr res+    nByte x = let n1 = foldl' (+) 0 $ map B.length $ fst $ unzip x+                  n2 = 16 * n3+                  n3 = length x+               in n1 + n2 + n3+    offset = 4 + 4 + 4 + nByte rowChr + nByte colChr + 2 + 1+    errMsg = error "createContactMap: Unknown chromosome"++saveContactMap :: (IOM.IOMatrix m t a, mat ~ m t a) => ContactMap mat -> IO ()+saveContactMap (ContactMap rowChr colChr res mat handle) = do+        hSeek handle AbsoluteSeek 0+        L.hPutStr handle . runPut . putWord32le $ contact_map_magic+        L.hPutStr handle . runPut . putWord32le $ offset+        L.hPutStr handle . runPut . putWord32le . fromIntegral $ res+        L.hPutStr handle rowAndcol+        L.hPutStr handle . runPut . putWord8 $ 0+        IOM.hSaveMatrix handle mat+      where+        rowAndcol = L.concat [rows, "\0", cols, "\0"]+        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+        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+    guard $ magic == contact_map_magic+    _ <- runGet getWord32le <$> L.hGet h 4+    res <- 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+    return $ ContactMap rows cols res mat h+  where+    getChrs acc h = do+        chr <- getByteStringNul h+        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+               getChrs ((chr, (a, b)) : acc) h+    getByteStringNul h = B.concat <$> go []+      where+        go acc = do+            x <- B.hGet h 1+            case x of+                "\0" -> return acc+                _ -> go $ acc ++ [x]++closeContactMap :: ContactMap mat -> IO ()+closeContactMap cm = hClose $ _handle cm+++-------------------------------------------------------------------------------+-- Helper functions+-------------------------------------------------------------------------------++mkLabels :: [(B.ByteString, Int)] -> Int -> M.HashMap B.ByteString (Int, Int)+mkLabels xs step = M.fromList $ foldr f [] xs+  where+    f (chr, size) [] = [(chr, (0, size))]+    f (chr, size) acc@((_, (a, b)) : _) = (chr, (a + ((b - 1) `div` step + 1), size)) : acc+{-# INLINE mkLabels #-}
+ src/BCM/DiskMatrix.hs view
@@ -0,0 +1,212 @@+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE ScopedTypeVariables #-}+module BCM.DiskMatrix+    ( Offset+    , DiskData(..)+    , DiskMatrix(..)+    , read+    , write+    , DMatrix(..)+    , DSMatrix(..)+    ) where++import Prelude hiding (replicate)+import Control.Monad (replicateM_)+import Control.Monad.IO.Class (MonadIO(..))+import Control.Applicative ((<$>))+import qualified Data.ByteString.Lazy as L+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 qualified Data.Vector.Generic as G+import System.IO+import Numeric (showHex)++import BCM.Matrix.Instances (d_matrix_magic, ds_matrix_magic)++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+    elemType :: m a -> a+    elemType _ = undefined+    {-# INLINE elemType #-}++    hReadMatrixEither :: MonadIO io => Handle -> io (Either String (m a))++    dim :: m a -> (Int, Int)++    replicate :: MonadIO io => Handle -> (Int, Int) -> a -> io (m a)++    unsafeRead :: MonadIO io => m a -> (Int, Int) -> io a++    unsafeWrite :: MonadIO io => m a -> (Int, Int) -> a -> io ()++    unsafeReadRow :: (G.Vector v a, MonadIO io) => m a -> Int -> io (v a)+    unsafeReadRow mat i = G.generateM c $ \j -> unsafeRead mat (i,j)+      where+        (_,c) = dim mat+    {-# INLINE unsafeReadRow #-}++    close :: MonadIO io => m a -> io ()++-- Derived methods++write :: (MonadIO io, DiskMatrix m a) => m a -> (Int, Int) -> a -> io ()+write mat (i,j) x | i >= r || j >= c = error "Index out of bounds"+                  | otherwise = unsafeWrite mat (i,j) x+  where+    (r,c) = dim mat+{-# INLINE write #-}++data DMatrix a = DMatrix !Int  -- ^ rows+                         !Int  -- ^ cols+                         !Offset -- offset+                         !Handle  -- ^ file handle++instance DiskData 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 ""+    {-# INLINE hReadMatrixEither #-}++    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+    {-# 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+    {-# 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+    {-# INLINE unsafeWrite #-}++    unsafeReadRow mat@(DMatrix _ c offset h) i = liftIO $ do+        hSeek h AbsoluteSeek $ offset + fromIntegral (size * c * i)+        G.replicateM c $ hRead1 h+      where+        size = sizeOf $ elemType mat+    {-# INLINE unsafeReadRow #-}++    close (DMatrix _ _ _ h) = liftIO $ hClose h++-- | Symmetric matrix+data DSMatrix a = DSMatrix !Int  -- ^ size+                           !Offset  -- ^ offset+                           !Handle  -- ^ file handle++instance DiskData 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 ""+    {-# INLINE hReadMatrixEither #-}++    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+    {-# 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+    {-# 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+    {-# INLINE unsafeWrite #-}++    close (DSMatrix _ _ h) = liftIO $ hClose h++------------------------------------------------------------------------------+-- helper functions+------------------------------------------------------------------------------++-- normal matrix indexing+idx :: Int -> Int -> Int -> Int+idx c i j = i * c + j+{-# INLINE idx #-}++-- upper triangular matrix indexing+idx' :: Int -> Int -> Int -> Int+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
@@ -0,0 +1,172 @@+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE KindSignatures #-}+module BCM.IOMatrix+    ( IOMatrix(..)+    , DMatrix+    , DSMatrix+    , MMatrix+    , MSMatrix+    , MCSR+    ) where++import Control.Monad (when, forM_)+import Control.Applicative ((<$>))+import Control.Monad.IO.Class (MonadIO(..))+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+import Data.Matrix.Dense.Generic (Matrix(..))+import Data.Matrix.Sparse.Generic (CSR(..), Zero(..))+import Data.Matrix.Symmetric (SymMatrix(..))++import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Unboxed.Mutable as UM++import Data.Conduit (Sink)+import qualified Data.Conduit.List as CL+import Text.Printf (printf)+import System.IO++import qualified BCM.DiskMatrix as DM++type MMatrix = IOMat (Matrix U.Vector) Double+type MSMatrix = IOMat (SymMatrix U.Vector) Double+type MCSR = IOSMat (CSR U.Vector) Double++type DMatrix = IODMat DM.DMatrix Double+type DSMatrix = IODMat DM.DSMatrix Double++class IOMatrix mat (t :: * -> *) a where+    dim :: mat t a -> (Int, Int)++    unsafeIndexM :: MonadIO io => mat t a -> (Int, Int) -> io a+    +    unsafeTakeRowM :: (U.Unbox a, MonadIO io) => mat t a -> Int -> io (U.Vector a)+    +    -- | Read a matrix from file handle+    hReadMatrix :: MonadIO io+                => Handle+                -> io (mat t a)++    hSaveMatrix :: MonadIO io => Handle -> mat t a -> io ()+    hSaveMatrix _ _ = return ()++    hCreateMatrix :: MonadIO io+                  => Handle  -- ^ file handle+                  -> (Int, Int)  -- ^ matrix dimension+                  -> Maybe Int  -- ^ number of non-zero elements+                  -> Sink ((Int,Int), a) io (mat t a)++-- | Just a wrapper+newtype IODMat m a = IODMat { unwrapD :: m a }++instance (DM.DiskData a, DM.DiskMatrix m a) => IOMatrix IODMat m a where+    dim = DM.dim . unwrapD+    {-# INLINE dim #-}++    unsafeIndexM = DM.unsafeRead . unwrapD+    {-# INLINE unsafeIndexM #-}++    unsafeTakeRowM = DM.unsafeReadRow . unwrapD+    {-# INLINE unsafeTakeRowM #-}++    hReadMatrix handle = do +        r <- DM.hReadMatrixEither handle+        case r of+            Left e -> error e+            Right m -> return $ IODMat m+    {-# INLINE hReadMatrix #-}++    hCreateMatrix handle (r,c) _ = do+        mat <- DM.replicate handle (r,c) DM.zero+        CL.mapM_ $ \((i,j), x) -> DM.unsafeWrite mat (i,j) x+        return $ IODMat mat+    {-# INLINE hCreateMatrix #-}++-- | Just a wrapper+newtype IOMat m a = IOMat { unwrap :: m a }++instance ( U.Unbox a+         , DM.DiskData a+         , Binary (m U.Vector a)+         , MG.Matrix m U.Vector a+         ) => IOMatrix IOMat (m U.Vector) a where+    dim = MG.dim . unwrap+    {-# INLINE dim #-}++    unsafeIndexM (IOMat mat) (i,j) = return $ MG.unsafeIndex mat (i,j)+    {-# INLINE unsafeIndexM #-}++    unsafeTakeRowM (IOMat mat) i = return $ MG.takeRow mat i+    {-# INLINE unsafeTakeRowM #-}++    hReadMatrix handle = liftIO $ do+        mat <- decode <$> L.hGetContents handle+        return $ IOMat mat+    {-# INLINE hReadMatrix #-}++    hSaveMatrix handle (IOMat mat) = liftIO . L.hPutStr handle . encode $ mat+    {-# INLINE hSaveMatrix #-}++    hCreateMatrix _ (r,c) _ = do+        mat <- liftIO $ MGM.replicate (r,c) DM.zero+        CL.mapM_ $ \((i,j), x) -> liftIO $ MGM.unsafeWrite mat (i,j) x+        mat' <- liftIO $ MG.unsafeFreeze mat+        return $ IOMat mat'+    {-# INLINE hCreateMatrix #-}++-- | Just a wrapper+newtype IOSMat m a = IOSMat { unwrapS :: m a }++instance ( Zero a+         , U.Unbox a+         , DM.DiskData a+         , Binary (CSR U.Vector a)+         ) => IOMatrix IOSMat (CSR U.Vector) a where+    dim = MG.dim . unwrapS+    {-# INLINE dim #-}++    unsafeIndexM (IOSMat mat) (i,j) = return $ (MG.!) mat (i,j)+    {-# INLINE unsafeIndexM #-}++    unsafeTakeRowM (IOSMat mat) i = return $ MG.takeRow mat i+    {-# INLINE unsafeTakeRowM #-}++    hReadMatrix handle = liftIO $ do+        mat <- decode <$> L.hGetContents handle+        return $ IOSMat mat+    {-# INLINE hReadMatrix #-}++    hSaveMatrix handle (IOSMat mat) = liftIO . L.hPutStr handle . encode $ mat+    {-# INLINE hSaveMatrix #-}++    hCreateMatrix _ (r,c) (Just n) = do+        v <- liftIO $ UM.new n+        col <- liftIO $ UM.new n+        row <- liftIO $ UM.new (r+1)++        ((i,_),_) <- flip CL.foldM ((-1,-1),0) $ \((i',j'), acc) ((i,j),x) ->+            if i > i' || (i == i' && j > j')+               then liftIO $ do+                   UM.write v acc x+                   UM.write col acc j+                   let stride = i - i'+                   when (stride > 0) $ forM_ [0..stride-1] $ \s -> UM.write row (i-s) acc+            +                   return ((i,j), acc+1)+               else error $ printf "Input must be sorted by row and then by column: (%d,%d) >= (%d,%d)" i' j' i j++        let stride = r - i+        liftIO $ forM_ [0..stride-1] $ \s -> UM.write row (r-s) n++        v' <- liftIO $ U.unsafeFreeze v+        col' <- liftIO $ U.unsafeFreeze col+        row' <- liftIO $ U.unsafeFreeze row++        return $ IOSMat $ CSR r c v' col' row'+    hCreateMatrix _ _ _ = error "no length info available"+    {-# INLINE hCreateMatrix #-}
+ src/BCM/Matrix/Instances.hs view
@@ -0,0 +1,112 @@+{-# 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
@@ -0,0 +1,74 @@+{-# LANGUAGE FlexibleContexts #-}+module BCM.Visualize+    ( DrawOpt(..)+    , reds+    , blueRed+    , drawMatrix+    ) where++import Control.Monad.IO.Class+import qualified Data.ByteString.Lazy as L+import Data.Binary (encode)+import Data.Default.Class+import Data.Colour+import Data.Colour.Names+import Data.Conduit+import qualified Data.Conduit.List as CL+import qualified Data.Vector.Unboxed as U++import BCM.Visualize.Internal+import BCM.Visualize.Internal.Types+import qualified BCM.IOMatrix as IOM++data DrawOpt = DrawOpt +    { _range :: !(Double, Double)+    , _palette :: ![Colour Double]+    }++instance Default DrawOpt where+    def = DrawOpt+        { _range = (0,10)+        , _palette = reds+        }++reds :: [Colour Double]+reds = interpolate 62 white red++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 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 endChunk+  where+    loop m i+      | i < h = do+          row <- IOM.unsafeTakeRowM m i+          yield $ U.toList $ U.map drawPixel row+          loop m (i+1)+      | otherwise = return ()++    drawPixel x | x <= lo = 0+                | x >= hi = fromIntegral $ n - 1+                | otherwise = truncate $ (x - lo) / step+    (w,h) = IOM.dim mat+    (lo,hi) = _range opt+    step = (hi - lo) / fromIntegral n+    n = length $ _palette opt+    header = preparePngHeader w h PngIndexedColor 8+{-# INLINE drawMatrix #-}++interpolate :: Int -> Colour Double -> Colour Double -> [Colour Double]+interpolate n c1 c2 = loop 1+  where+    loop i | i > n = []+           | otherwise = blend (fromIntegral i * step) c2 c1 : loop (i+1)+    step = 1 / fromIntegral (n+1)+{-# INLINE interpolate #-}
+ src/BCM/Visualize/Internal.hs view
@@ -0,0 +1,81 @@+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE CPP #-}+-- most of the codes in this file are directly copied from JuicyPixel++module BCM.Visualize.Internal where++#if !MIN_VERSION_base(4,8,0)+import Foreign.ForeignPtr.Safe( ForeignPtr, castForeignPtr )+#else+import Foreign.ForeignPtr( ForeignPtr, castForeignPtr )+#endif++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+import Data.Colour+import Data.Colour.SRGB++import Data.Conduit.Zlib as Z++import Data.Conduit+import qualified Data.Conduit.List as CL++import BCM.Visualize.Internal.Types++preparePngHeader :: Int -> Int -> PngImageType -> Word8 -> PngIHdr+preparePngHeader w h imgType depth = PngIHdr+    { width             = fromIntegral w+    , height            = fromIntegral h+    , bitDepth          = depth+    , colourType        = imgType+    , compressionMethod = 0+    , filterMethod      = 0+    , interlaceMethod   = PngNoInterlace+    }++prepareIDatChunk :: L.ByteString -> PngRawChunk+prepareIDatChunk imgData = PngRawChunk+    { chunkLength = fromIntegral $ L.length imgData+    , chunkType   = iDATSignature+    , chunkCRC    = pngComputeCrc [iDATSignature, imgData]+    , chunkData   = imgData+    }++endChunk :: PngRawChunk+endChunk = PngRawChunk { chunkLength = 0+                       , chunkType = iENDSignature+                       , chunkCRC = pngComputeCrc [iENDSignature]+                       , chunkData = L.empty+                       }++preparePalette :: Palette -> PngRawChunk+preparePalette pal = PngRawChunk+  { chunkLength = fromIntegral $ G.length pal+  , chunkType   = pLTESignature+  , chunkCRC    = pngComputeCrc [pLTESignature, binaryData]+  , chunkData   = binaryData+  }+   where binaryData = L.fromChunks [toByteString pal]++toByteString :: forall a. (Storable a) => Vector a -> B.ByteString+toByteString vec = S.PS (castForeignPtr ptr) offset (len * size)+  where (ptr, offset, len) = unsafeToForeignPtr vec+        size = sizeOf (undefined :: a)+{-# INLINE toByteString #-}++coloursToPalette :: [Colour Double] -> Palette+coloursToPalette = G.fromList . concatMap f+  where+    f c = let RGB r g b = toSRGB24 c+          in [r,g,b]+{-# INLINE coloursToPalette #-}++toPngData :: Conduit [Word8] IO B.ByteString+toPngData = CL.map (B.pack . (0:)) $= Z.compress 5 Z.defaultWindowBits+{-# INLINE toPngData #-}
+ src/BCM/Visualize/Internal/Types.hs view
@@ -0,0 +1,212 @@+{-# LANGUAGE CPP #-}+module BCM.Visualize.Internal.Types where++#if !MIN_VERSION_base(4,8,0)+import Foreign.ForeignPtr.Safe( ForeignPtr, castForeignPtr )+#else+import Foreign.ForeignPtr( ForeignPtr, castForeignPtr )+#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.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 Data.Vector.Storable (Vector)++-- | Value used to identify a png chunk, must be 4 bytes long.+type ChunkSignature = L.ByteString++type Palette = Vector Word8++-- | Generic header used in PNG images.+data PngIHdr = PngIHdr+    { width             :: !Word32       -- ^ Image width in number of pixel+    , height            :: !Word32       -- ^ Image height in number of pixel+    , bitDepth          :: !Word8        -- ^ Number of bit per sample+    , colourType        :: !PngImageType -- ^ Kind of png image (greyscale, true color, indexed...)+    , compressionMethod :: !Word8        -- ^ Compression method used+    , filterMethod      :: !Word8        -- ^ Must be 0+    , interlaceMethod   :: !PngInterlaceMethod   -- ^ If the image is interlaced (for progressive rendering)+    }+    deriving Show++-- | Data structure during real png loading/parsing+data PngRawChunk = PngRawChunk+    { chunkLength :: Word32+    , chunkType   :: ChunkSignature+    , chunkCRC    :: Word32+    , chunkData   :: L.ByteString+    }++-- | What kind of information is encoded in the IDAT section+-- of the PngFile+data PngImageType =+      PngGreyscale+    | PngTrueColour+    | PngIndexedColor+    | PngGreyscaleWithAlpha+    | PngTrueColourWithAlpha+    deriving Show++-- | Different known interlace methods for PNG image+data PngInterlaceMethod =+      -- | No interlacing, basic data ordering, line by line+      -- from left to right.+      PngNoInterlace++      -- | Use the Adam7 ordering, see `adam7Reordering`+    | PngInterlaceAdam7+    deriving (Enum, Show)++instance Binary PngRawChunk where+    put chunk = do+        putWord32be $ chunkLength chunk+        putLazyByteString $ chunkType chunk+        when (chunkLength chunk /= 0)+             (putLazyByteString $ chunkData chunk)+        putWord32be $ chunkCRC chunk++    get = do+        size <- getWord32be+        chunkSig <- getLazyByteString (fromIntegral $ L.length iHDRSignature)+        imgData <- if size == 0+            then return L.empty+            else getLazyByteString (fromIntegral size)+        crc <- getWord32be++        let computedCrc = pngComputeCrc [chunkSig, imgData]+        when (computedCrc `xor` crc /= 0)+             (fail $ "Invalid CRC : " ++ show computedCrc ++ ", "+                                      ++ show crc)+        return PngRawChunk {+            chunkLength = size,+            chunkData = imgData,+            chunkCRC = crc,+            chunkType = chunkSig+        }++instance Binary PngImageType where+    put PngGreyscale = putWord8 0+    put PngTrueColour = putWord8 2+    put PngIndexedColor = putWord8 3+    put PngGreyscaleWithAlpha = putWord8 4+    put PngTrueColourWithAlpha = putWord8 6++    get = get >>= imageTypeOfCode++imageTypeOfCode :: Word8 -> Get PngImageType+imageTypeOfCode 0 = return PngGreyscale+imageTypeOfCode 2 = return PngTrueColour+imageTypeOfCode 3 = return PngIndexedColor+imageTypeOfCode 4 = return PngGreyscaleWithAlpha+imageTypeOfCode 6 = return PngTrueColourWithAlpha+imageTypeOfCode _ = fail "Invalid png color code"++instance Binary PngIHdr where+    put hdr = do+        putWord32be 13+        let inner = runPut $ do+                putLazyByteString iHDRSignature+                putWord32be $ width hdr+                putWord32be $ height hdr+                putWord8    $ bitDepth hdr+                put $ colourType hdr+                put $ compressionMethod hdr+                put $ filterMethod hdr+                put $ interlaceMethod hdr+            crc = pngComputeCrc [inner]+        putLazyByteString inner+        putWord32be crc++    get = do+        _size <- getWord32be+        ihdrSig <- getLazyByteString (L.length iHDRSignature)+        when (ihdrSig /= iHDRSignature)+             (fail "Invalid PNG file, wrong ihdr")+        w <- getWord32be+        h <- getWord32be+        depth <- get+        colorType <- get+        compression <- get+        filtermethod <- get+        interlace <- get+        _crc <- getWord32be+        return PngIHdr {+            width = w,+            height = h,+            bitDepth = depth,+            colourType = colorType,+            compressionMethod = compression,+            filterMethod = filtermethod,+            interlaceMethod = interlace+        }++instance Binary PngInterlaceMethod where+    get = getWord8 >>= \w -> case w of+        0 -> return PngNoInterlace+        1 -> return PngInterlaceAdam7+        _ -> fail "Invalid interlace method"++    put PngNoInterlace    = putWord8 0+    put PngInterlaceAdam7 = putWord8 1++-- signature++-- | 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]++-- | Helper function to help pack signatures.+signature :: String -> ChunkSignature+signature = LS.pack ++-- | Signature for the header chunk of png (must be the first)+iHDRSignature :: ChunkSignature +iHDRSignature = signature "IHDR"++-- | Signature for a palette chunk in the pgn file. Must+-- occure before iDAT.+pLTESignature :: ChunkSignature+pLTESignature = signature "PLTE"++-- | Signature for a data chuck (with image parts in it)+iDATSignature :: ChunkSignature+iDATSignature = signature "IDAT"++-- | Signature for the last chunk of a png image, telling+-- the end.+iENDSignature :: ChunkSignature+iENDSignature = signature "IEND"++----------------------------------------------------------------------------++-- | 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+    where updateCrc crc val =+              let u32Val = fromIntegral val+                  lutVal = pngCrcTable U.! fromIntegral ((crc `xor` u32Val) .&. 0xFF)+              in lutVal `xor` (crc `unsafeShiftR` 8)++-- | From the Annex D of the png specification.+pngCrcTable :: U.Vector Word32+pngCrcTable = U.fromListN 256 [ foldl' updateCrcConstant c [zero .. 7] | c <- [0 .. 255] ]+    where zero = 0 :: Int -- To avoid defaulting to Integer+          updateCrcConstant c _ | c .&. 1 /= 0 = magicConstant `xor` (c `unsafeShiftR` 1)+                                | otherwise = c `unsafeShiftR` 1+          magicConstant = 0xedb88320 :: Word32