diff --git a/BCMtools.cabal b/BCMtools.cabal
new file mode 100644
--- /dev/null
+++ b/BCMtools.cabal
@@ -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
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -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.
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/exe/BCMtools/Convert.hs b/exe/BCMtools/Convert.hs
new file mode 100644
--- /dev/null
+++ b/exe/BCMtools/Convert.hs
@@ -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)
+                  ]
diff --git a/exe/BCMtools/Types.hs b/exe/BCMtools/Types.hs
new file mode 100644
--- /dev/null
+++ b/exe/BCMtools/Types.hs
@@ -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
+    }
diff --git a/exe/BCMtools/View.hs b/exe/BCMtools/View.hs
new file mode 100644
--- /dev/null
+++ b/exe/BCMtools/View.hs
@@ -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 #-}
+
diff --git a/exe/Main.hs b/exe/Main.hs
new file mode 100644
--- /dev/null
+++ b/exe/Main.hs
@@ -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" )
diff --git a/src/BCM.hs b/src/BCM.hs
new file mode 100644
--- /dev/null
+++ b/src/BCM.hs
@@ -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 #-}
diff --git a/src/BCM/DiskMatrix.hs b/src/BCM/DiskMatrix.hs
new file mode 100644
--- /dev/null
+++ b/src/BCM/DiskMatrix.hs
@@ -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' #-}
diff --git a/src/BCM/IOMatrix.hs b/src/BCM/IOMatrix.hs
new file mode 100644
--- /dev/null
+++ b/src/BCM/IOMatrix.hs
@@ -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 #-}
diff --git a/src/BCM/Matrix/Instances.hs b/src/BCM/Matrix/Instances.hs
new file mode 100644
--- /dev/null
+++ b/src/BCM/Matrix/Instances.hs
@@ -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 #-}
diff --git a/src/BCM/Visualize.hs b/src/BCM/Visualize.hs
new file mode 100644
--- /dev/null
+++ b/src/BCM/Visualize.hs
@@ -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 #-}
diff --git a/src/BCM/Visualize/Internal.hs b/src/BCM/Visualize/Internal.hs
new file mode 100644
--- /dev/null
+++ b/src/BCM/Visualize/Internal.hs
@@ -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 #-}
diff --git a/src/BCM/Visualize/Internal/Types.hs b/src/BCM/Visualize/Internal/Types.hs
new file mode 100644
--- /dev/null
+++ b/src/BCM/Visualize/Internal/Types.hs
@@ -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