diff --git a/Numeric/SGD.hs b/Numeric/SGD.hs
deleted file mode 100644
--- a/Numeric/SGD.hs
+++ /dev/null
@@ -1,136 +0,0 @@
-{-# LANGUAGE RecordWildCards #-}
-
-
--- | Stochastic gradient descent implementation using mutable
--- vectors for efficient update of the parameters vector.
--- A user is provided with the immutable vector of parameters
--- so he is able to compute the gradient outside of the IO monad.
--- Currently only the Gaussian priors are implemented.
---
--- This is a preliminary version of the SGD library and API may change
--- in future versions.
-
-
-module Numeric.SGD
-( SgdArgs (..)
-, sgdArgsDefault
-, Para
-, sgd
-, module Numeric.SGD.Grad
-, module Numeric.SGD.Dataset
-) where
-
-
-import           Control.Monad (forM_)
-import qualified System.Random as R
-import qualified Data.Vector.Unboxed as U
-import qualified Data.Vector.Unboxed.Mutable as UM
-import qualified Control.Monad.Primitive as Prim
-
-import           Numeric.SGD.Grad
-import           Numeric.SGD.Dataset
-
-
--- | SGD parameters controlling the learning process.
-data SgdArgs = SgdArgs
-    { -- | Size of the batch
-      batchSize :: Int
-    -- | Regularization variance
-    , regVar    :: Double
-    -- | Number of iterations
-    , iterNum   :: Double
-    -- | Initial gain parameter
-    , gain0     :: Double
-    -- | After how many iterations over the entire dataset
-    -- the gain parameter is halved
-    , tau       :: Double }
-
-
--- | Default SGD parameter values.
-sgdArgsDefault :: SgdArgs
-sgdArgsDefault = SgdArgs
-    { batchSize = 30
-    , regVar    = 10
-    , iterNum   = 10
-    , gain0     = 1
-    , tau       = 5 }
-
-
--- | Vector of parameters.
-type Para       = U.Vector Double 
-
-
--- | Type synonym for mutable vector with Double values.
-type MVect      = UM.MVector (Prim.PrimState IO) Double
-
-
--- | A stochastic gradient descent method.
--- A notification function can be used to provide user with
--- information about the progress of the learning.
-sgd
-    :: SgdArgs                  -- ^ SGD parameter values
-    -> (Para -> Int -> IO ())   -- ^ Notification run every update
-    -> (Para -> x -> Grad)      -- ^ Gradient for dataset element
-    -> Dataset x                -- ^ Dataset
-    -> Para                     -- ^ Starting point
-    -> IO Para                  -- ^ SGD result
-sgd SgdArgs{..} notify mkGrad dataset x0 = do
-    u <- UM.new (U.length x0)
-    doIt u 0 (R.mkStdGen 0) =<< U.thaw x0
-  where
-    -- Gain in k-th iteration.
-    gain k = (gain0 * tau) / (tau + done k)
-
-    -- Number of completed iterations over the full dataset.
-    done k
-        = fromIntegral (k * batchSize)
-        / fromIntegral (size dataset)
-
-    doIt u k stdGen x
-      | done k > iterNum = do
-        frozen <- U.unsafeFreeze x
-        notify frozen k
-        return frozen
-      | otherwise = do
-        (batch, stdGen') <- sample stdGen batchSize dataset
-
-        -- Freeze mutable vector of parameters. The frozen version is
-        -- then supplied to external mkGrad function provided by user.
-        frozen <- U.unsafeFreeze x
-        notify frozen k
-
-        -- let grad = M.unions (map (mkGrad frozen) batch)
-        let grad = parUnions (map (mkGrad frozen) batch)
-        addUp grad u
-        scale (gain k) u
-
-        x' <- U.unsafeThaw frozen
-        apply u x'
-        doIt u (k+1) stdGen' x'
-
-
--- | Add up all gradients and store results in normal domain.
-addUp :: Grad -> MVect -> IO ()
-addUp grad v = do
-    UM.set v 0
-    forM_ (toList grad) $ \(i, x) -> do
-        y <- UM.unsafeRead v i
-        UM.unsafeWrite v i (x + y)
-
-
--- | Scale the vector by the given value.
-scale :: Double -> MVect -> IO ()
-scale c v = do
-    forM_ [0 .. UM.length v - 1] $ \i -> do
-        y <- UM.unsafeRead v i
-        UM.unsafeWrite v i (c * y)
-
-
--- | Apply gradient to the parameters vector, that is add the first vector to
--- the second one.
-apply :: MVect -> MVect -> IO ()
-apply w v = do 
-    forM_ [0 .. UM.length v - 1] $ \i -> do
-        x <- UM.unsafeRead v i
-        y <- UM.unsafeRead w i
-        UM.unsafeWrite v i (x + y)
diff --git a/Numeric/SGD/Dataset.hs b/Numeric/SGD/Dataset.hs
deleted file mode 100644
--- a/Numeric/SGD/Dataset.hs
+++ /dev/null
@@ -1,102 +0,0 @@
-{-# LANGUAGE RecordWildCards #-}
-
-
--- | Dataset abstraction.
-
-
-module Numeric.SGD.Dataset
-( 
--- * Dataset
-  Dataset (..)
--- * Reading
-, loadData
-, sample
--- * Construction
-, withVect
-, withDisk
-, withData
-) where
-
-
-import           Control.Monad (forM_)
-import           Data.Binary (Binary, encodeFile, decodeFile)
-import           System.IO.Unsafe (unsafeInterleaveIO)
-import           System.IO.Temp (withTempDirectory)
-import           System.FilePath ((</>))
-import qualified System.Random as R
-import qualified Data.Vector as V
-import qualified Control.Monad.State.Strict as S
-
-
--- | A dataset with elements of type @a@.
-data Dataset a = Dataset {
-    -- | A size of the dataset.
-      size      :: Int
-    -- | Get dataset element with a given index.  The set of indices
-    -- is of a {0, 1, .., size - 1} form.
-    , elemAt    :: Int -> IO a }
-
-
--------------------------------------------
--- Reading
--------------------------------------------
-
-
--- | Lazily load dataset from a disk.
-loadData :: Dataset a -> IO [a]
-loadData Dataset{..} = lazyMapM elemAt [0 .. size - 1]
-
-
--- | A dataset sample of the given size.
-sample :: R.RandomGen g => g -> Int -> Dataset a -> IO ([a], g)
-sample g 0 _       = return ([], g)
-sample g n dataset = do
-    (xs, g') <- sample g (n-1) dataset
-    let (i, g'') = R.next g'
-    x <- dataset `elemAt` (i `mod` size dataset)
-    return (x:xs, g'')
-
-
-lazyMapM :: (a -> IO b) -> [a] -> IO [b]
-lazyMapM f (x:xs) = do
-    y <- f x
-    ys <- unsafeInterleaveIO $ lazyMapM f xs
-    return (y:ys)
-lazyMapM _ [] = return []
-
-
--------------------------------------------
--- Construction
--------------------------------------------
-
-
--- | Construct dataset from a vector of elements and run the
--- given handler.
-withVect :: [a] -> (Dataset a -> IO b) -> IO b
-withVect xs handler =
-    handler dataset
-  where
-    v = V.fromList xs
-    dataset = Dataset
-        { size      = V.length v
-        , elemAt    = \k -> return (v V.! k) }
-
-
--- | Construct dataset from a list of elements, store it on a disk
--- and run the given handler.
-withDisk :: Binary a => [a] -> (Dataset a -> IO b) -> IO b
-withDisk xs handler = withTempDirectory "." ".sgd" $ \tmpDir -> do
-    n <- flip S.execStateT 0 $ forM_ (zip xs [0 :: Int ..]) $ \(x, ix) -> do
-        S.lift $ encodeFile (tmpDir </> show ix) x
-        S.modify (+1)
-    let at ix = decodeFile (tmpDir </> show ix)
-    handler $ Dataset {size = n, elemAt = at}
-
-
--- | Use disk or vector dataset representation depending on
--- the first argument: when `True`, use `withDisk`, otherwise
--- use `withVect`.
-withData :: Binary a => Bool -> [a] -> (Dataset a -> IO b) -> IO b
-withData x = case x of
-    True    -> withDisk
-    False   -> withVect
diff --git a/Numeric/SGD/Grad.hs b/Numeric/SGD/Grad.hs
deleted file mode 100644
--- a/Numeric/SGD/Grad.hs
+++ /dev/null
@@ -1,132 +0,0 @@
-{-# LANGUAGE CPP #-}
-
--- | A gradient is represented by an IntMap from gradient indices
--- to values. Elements with no associated values in the gradient
--- are assumed to have a 0 value assigned. Such elements are
--- not interesting: when adding the gradient to the vector of
--- parameters, only nonzero elements are taken into account.
--- 
--- Each value associated with a gradient position is a pair of
--- positive and negative components. They are stored separately
--- to ensure high accuracy of computation results.
--- Besides, both positive and negative components are stored
--- in a logarithmic domain.
-
-module Numeric.SGD.Grad
-( Grad
-, empty
-, add
-, addL
-, fromList
-, fromLogList
-, toList
-, parUnions
-) where
-
-import Data.List (foldl')
-import Control.Applicative ((<$>), (<*>))
-import Control.Monad.Par.Scheds.Direct (Par, runPar, get)
-#if MIN_VERSION_containers(0,4,2)
-import Control.Monad.Par.Scheds.Direct (spawn)
-#else
-import Control.DeepSeq (deepseq)
-import Control.Monad.Par.Scheds.Direct (spawn_)
-#endif
-import qualified Data.IntMap as M
-
-import Numeric.SGD.LogSigned
-
--- | Gradient with nonzero values stored in a logarithmic domain.
--- Since values equal to zero have no impact on the update phase
--- of the SGD method, it is more efficient to not to store those
--- components in the gradient.
-type Grad = M.IntMap LogSigned
-
-{-# INLINE insertWith #-}
-insertWith :: (a -> a -> a) -> M.Key -> a -> M.IntMap a -> M.IntMap a
-#if MIN_VERSION_containers(0,4,1)
-insertWith = M.insertWith'
-#else
-insertWith f k x m = 
-    M.alter g k m
-  where
-    g my = case my of
-        Nothing -> Just x
-        Just y  ->
-            let z = f x y
-            in  z `seq` Just z
--- insertWith f k x m = case M.lookup k m of
---     Just y  ->
---         let x' = f x y
---         in  x' `seq` M.insert k x' m
---     Nothing -> x `seq` M.insert k x m
-#endif
-
--- | Add normal-domain double to the gradient at the given position.
-{-# INLINE add #-}
-add :: Grad -> Int -> Double -> Grad
-add grad i y = insertWith (+) i (logSigned y) grad 
-
-
--- | Add log-domain, singed number to the gradient at the given position.
-{-# INLINE addL #-}
-addL :: Grad -> Int -> LogSigned -> Grad
-addL grad i y = insertWith (+) i y grad 
-
--- | Construct gradient from a list of (index, value) pairs.
--- All values from the list are added at respective gradient
--- positions.
-{-# INLINE fromList #-}
-fromList :: [(Int, Double)] -> Grad
-fromList =
-    let ins grad (i, y) = add grad i y
-    in  foldl' ins empty
-
--- | Construct gradient from a list of (index, signed, log-domain number)
--- pairs.  All values from the list are added at respective gradient
--- positions.
-{-# INLINE fromLogList #-}
-fromLogList :: [(Int, LogSigned)] -> Grad
-fromLogList =
-    let ins grad (i, y) = addL grad i y
-    in  foldl' ins empty
-
--- | Collect gradient components with values in normal domain.
-{-# INLINE toList #-}
-toList :: Grad -> [(Int, Double)]
-toList =
-    let unLog (i, x) = (i, toNorm x)
-    in  map unLog . M.assocs
-
--- | Empty gradient, i.e. with all elements set to 0.
-{-# INLINE empty #-}
-empty :: Grad
-empty = M.empty
-
--- | Perform parallel unions operation on gradient list. 
--- Experimental version.
-parUnions :: [Grad] -> Grad
-parUnions [] = error "parUnions: empty list"
-parUnions xs = runPar (parUnionsP xs)
-
--- | Parallel unions in the Par monad.
-parUnionsP :: [Grad] -> Par Grad
-parUnionsP [x] = return x
-parUnionsP zs  = do
-    let (xs, ys) = split zs
-#if MIN_VERSION_containers(0,4,2)
-    xsP <- spawn (parUnionsP xs)
-    ysP <- spawn (parUnionsP ys)
-    M.unionWith (+) <$> get xsP <*> get ysP
-#else
-    xsP <- spawn_ (parUnionsP xs)
-    ysP <- spawn_ (parUnionsP ys)
-    x <- M.unionWith (+) <$> get xsP <*> get ysP
-    M.elems x `deepseq` return x
-#endif
-  where
-    split []        = ([], [])
-    split (x:[])    = ([x], [])
-    split (x:y:rest)  =
-        let (xs, ys) = split rest
-        in  (x:xs, y:ys)
diff --git a/Numeric/SGD/LogSigned.hs b/Numeric/SGD/LogSigned.hs
deleted file mode 100644
--- a/Numeric/SGD/LogSigned.hs
+++ /dev/null
@@ -1,85 +0,0 @@
-{-# LANGUAGE BangPatterns #-}
-{-# LANGUAGE GeneralizedNewtypeDeriving #-}
-
--- | Module provides data type for signed log-domain calculations.
-
-module Numeric.SGD.LogSigned
-( LogSigned (..)
-, logSigned
-, fromPos
-, fromNeg
-, toNorm
-, toLogFloat
-) where
-
-import qualified Data.Number.LogFloat as L
-import Data.Function (on)
-import Control.DeepSeq (NFData(..))
-
--- | Signed real value in the logarithmic domain.
-data LogSigned = LogSigned
-    { pos :: {-# UNPACK #-} !L.LogFloat     -- ^ Positive component
-    , neg :: {-# UNPACK #-} !L.LogFloat     -- ^ Negative component
-    } deriving Show
-
-instance Eq LogSigned where
-    (==) = (==) `on` toLogFloat
-
-instance Ord LogSigned where
-    compare = compare `on` toLogFloat
-
--- All fields are strict and unpacked.
-instance NFData LogSigned where
-    rnf (LogSigned p q) = p `seq` q `seq` ()
-
--- | Smart LogSigned constructor.
-{-# INLINE logSigned #-}
-logSigned :: Double -> LogSigned
-logSigned x
-    | x > 0     = LogSigned (L.logFloat x) zero
-    | x < 0     = LogSigned zero (L.logFloat (-x))
-    | otherwise = LogSigned zero zero
-
--- | Make LogSigned from a positive, log-domain number.
-{-# INLINE fromPos #-}
-fromPos :: L.LogFloat -> LogSigned
-fromPos x = LogSigned x zero
-
--- | Make LogSigned from a negative, log-domain number.
-{-# INLINE fromNeg #-}
-fromNeg :: L.LogFloat -> LogSigned
-fromNeg x = LogSigned zero x
-
--- | Shift LogSigned to a normal domain.
-{-# INLINE toNorm #-}
-toNorm :: LogSigned -> Double
-toNorm (LogSigned x y) = L.fromLogFloat x - L.fromLogFloat y
-
--- | Change the 'LogSigned' to either negative 'Left' 'L.LogFloat'
--- or positive 'Right' 'L.LogFloat'.
-toLogFloat :: LogSigned -> Either L.LogFloat L.LogFloat
-toLogFloat x = case signum x of
-    -1  -> Left  $ neg x - pos x
-    1   -> Right $ pos x - neg x
-    _   -> Right $ L.logFloat (0 :: Double)
-
-instance Num LogSigned where
-    LogSigned x y + LogSigned x' y' =
-        LogSigned (x + x') (y + y')
-    LogSigned x y * LogSigned x' y' =
-        LogSigned (x*x' + y*y') (x*y' + y*x')
-    LogSigned x y - LogSigned x' y' =
-        LogSigned (x + y') (y + x')
-    negate  (LogSigned x y) = LogSigned y x
-    abs     (LogSigned x y)
-        | x >= y    = LogSigned x y
-        | otherwise = LogSigned y x
-    signum (LogSigned x y)
-        | x > y     =  1
-        | x < y     = -1
-        | otherwise =  0
-    fromInteger = logSigned . fromInteger
-
-{-# INLINE zero #-}
-zero :: L.LogFloat
-zero = L.logFloat (0 :: Double)
diff --git a/sgd.cabal b/sgd.cabal
--- a/sgd.cabal
+++ b/sgd.cabal
@@ -1,5 +1,5 @@
 name:               sgd
-version:            0.3
+version:            0.3.1
 synopsis:           Stochastic gradient descent
 description:
     Implementation of a Stochastic Gradient Descent optimization method.
@@ -21,6 +21,8 @@
 extra-source-files: examples/example1.hs
 
 library
+    hs-source-dirs: src
+
     build-depends:
         base >= 4 && < 5
       , containers
@@ -31,6 +33,7 @@
       , monad-par
       , deepseq
       , binary
+      , bytestring
       , mtl
       , filepath
       , temporary
diff --git a/src/Numeric/SGD.hs b/src/Numeric/SGD.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/SGD.hs
@@ -0,0 +1,136 @@
+{-# LANGUAGE RecordWildCards #-}
+
+
+-- | Stochastic gradient descent implementation using mutable
+-- vectors for efficient update of the parameters vector.
+-- A user is provided with the immutable vector of parameters
+-- so he is able to compute the gradient outside of the IO monad.
+-- Currently only the Gaussian priors are implemented.
+--
+-- This is a preliminary version of the SGD library and API may change
+-- in future versions.
+
+
+module Numeric.SGD
+( SgdArgs (..)
+, sgdArgsDefault
+, Para
+, sgd
+, module Numeric.SGD.Grad
+, module Numeric.SGD.Dataset
+) where
+
+
+import           Control.Monad (forM_)
+import qualified System.Random as R
+import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector.Unboxed.Mutable as UM
+import qualified Control.Monad.Primitive as Prim
+
+import           Numeric.SGD.Grad
+import           Numeric.SGD.Dataset
+
+
+-- | SGD parameters controlling the learning process.
+data SgdArgs = SgdArgs
+    { -- | Size of the batch
+      batchSize :: Int
+    -- | Regularization variance
+    , regVar    :: Double
+    -- | Number of iterations
+    , iterNum   :: Double
+    -- | Initial gain parameter
+    , gain0     :: Double
+    -- | After how many iterations over the entire dataset
+    -- the gain parameter is halved
+    , tau       :: Double }
+
+
+-- | Default SGD parameter values.
+sgdArgsDefault :: SgdArgs
+sgdArgsDefault = SgdArgs
+    { batchSize = 30
+    , regVar    = 10
+    , iterNum   = 10
+    , gain0     = 1
+    , tau       = 5 }
+
+
+-- | Vector of parameters.
+type Para       = U.Vector Double 
+
+
+-- | Type synonym for mutable vector with Double values.
+type MVect      = UM.MVector (Prim.PrimState IO) Double
+
+
+-- | A stochastic gradient descent method.
+-- A notification function can be used to provide user with
+-- information about the progress of the learning.
+sgd
+    :: SgdArgs                  -- ^ SGD parameter values
+    -> (Para -> Int -> IO ())   -- ^ Notification run every update
+    -> (Para -> x -> Grad)      -- ^ Gradient for dataset element
+    -> Dataset x                -- ^ Dataset
+    -> Para                     -- ^ Starting point
+    -> IO Para                  -- ^ SGD result
+sgd SgdArgs{..} notify mkGrad dataset x0 = do
+    u <- UM.new (U.length x0)
+    doIt u 0 (R.mkStdGen 0) =<< U.thaw x0
+  where
+    -- Gain in k-th iteration.
+    gain k = (gain0 * tau) / (tau + done k)
+
+    -- Number of completed iterations over the full dataset.
+    done k
+        = fromIntegral (k * batchSize)
+        / fromIntegral (size dataset)
+
+    doIt u k stdGen x
+      | done k > iterNum = do
+        frozen <- U.unsafeFreeze x
+        notify frozen k
+        return frozen
+      | otherwise = do
+        (batch, stdGen') <- sample stdGen batchSize dataset
+
+        -- Freeze mutable vector of parameters. The frozen version is
+        -- then supplied to external mkGrad function provided by user.
+        frozen <- U.unsafeFreeze x
+        notify frozen k
+
+        -- let grad = M.unions (map (mkGrad frozen) batch)
+        let grad = parUnions (map (mkGrad frozen) batch)
+        addUp grad u
+        scale (gain k) u
+
+        x' <- U.unsafeThaw frozen
+        apply u x'
+        doIt u (k+1) stdGen' x'
+
+
+-- | Add up all gradients and store results in normal domain.
+addUp :: Grad -> MVect -> IO ()
+addUp grad v = do
+    UM.set v 0
+    forM_ (toList grad) $ \(i, x) -> do
+        y <- UM.unsafeRead v i
+        UM.unsafeWrite v i (x + y)
+
+
+-- | Scale the vector by the given value.
+scale :: Double -> MVect -> IO ()
+scale c v = do
+    forM_ [0 .. UM.length v - 1] $ \i -> do
+        y <- UM.unsafeRead v i
+        UM.unsafeWrite v i (c * y)
+
+
+-- | Apply gradient to the parameters vector, that is add the first vector to
+-- the second one.
+apply :: MVect -> MVect -> IO ()
+apply w v = do 
+    forM_ [0 .. UM.length v - 1] $ \i -> do
+        x <- UM.unsafeRead v i
+        y <- UM.unsafeRead w i
+        UM.unsafeWrite v i (x + y)
diff --git a/src/Numeric/SGD/Dataset.hs b/src/Numeric/SGD/Dataset.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/SGD/Dataset.hs
@@ -0,0 +1,111 @@
+{-# LANGUAGE RecordWildCards #-}
+
+
+-- | Dataset abstraction.
+
+
+module Numeric.SGD.Dataset
+( 
+-- * Dataset
+  Dataset (..)
+-- * Reading
+, loadData
+, sample
+-- * Construction
+, withVect
+, withDisk
+, withData
+) where
+
+
+import           Control.Monad (forM_)
+import           Data.Binary (Binary, encodeFile, decode)
+import qualified Data.ByteString as B
+import qualified Data.ByteString.Lazy as BL
+import           System.IO.Unsafe (unsafeInterleaveIO)
+import           System.IO.Temp (withTempDirectory)
+import           System.FilePath ((</>))
+import qualified System.Random as R
+import qualified Data.Vector as V
+import qualified Control.Monad.State.Strict as S
+
+
+-- | A dataset with elements of type @a@.
+data Dataset a = Dataset {
+    -- | A size of the dataset.
+      size      :: Int
+    -- | Get dataset element with a given index.  The set of indices
+    -- is of a {0, 1, .., size - 1} form.
+    , elemAt    :: Int -> IO a }
+
+
+-------------------------------------------
+-- Reading
+-------------------------------------------
+
+
+-- | Lazily load dataset from a disk.
+loadData :: Dataset a -> IO [a]
+loadData Dataset{..} = lazyMapM elemAt [0 .. size - 1]
+
+
+-- | A dataset sample of the given size.
+sample :: R.RandomGen g => g -> Int -> Dataset a -> IO ([a], g)
+sample g 0 _       = return ([], g)
+sample g n dataset = do
+    (xs, g') <- sample g (n-1) dataset
+    let (i, g'') = R.next g'
+    x <- dataset `elemAt` (i `mod` size dataset)
+    return (x:xs, g'')
+
+
+lazyMapM :: (a -> IO b) -> [a] -> IO [b]
+lazyMapM f (x:xs) = do
+    y <- f x
+    ys <- unsafeInterleaveIO $ lazyMapM f xs
+    return (y:ys)
+lazyMapM _ [] = return []
+
+
+-------------------------------------------
+-- Construction
+-------------------------------------------
+
+
+-- | Construct dataset from a vector of elements and run the
+-- given handler.
+withVect :: [a] -> (Dataset a -> IO b) -> IO b
+withVect xs handler =
+    handler dataset
+  where
+    v = V.fromList xs
+    dataset = Dataset
+        { size      = V.length v
+        , elemAt    = \k -> return (v V.! k) }
+
+
+-- | Construct dataset from a list of elements, store it on a disk
+-- and run the given handler.
+withDisk :: Binary a => [a] -> (Dataset a -> IO b) -> IO b
+withDisk xs handler = withTempDirectory "." ".sgd" $ \tmpDir -> do
+    -- We use state monad to compute the number of dataset elements. 
+    n <- flip S.execStateT 0 $ forM_ (zip xs [0 :: Int ..]) $ \(x, ix) -> do
+        S.lift $ encodeFile (tmpDir </> show ix) x
+        S.modify (+1)
+
+    -- We need to avoid decodeFile laziness when using some older
+    -- versions of the binary library.
+    let at ix = do
+        cs <- B.readFile (tmpDir </> show ix)
+        return . decode $ BL.fromChunks [cs]
+
+    handler $ Dataset {size = n, elemAt = at}
+
+
+-- | Use disk or vector dataset representation depending on
+-- the first argument: when `True`, use `withDisk`, otherwise
+-- use `withVect`.
+withData :: Binary a => Bool -> [a] -> (Dataset a -> IO b) -> IO b
+withData x = case x of
+    True    -> withDisk
+    False   -> withVect
diff --git a/src/Numeric/SGD/Grad.hs b/src/Numeric/SGD/Grad.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/SGD/Grad.hs
@@ -0,0 +1,132 @@
+{-# LANGUAGE CPP #-}
+
+-- | A gradient is represented by an IntMap from gradient indices
+-- to values. Elements with no associated values in the gradient
+-- are assumed to have a 0 value assigned. Such elements are
+-- not interesting: when adding the gradient to the vector of
+-- parameters, only nonzero elements are taken into account.
+-- 
+-- Each value associated with a gradient position is a pair of
+-- positive and negative components. They are stored separately
+-- to ensure high accuracy of computation results.
+-- Besides, both positive and negative components are stored
+-- in a logarithmic domain.
+
+module Numeric.SGD.Grad
+( Grad
+, empty
+, add
+, addL
+, fromList
+, fromLogList
+, toList
+, parUnions
+) where
+
+import Data.List (foldl')
+import Control.Applicative ((<$>), (<*>))
+import Control.Monad.Par.Scheds.Direct (Par, runPar, get)
+#if MIN_VERSION_containers(0,4,2)
+import Control.Monad.Par.Scheds.Direct (spawn)
+#else
+import Control.DeepSeq (deepseq)
+import Control.Monad.Par.Scheds.Direct (spawn_)
+#endif
+import qualified Data.IntMap as M
+
+import Numeric.SGD.LogSigned
+
+-- | Gradient with nonzero values stored in a logarithmic domain.
+-- Since values equal to zero have no impact on the update phase
+-- of the SGD method, it is more efficient to not to store those
+-- components in the gradient.
+type Grad = M.IntMap LogSigned
+
+{-# INLINE insertWith #-}
+insertWith :: (a -> a -> a) -> M.Key -> a -> M.IntMap a -> M.IntMap a
+#if MIN_VERSION_containers(0,4,1)
+insertWith = M.insertWith'
+#else
+insertWith f k x m = 
+    M.alter g k m
+  where
+    g my = case my of
+        Nothing -> Just x
+        Just y  ->
+            let z = f x y
+            in  z `seq` Just z
+-- insertWith f k x m = case M.lookup k m of
+--     Just y  ->
+--         let x' = f x y
+--         in  x' `seq` M.insert k x' m
+--     Nothing -> x `seq` M.insert k x m
+#endif
+
+-- | Add normal-domain double to the gradient at the given position.
+{-# INLINE add #-}
+add :: Grad -> Int -> Double -> Grad
+add grad i y = insertWith (+) i (logSigned y) grad 
+
+
+-- | Add log-domain, singed number to the gradient at the given position.
+{-# INLINE addL #-}
+addL :: Grad -> Int -> LogSigned -> Grad
+addL grad i y = insertWith (+) i y grad 
+
+-- | Construct gradient from a list of (index, value) pairs.
+-- All values from the list are added at respective gradient
+-- positions.
+{-# INLINE fromList #-}
+fromList :: [(Int, Double)] -> Grad
+fromList =
+    let ins grad (i, y) = add grad i y
+    in  foldl' ins empty
+
+-- | Construct gradient from a list of (index, signed, log-domain number)
+-- pairs.  All values from the list are added at respective gradient
+-- positions.
+{-# INLINE fromLogList #-}
+fromLogList :: [(Int, LogSigned)] -> Grad
+fromLogList =
+    let ins grad (i, y) = addL grad i y
+    in  foldl' ins empty
+
+-- | Collect gradient components with values in normal domain.
+{-# INLINE toList #-}
+toList :: Grad -> [(Int, Double)]
+toList =
+    let unLog (i, x) = (i, toNorm x)
+    in  map unLog . M.assocs
+
+-- | Empty gradient, i.e. with all elements set to 0.
+{-# INLINE empty #-}
+empty :: Grad
+empty = M.empty
+
+-- | Perform parallel unions operation on gradient list. 
+-- Experimental version.
+parUnions :: [Grad] -> Grad
+parUnions [] = error "parUnions: empty list"
+parUnions xs = runPar (parUnionsP xs)
+
+-- | Parallel unions in the Par monad.
+parUnionsP :: [Grad] -> Par Grad
+parUnionsP [x] = return x
+parUnionsP zs  = do
+    let (xs, ys) = split zs
+#if MIN_VERSION_containers(0,4,2)
+    xsP <- spawn (parUnionsP xs)
+    ysP <- spawn (parUnionsP ys)
+    M.unionWith (+) <$> get xsP <*> get ysP
+#else
+    xsP <- spawn_ (parUnionsP xs)
+    ysP <- spawn_ (parUnionsP ys)
+    x <- M.unionWith (+) <$> get xsP <*> get ysP
+    M.elems x `deepseq` return x
+#endif
+  where
+    split []        = ([], [])
+    split (x:[])    = ([x], [])
+    split (x:y:rest)  =
+        let (xs, ys) = split rest
+        in  (x:xs, y:ys)
diff --git a/src/Numeric/SGD/LogSigned.hs b/src/Numeric/SGD/LogSigned.hs
new file mode 100644
--- /dev/null
+++ b/src/Numeric/SGD/LogSigned.hs
@@ -0,0 +1,85 @@
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+
+-- | Module provides data type for signed log-domain calculations.
+
+module Numeric.SGD.LogSigned
+( LogSigned (..)
+, logSigned
+, fromPos
+, fromNeg
+, toNorm
+, toLogFloat
+) where
+
+import qualified Data.Number.LogFloat as L
+import Data.Function (on)
+import Control.DeepSeq (NFData(..))
+
+-- | Signed real value in the logarithmic domain.
+data LogSigned = LogSigned
+    { pos :: {-# UNPACK #-} !L.LogFloat     -- ^ Positive component
+    , neg :: {-# UNPACK #-} !L.LogFloat     -- ^ Negative component
+    } deriving Show
+
+instance Eq LogSigned where
+    (==) = (==) `on` toLogFloat
+
+instance Ord LogSigned where
+    compare = compare `on` toLogFloat
+
+-- All fields are strict and unpacked.
+instance NFData LogSigned where
+    rnf (LogSigned p q) = p `seq` q `seq` ()
+
+-- | Smart LogSigned constructor.
+{-# INLINE logSigned #-}
+logSigned :: Double -> LogSigned
+logSigned x
+    | x > 0     = LogSigned (L.logFloat x) zero
+    | x < 0     = LogSigned zero (L.logFloat (-x))
+    | otherwise = LogSigned zero zero
+
+-- | Make LogSigned from a positive, log-domain number.
+{-# INLINE fromPos #-}
+fromPos :: L.LogFloat -> LogSigned
+fromPos x = LogSigned x zero
+
+-- | Make LogSigned from a negative, log-domain number.
+{-# INLINE fromNeg #-}
+fromNeg :: L.LogFloat -> LogSigned
+fromNeg x = LogSigned zero x
+
+-- | Shift LogSigned to a normal domain.
+{-# INLINE toNorm #-}
+toNorm :: LogSigned -> Double
+toNorm (LogSigned x y) = L.fromLogFloat x - L.fromLogFloat y
+
+-- | Change the 'LogSigned' to either negative 'Left' 'L.LogFloat'
+-- or positive 'Right' 'L.LogFloat'.
+toLogFloat :: LogSigned -> Either L.LogFloat L.LogFloat
+toLogFloat x = case signum x of
+    -1  -> Left  $ neg x - pos x
+    1   -> Right $ pos x - neg x
+    _   -> Right $ L.logFloat (0 :: Double)
+
+instance Num LogSigned where
+    LogSigned x y + LogSigned x' y' =
+        LogSigned (x + x') (y + y')
+    LogSigned x y * LogSigned x' y' =
+        LogSigned (x*x' + y*y') (x*y' + y*x')
+    LogSigned x y - LogSigned x' y' =
+        LogSigned (x + y') (y + x')
+    negate  (LogSigned x y) = LogSigned y x
+    abs     (LogSigned x y)
+        | x >= y    = LogSigned x y
+        | otherwise = LogSigned y x
+    signum (LogSigned x y)
+        | x > y     =  1
+        | x < y     = -1
+        | otherwise =  0
+    fromInteger = logSigned . fromInteger
+
+{-# INLINE zero #-}
+zero :: L.LogFloat
+zero = L.logFloat (0 :: Double)