sgd 0.7.0.1 → 0.8.0.0

raw patch · 11 files changed

+651/−110 lines, 11 filesdep +addep +sgddep +tastydep ~basePVP ok

version bump matches the API change (PVP)

Dependencies added: ad, sgd, tasty, tasty-hunit

Dependency ranges changed: base

API changes (from Hackage documentation)

Files

• README.md +26/−0
• changelog +5/−0
• sgd.cabal +96/−54
• src/Numeric/SGD.hs +90/−41
• src/Numeric/SGD/Adam.hs +28/−5
• src/Numeric/SGD/DataSet.hs +2/−2
• src/Numeric/SGD/Momentum.hs +15/−8
• src/Numeric/SGD/ParamSet.hs +11/−0
• src/Numeric/SGD/Pipe.hs +169/−0
• src/Numeric/SGD/Sparse/Dataset.hs +124/−0
• test/Spec.hs +85/−0

@@ -0,0 +1,26 @@+# Haskell stochastic gradient descent library++Stochastic gradient descent (SGD) is a method for optimizing a global objective+function defined as a sum of smaller, differentiable functions.  In each+iteration of SGD the gradient is calculated based on a subset of the training+dataset.  In Haskell, this process can be simply represented as a [fold over a+of subsequent dataset+subsets](https://blog.jle.im/entry/purely-functional-typed-models-1.html)+(singleton elements in the extreme).++However, it can be beneficial to select the subsequent subsets randomly (e.g.,+shuffle the entire dataset before each pass).  Moreover, the dataset can be+large enough to make it impractical to store it all in memory.  Hence, the+`sgd` library adopts a [pipe](http://hackage.haskell.org/package/pipes)-based+interface in which SGD takes the form of a process consuming dataset subsets+(the so-called mini-batches) and producing a stream of output parameter values.++The `sgd` library implements several SGD variants (SGD with momentum, AdaDelta,+Adam) and handles heterogeneous parameter representations (vectors, maps, custom+records, etc.).  It can be used in combination with automatic differentiation+libraries ([ad](http://hackage.haskell.org/package/ad),+[backprop](http://hackage.haskell.org/package/backprop)), which can be used to+automatically calculate the gradient of the objective function.++Look at [the hackage repository](http://hackage.haskell.org/package/sgd) for a+library documentation.

@@ -0,0 +1,5 @@+-*-change-log-*-++0.8.0   Apr 2019+	* Allow step size decay in Adam+	* Parallel calculation of the objective function

