{-
mkContactMap' :: FilePath -> [BED3] -> Int -> Int -> Source IO ..
mkContactMap' bamFl regions w extend = do
    handle <- liftIO $ BI.open bamFl
    forM_ [0..n-1] 
  where
    positions = V.scanl f 0 regions'
    n = V.last positions
    f acc (BED3 chr s e) = acc + (e - s) `div` w
    regions' = V.fromList regions
    g i = binarySearch positions i
    -}

-- | O(n * (logN + k)). n = number of bins, N = number of tags. Generate contanct
-- map using constant memory
mkContactMap :: FilePath -> BED3 -> Int -> Int -> Source IO ((Int, Int), Int)
mkContactMap bamFl (BED3 chr s e) w extend = do
    handle <- lift $ BI.open bamFl
    forM_ [0..n-1] $ \i ->
        forM_ [i..n-1] $ \j -> do
            let s1 = s + i * w
                c1 = s1 + w'
                s2 = s + j * w
                c2 = s2 + w'
            r <- lift $ readCount handle (w'+extend) ((chr, c1), (chr, c2))
            if i == j
               then yield ((s1, s2), r `div` 2)
               else yield ((s1, s2), r)
  where
    n = (e - s) `div` w
    w' = w `div` 2

-- | O(N + n). Store matrix in memory, faster when region is big
mkContactMap :: FilePath -> BED3 -> Int -> Int -> Source IO ((Int, Int), Int)
mkContactMap bamFl (BED3 chr s e) w extend = do
    handle <- lift $ BI.open bamFl
    forM_ [0..n-1] $ \i ->
        forM_ [i..n-1] $ \j -> do
            let s1 = s + i * w
                c1 = s1 + w'
                s2 = s + j * w
                c2 = s2 + w'
            r <- lift $ readCount handle (w'+extend) ((chr, c1), (chr, c2))
            if i == j
               then yield ((s1, s2), r `div` 2)
               else yield ((s1, s2), r)
  where
    n = (e - s) `div` w
    w' = w `div` 2

-- | O(N + n). Store matrix in memory, faster when region is big
mkContactMap bamFl (BED3 chr s e) w extend = do
    handle <- lift $ BI.open bamFl
    forM_ [0..n-1] $ \i ->
        forM_ [i..n-1] $ \j -> do
            let s1 = s + i * w
                c1 = s1 + w'
                s2 = s + j * w
                c2 = s2 + w'
            r <- lift $ readCount handle (w'+extend) ((chr, c1), (chr, c2))
            if i == j
               then yield ((s1, s2), r `div` 2)
               else yield ((s1, s2), r)
  where
    n = (e - s) `div` w
    w' = w `div` 2

-- | O(N + n). Store matrix in memory, faster when region is big