@@ -1,59 +1,101 @@-name:               sgd-version:            0.7.0.1-synopsis:           Stochastic gradient descent-description:-    Stochastic gradient descent library.-    .-    Import "Numeric.SGD" to use the library.-license:            BSD3-license-file:       LICENSE-cabal-version:      >= 1.6-copyright:          Copyright (c) 2012-2019 Jakub Waszczuk-author:             Jakub Waszczuk-maintainer:         waszczuk.kuba@gmail.com-stability:          experimental-category:           Math-homepage:           https://github.com/kawu/sgd-build-type:         Simple+cabal-version: 1.12 -library-    hs-source-dirs: src+-- This file has been generated from package.yaml by hpack version 0.31.1.+--+-- see: https://github.com/sol/hpack+--+-- hash: c76dc64057b0ad62fdd40b4496b01d9bdbf827e48db5996639347df6c4209224 -    build-depends:-        base            >= 4        && < 5-      , containers      >= 0.4      && < 0.7-      , pipes           >= 4.3      && < 4.4-      , vector          >= 0.10     && < 0.13-      , random          >= 1.0      && < 1.2-      , random-shuffle  >= 0.0.4    && < 0.1-      , primitive       >= 0.5      && < 0.7-      , logfloat        >= 0.12     && < 0.14-      , monad-par       >= 0.3.4    && < 0.4-      , deepseq         >= 1.3      && < 1.5-      , binary          >= 0.5      && < 0.9-      , bytestring      >= 0.9      && < 0.11-      , mtl             >= 2.0      && < 2.3-      , filepath        >= 1.3      && < 1.5-      , temporary       >= 1.1      && < 1.4-      , hmatrix         >= 0.19     && < 0.20-      , data-default    >= 0.7      && < 0.8-      , parallel        >= 3.2      && < 3.3+name:           sgd+version:        0.8.0.0+synopsis:       Stochastic gradient descent library+description:    Please see the README on GitHub at <https://github.com/kawu/sgd#readme>+category:       Math+homepage:       https://github.com/kawu/sgd#readme+bug-reports:    https://github.com/kawu/sgd/issues+author:         Jakub Waszczuk+maintainer:     waszczuk.kuba@gmail.com+copyright:      2012-2019 Jakub Waszczuk+license:        BSD3+license-file:   LICENSE+build-type:     Simple+extra-source-files:+    README.md+    changelog -    exposed-modules:-        Numeric.SGD-      , Numeric.SGD.Type-      , Numeric.SGD.DataSet-      , Numeric.SGD.ParamSet-      , Numeric.SGD.Momentum-      , Numeric.SGD.AdaDelta-      , Numeric.SGD.Adam-      , Numeric.SGD.Sparse-      , Numeric.SGD.Sparse.Momentum-      , Numeric.SGD.Sparse.LogSigned-      , Numeric.SGD.Sparse.Grad+source-repository head+  type: git+  location: https://github.com/kawu/sgd -    ghc-options: -Wall -O2+library+  exposed-modules:+      Numeric.SGD+      Numeric.SGD.AdaDelta+      Numeric.SGD.Adam+      Numeric.SGD.DataSet+      Numeric.SGD.Momentum+      Numeric.SGD.ParamSet+      Numeric.SGD.Pipe+      Numeric.SGD.Sparse+      Numeric.SGD.Sparse.Dataset+      Numeric.SGD.Sparse.Grad+      Numeric.SGD.Sparse.LogSigned+      Numeric.SGD.Sparse.Momentum+      Numeric.SGD.Type+  other-modules:+      Paths_sgd+  hs-source-dirs:+      src+  build-depends:+      base >=4.7 && <5+    , binary >=0.5 && <0.9+    , bytestring >=0.9 && <0.11+    , containers >=0.4 && <0.7+    , data-default >=0.7 && <0.8+    , deepseq >=1.3 && <1.5+    , filepath >=1.3 && <1.5+    , hmatrix >=0.19 && <0.20+    , logfloat >=0.12 && <0.14+    , monad-par >=0.3.4 && <0.4+    , mtl >=2.0 && <2.3+    , parallel >=3.2 && <3.3+    , pipes >=4.3 && <4.4+    , primitive >=0.5 && <0.7+    , random >=1.0 && <1.2+    , random-shuffle >=0.0.4 && <0.1+    , temporary >=1.1 && <1.4+    , vector >=0.10 && <0.13+  default-language: Haskell2010 -source-repository head-    type: git-    location: git://github.com/kawu/sgd.git+test-suite vine-test+  type: exitcode-stdio-1.0+  main-is: Spec.hs+  other-modules:+      Paths_sgd+  hs-source-dirs:+      test+  ghc-options: -threaded -rtsopts -with-rtsopts=-N+  build-depends:+      ad >=4.3 && <4.4+    , base >=4.7 && <5+    , binary >=0.5 && <0.9+    , bytestring >=0.9 && <0.11+    , containers >=0.4 && <0.7+    , data-default >=0.7 && <0.8+    , deepseq >=1.3 && <1.5+    , filepath >=1.3 && <1.5+    , hmatrix >=0.19 && <0.20+    , logfloat >=0.12 && <0.14+    , monad-par >=0.3.4 && <0.4+    , mtl >=2.0 && <2.3+    , parallel >=3.2 && <3.3+    , pipes >=4.3 && <4.4+    , primitive >=0.5 && <0.7+    , random >=1.0 && <1.2+    , random-shuffle >=0.0.4 && <0.1+    , sgd+    , tasty >=1.2 && <1.3+    , tasty-hunit >=0.10 && <0.11+    , temporary >=1.1 && <1.4+    , vector >=0.10 && <0.13+  default-language: Haskell2010

@@ -39,6 +39,9 @@    -- * IO-based SGD   , Config (..)+  , iterNumPerEpoch+  , reportObjective+  , objectiveWith   , runIO    -- * Combinators@@ -62,20 +65,21 @@   import           GHC.Generics (Generic)--- import           GHC.Conc (numCapabilities)+import           GHC.Conc (numCapabilities)+ import           Numeric.Natural (Natural)  -- import qualified System.Random as R -import           Control.Monad (when, forM_, forever)+import           Control.Monad (when, forever) -- forM_, import           Control.Parallel.Strategies (parMap, rseq, rdeepseq, Strategy) import           Control.DeepSeq (NFData) import qualified Control.Monad.State.Strict as State  import           Data.Functor.Identity (Identity(..))-import           Data.List (foldl1') -- , transpose)+import           Data.List (foldl1', transpose) --foldl', -import qualified Data.IORef as IO+-- import qualified Data.IORef as IO import           Data.Default  import qualified Pipes as P@@ -107,9 +111,12 @@   or use an automatic differentiation library, for instance:  > import qualified Numeric.AD as AD-> derivs = map->   (\k -> AD.diff (funs !! k))->   [0..length funs-1]+> derivs =+>   [ AD.diff $ \x -> 0.3*x^2+>   , AD.diff $ \x -> -2*x+>   , AD.diff $ const 3+>   , AD.diff $ sin+>   ]    Finally, `run` allows to approach a (potentially local) minimum of the   global objective function:@@ -184,6 +191,67 @@     }  +-- | Number of new elements in each new batch+batchNew :: Config -> Int+batchNew cfg = max 1+  ( fromIntegral (batchSize cfg)+  - fromIntegral (batchOverlap cfg)+  )+++-- | Calculate the effective number of SGD iterations (and gradient+-- calculations) performed per epoch.+iterNumPerEpoch+  :: (Integral a)+  => Config+  -> a -- ^ Dataset size+  -> Double+iterNumPerEpoch cfg size =+  fromIntegral size / fromIntegral (batchNew cfg)+++-- | Report the total objective value on stdout.+reportObjective+  :: (ParamSet p)+  => (e -> p -> Double)+    -- ^ Value of the objective function on a dataset element+  -> DataSet e+    -- ^ Training dataset+  -> p -> IO Double+reportObjective objAt dataSet net = do+  q <- objectiveWith objAt dataSet net+  putStr $ show q+  putStrLn $ " (norm_2 = " ++ show (norm_2 net) ++ ")"+  return q+++-- | Value of the objective function over the entire dataset (i.e. the sum of+-- the objectives on all dataset elements).+objectiveWith+  :: (e -> p -> Double)+    -- ^ Value of the objective function on a dataset element+  -> DataSet e+    -- ^ Training dataset+  -> p -> IO Double+objectiveWith objAt dataSet net = do++--   elems <- loadData dataSet+--   let scores = parMap rseq (flip objAt net) elems+--   return $ sum scores++  parts <- partition numCapabilities <$> loadData dataSet+  let scores = parMap rseq groupScore parts+  return $ sum scores+  where+    groupScore = sum . map (flip objAt net)++--   res <- IO.newIORef 0.0+--   forM_ [0 .. size dataSet - 1] $ \ix -> do+--     x <- elemAt dataSet ix+--     IO.modifyIORef' res (+ objAt x net)+--   IO.readIORef res++ -- | Perform SGD in the IO monad, regularly reporting the value of the -- objective function on the entire dataset.  A higher-level wrapper which -- should be convenient to use when the training dataset is large.@@ -197,55 +265,36 @@     -- ^ SGD configuration   -> SGD IO [e] p     -- ^ SGD pipe consuming mini-batches of dataset elements-  -> (e -> p -> Double)-    -- ^ Value of the objective function on a dataset element (used for model-    -- quality reporting)+  -> (p -> IO Double)+    -- ^ Quality reporting function (the reporting frequency is specified+    -- via `reportEvery`)   -> DataSet e     -- ^ Training dataset   -> p     -- ^ Initial parameter values   -> IO p-runIO Config{..} sgd quality0 dataSet net0 = do-  _ <- report net0+runIO cfg@Config{..} sgd reportObj dataSet net0 = do+  _ <- reportObj net0   result net0 $ pipeData dataSet     >-> batch (fromIntegral batchSize)     >-> batchFilter     >-> sgd net0     >-> keepEvery realReportPeriod     >-> P.take (fromIntegral iterNum)-    >-> decreasingBy report+    >-> decreasingBy reportObj   where     -- Data streaming function     pipeData = forever .       if batchRandom          then pipeRan          else pipeSeq-    -- Number of new elements in each subsequent batch-    batchNew = max 1-      ( fromIntegral batchSize-      - fromIntegral batchOverlap )     -- Batch stream filter     batchFilter = do       P.await >>= P.yield-      keepEvery batchNew+      keepEvery (batchNew cfg)     -- Iteration (epoch) scaling-    iterScale x = fromIntegral (size dataSet) * x / fromIntegral batchNew-    -- Number of iterations and reporting period-    -- realIterNum = ceiling $ iterScale (fromIntegral iterNum :: Double)-    -- realReportPeriod = floor $ iterScale reportEvery-    realReportPeriod = ceiling $ iterScale reportEvery-    -- Network quality over the entire training dataset-    report net = do-      q <- quality net-      putStr $ show q-      putStrLn $ " (norm_2 = " ++ show (norm_2 net) ++ ")"-      return q-    quality net = do-      res <- IO.newIORef 0.0-      forM_ [0 .. size dataSet - 1] $ \ix -> do-        x <- elemAt dataSet ix-        IO.modifyIORef' res (+ quality0 x net)-      IO.readIORef res+    realReportPeriod = ceiling $+      reportEvery * iterNumPerEpoch cfg (size dataSet)   ------------------------------- @@ -409,9 +458,9 @@ -------------------------------  --- partition :: Int -> [a] -> [[a]]--- partition n =---     transpose . group n---   where---     group _ [] = []---     group k xs = take k xs : (group k $ drop k xs)+partition :: Int -> [a] -> [[a]]+partition n =+    transpose . group n+  where+    group _ [] = []+    group k xs = take k xs : (group k $ drop k xs)

@@ -10,6 +10,7 @@  module Numeric.SGD.Adam   ( Config(..)+  , scaleTau   , adam   ) where @@ -26,11 +27,15 @@ import           Numeric.SGD.Type import           Numeric.SGD.ParamSet +-- import Debug.Trace (trace) + -- | AdaDelta configuration data Config = Config-  { alpha :: Double-    -- ^ Step size+  { alpha0 :: Double+    -- ^ Initial step size+  , tau :: Double+    -- ^ The step size after k * `tau` iterations = `alpha0` / (k + 1)   , beta1 :: Double     -- ^ 1st exponential moment decay   , beta2 :: Double@@ -41,13 +46,20 @@  instance Default Config where   def = Config-    { alpha = 0.001+    { alpha0 = 0.001+    , tau = 10000     , beta1 = 0.9     , beta2 = 0.999     , eps = 1.0e-8     }  +-- | Scale the `tau` parameter.  Useful e.g. to account for the size of the+-- training dataset.+scaleTau :: Double -> Config -> Config+scaleTau coef cfg = cfg {tau = coef * tau cfg}++ -- | Perform gradient descent using the Adam algorithm.   -- See "Numeric.SGD.Adam" for more information. adam@@ -64,7 +76,18 @@    where +    -- Gain in the k-th iteration+    alpha k+      = (alpha0 * tau)+      / (tau + fromIntegral k)++--     report t action =+--       if t `mod` 25 == 0+--          then trace (show (tau, t, alpha t)) action+--          else action+     go t m v net = do+      -- x <- report t (P.await)       x <- P.await       let g = gradient x net           m' = pmap (*beta1) m `add` pmap (*(1-beta1)) g@@ -73,10 +96,10 @@           mb = pmap (/(1-beta1^t)) m'           vb = pmap (/(1-beta2^t)) v'           newNet = net `sub`-            ( pmap (*alpha) mb `div`+            ( pmap (*alpha t) mb `div`               (pmap (+eps) (pmap sqrt vb))             )-      P.yield newNet+      newNet `seq` P.yield newNet       go (t+1) m' v' newNet  

@@ -123,8 +123,8 @@     -- Avoid decodeFile laziness when using some older versions of the binary     -- library (as of year 2019, this could be probably simplified)     let at ix = do-        cs <- B.readFile (tmpDir </> show ix)-        return . decode $ BL.fromChunks [cs]+          cs <- B.readFile (tmpDir </> show ix)+          return . decode $ BL.fromChunks [cs]      handler $ DataSet {size = n, elemAt = at} 

@@ -10,6 +10,7 @@  module Numeric.SGD.Momentum   ( Config(..)+  , scaleTau   , momentum   ) where @@ -26,22 +27,28 @@  -- | Momentum configuration data Config = Config-  { gain0 :: Double-  -- ^ Initial gain parameter, used to scale the gradient+  { alpha0 :: Double+    -- ^ Initial step size, used to scale the gradient   , tau :: Double-  -- ^ After how many gradient calculations the gain parameter is halved+    -- ^ The step size after k * `tau` iterations = `alpha0` / (k + 1)   , gamma :: Double-  -- ^ Momentum term+    -- ^ Momentum term   } deriving (Show, Eq, Ord, Generic)  instance Default Config where   def = Config-    { gain0 = 0.01+    { alpha0 = 0.01     , gamma = 0.9     , tau = 1000     }  +-- | Scale the `tau` parameter.  Useful e.g. to account for the size of the+-- training dataset.+scaleTau :: Double -> Config -> Config+scaleTau coef cfg = cfg {tau = coef * tau cfg}++ -- | Stochastic gradient descent with momentum. See "Numeric.SGD.Momentum" for -- more information. momentum@@ -58,13 +65,13 @@   where      -- Gain in the k-th iteration-    gain k-      = (gain0 * tau)+    alpha k+      = (alpha0 * tau)       / (tau + fromIntegral k)      go k moment net = do       x <- P.await-      let grad = scale (gain k) (gradient x net)+      let grad = scale (alpha k) (gradient x net)           moment' = scale gamma moment `add` grad           newNet = net `sub` moment'       P.yield newNet

@@ -431,6 +431,17 @@   norm_2 = abs  +instance (ParamSet a, ParamSet b) => ParamSet (a, b) where+  pmap f (x, y) = (pmap f x, pmap f y)+  add (x1, y1) (x2, y2) = (x1 `add` x2, y1 `add` y2)+  sub (x1, y1) (x2, y2) = (x1 `sub` x2, y1 `sub` y2)+  mul (x1, y1) (x2, y2) = (x1 `mul` x2, y1 `mul` y2)+  div (x1, y1) (x2, y2) = (x1 `div` x2, y1 `div` y2)+  norm_2 (x, y)+    = sqrt . sum . map ((^(2::Int)))+    $ [norm_2 x, norm_2 y]++ instance (KnownNat n) => ParamSet (LA.R n) where   zero = const 0   pmap = LA.dvmap

@@ -0,0 +1,169 @@+{-# LANGUAGE RecordWildCards #-}+{-# LANGUAGE DeriveGeneric #-}+++-- | Pipe-based interface+++module Numeric.SGD.Pipe where+--   ( Config(..)+--   , Method(..)+--   , sgd+--   , result+--   , every+--   , pipeSeq+--   , pipeRan+--   ) where+++-- import           GHC.Generics (Generic)+-- import           Numeric.Natural (Natural)+-- +-- import           Control.Monad (when, forM_)+-- +-- import qualified System.Random as R+-- +-- import qualified Data.IORef as IO+-- +-- import qualified Pipes as P+-- import qualified Pipes.Prelude as P+-- import           Pipes ((>->))+-- +-- import           Numeric.SGD.DataSet+-- import           Numeric.SGD.ParamSet+-- import qualified Numeric.SGD.AdaDelta as Ada+-- import qualified Numeric.SGD.Momentum as Mom+-- +-- +-- ------------------------------- +-- -- Data+-- -------------------------------+-- +-- +-- -- | Top-level SGD configuration+-- data Config = Config+--   { iterNum :: Natural+--     -- ^ Number of iteration over the entire training dataset+--   , batchRandom :: Bool+--     -- ^ Should the mini-batch be selected at random?  If not, the subsequent+--     -- training elements will be picked sequentially.  Random selection gives+--     -- no guarantee of seeing each training sample in every epoch.  Use `False`+--     -- if unsure.+--   , method :: Method+--     -- ^ Selected SGD method+--   , reportPeriod :: Double+--     -- ^ How often the quality should be reported (with `1` meaning once per+--     -- pass over the training data)+--   } deriving (Show, Eq, Ord, Generic)+-- +-- +-- -- | Different SGD methods, together with the corresponding configurations+-- data Method+--   = AdaDelta {adaDeltaCfg :: Ada.Config}+--   | Momentum {momentumCfg :: Mom.Config}+--   deriving (Show, Eq, Ord, Generic)+-- +-- +-- ------------------------------- +-- -- SGD+-- -------------------------------+-- +-- +-- -- | Pipe-based SGD.+-- sgd+--   :: (ParamSet p)+--   => Config+--   -> DataSet e+--   -> (e -> p -> p)+--     -- ^ Network gradient on a sample element+--   -> (e -> p -> Double)+--     -- ^ Value of the objective function on a sample element+--   -> p+--     -- ^ Initial parameter values+--   -> IO p+-- sgd Config{..} dataSet grad0 quality0 net0 = do+--   let sgdPipe =+--         case method of+--           Momentum cfg -> Mom.momentum cfg+--             -- (cfg {Mom.tau = iterScale (Mom.tau cfg)})+--             grad0 net0+--           AdaDelta cfg -> Ada.adaDelta cfg grad0 net0+--   report net0+--   result net0 $ pipeSeq dataSet+--     >-> sgdPipe+--     >-> P.take realIterNum+--     >-> every realReportPeriod report+--   where+--     -- Iteration scaling+--     iterScale x = fromIntegral (size dataSet) * x+--     -- Number of iterations and reporting period+--     realIterNum = ceiling $ iterScale (fromIntegral iterNum :: Double)+--     realReportPeriod = ceiling $ iterScale reportPeriod+--     -- Network quality over the entire training dataset+--     report net = do+--       putStr . show =<< quality net+--       putStrLn $ " (norm_2 = " ++ show (norm_2 net) ++ ")"+--     quality net = do+--       res <- IO.newIORef 0.0+--       forM_ [0 .. size dataSet - 1] $ \ix -> do+--         x <- elemAt dataSet ix+--         IO.modifyIORef' res (+ quality0 x net)+--       IO.readIORef res+-- +-- +-- ------------------------------- +-- -- Dataset producers+-- -------------------------------+-- +-- +-- -- | Pipe the dataset sequentially in a loop.+-- pipeSeq :: DataSet e -> P.Producer e IO ()+-- pipeSeq dataSet = do+--   go (0 :: Int)+--   where+--     go k+--       | k >= size dataSet = go 0+--       | otherwise = do+--           x <- P.lift $ elemAt dataSet k+--           P.yield x+--           go (k+1)+-- +-- +-- -- | Pipe the dataset randomly in a loop.+-- pipeRan :: DataSet e -> P.Producer e IO ()+-- pipeRan dataSet = do+--   x <- P.lift $ do+--     ix <- R.randomRIO (0, size dataSet - 1)+--     elemAt dataSet ix+--   P.yield x+--   pipeRan dataSet+-- +-- +-- ------------------------------- +-- -- Utils+-- -------------------------------+-- +-- +-- -- | Extract the result of the SGD calculation (the last parameter+-- -- set flowing downstream).+-- result+--   :: (Monad m)+--   => p     +--     -- ^ Default value (in case the stream is empty)+--   -> P.Producer p m ()+--     -- ^ Stream of parameter sets+--   -> m p+-- result pDef = fmap (maybe pDef id) . P.last+-- +-- +-- -- | Report every `k`-th parameter set flowing downstream.+-- every :: (Monad m) => Int -> (p -> m ()) -> P.Pipe p p m x+-- every k report = do+--   go (1 `mod` k)+--   where+--     go i = do+--       paramSet <- P.await+--       when (i == 0) $ do+--         P.lift $ report paramSet+--       P.yield paramSet+--       go $ (i+1) `mod` k

@@ -0,0 +1,124 @@+{-# LANGUAGE RecordWildCards #-}+++-- | Dataset abstraction.+++module Numeric.SGD.Sparse.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.Temp (withTempDirectory)+import           System.IO.Unsafe (unsafeInterleaveIO)+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'')+++-------------------------------------------+-- 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+++-------------------------------------------+-- Lazy IO Utils+-------------------------------------------+++-- | Lazily evaluate each action in the sequence from left to right,+-- and collect the results.+lazySequence :: [IO a] -> IO [a]+lazySequence (mx:mxs) = do+    x   <- mx+    xs  <- unsafeInterleaveIO (lazySequence mxs)+    return (x : xs)+lazySequence [] = return []+++-- | `lazyMapM` f is equivalent to `lazySequence` . `map` f.+lazyMapM :: (a -> IO b) -> [a] -> IO [b]+lazyMapM f = lazySequence . map f

@@ -0,0 +1,85 @@+{-# LANGUAGE NoMonomorphismRestriction #-}+++import           Data.Ord (Ord(..))++import qualified Numeric.AD as AD++import           Test.Tasty (TestTree, testGroup)+import qualified Test.Tasty as Tasty+import qualified Test.Tasty.HUnit as U+import           Test.Tasty.HUnit ((@?=))+-- import qualified Test.Tasty.SmallCheck as SC++import qualified Numeric.SGD as SGD+import qualified Numeric.SGD.Adam as Adam+-- import qualified Numeric.SGD.AdaDelta as Ada+++main :: IO ()+main = Tasty.defaultMain tests+++tests :: TestTree+tests = testGroup "Tests" [unitTests]+++unitTests = testGroup "Unit tests"+  [ U.testCase "Simple optimization" $ do+      U.assertBool "Momentum" . approxWith 0.01 4.18 $+        sgdWith (SGD.momentum SGD.def id)+      U.assertBool "Adam" . approxWith 0.01 4.18 $+        let cfg = SGD.def {Adam.alpha0 = 0.1}+         in sgdWith (SGD.adam cfg id)+      U.assertBool "AdaDelta" . approxWith 0.25 4.18 $+         sgdWith (SGD.adaDelta SGD.def id)+--   -- the following test does not hold+--   , U.testCase "List comparison (same length)" $+--       [1, 2, 3] `compare` [1,2,2] @?= LT+  ]+++--------------------------------------------------+-- Main unit test+--------------------------------------------------+++-- | The component objective functions+funs :: [Double -> Double]+funs =+  [ \x -> 0.3*x^2+  , \x -> -2*x+  , const 3+  , sin+  ]+++-- | The corresponding derivatives+derivs :: [Double -> Double]+derivs =+  [ AD.diff $ \x -> 0.3*x^2+  , AD.diff $ \x -> -2*x+  , AD.diff $ const 3+  , AD.diff $ sin+  ]+++-- | The total objective is the sum of the objective component functions+objective :: Double -> Double+objective x = sum $ map ($x) funs+++-- | Perform SGD with the given SGD variant.+sgdWith typ = SGD.run typ (take 10000 $ cycle derivs) (0.0 :: Double)+++--------------------------------------------------+-- Utils+--------------------------------------------------+++-- | Is the second argument approximately equaly to the third one?+-- The first argument is the epsilon.+approxWith :: (Ord a, Floating a) => a -> a -> a -> Bool+approxWith eps x y =+  x >= y - eps && x <= y + eps