diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -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. diff --git a/changelog b/changelog new file mode 100644 --- /dev/null +++ b/changelog @@ -0,0 +1,5 @@ +-*-change-log-*- + +0.8.0 Apr 2019 + * Allow step size decay in Adam + * Parallel calculation of the objective function diff --git a/sgd.cabal b/sgd.cabal --- a/sgd.cabal +++ b/sgd.cabal @@ -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 +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 diff --git a/src/Numeric/SGD.hs b/src/Numeric/SGD.hs --- a/src/Numeric/SGD.hs +++ b/src/Numeric/SGD.hs @@ -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) diff --git a/src/Numeric/SGD/Adam.hs b/src/Numeric/SGD/Adam.hs --- a/src/Numeric/SGD/Adam.hs +++ b/src/Numeric/SGD/Adam.hs @@ -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 diff --git a/src/Numeric/SGD/DataSet.hs b/src/Numeric/SGD/DataSet.hs --- a/src/Numeric/SGD/DataSet.hs +++ b/src/Numeric/SGD/DataSet.hs @@ -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} diff --git a/src/Numeric/SGD/Momentum.hs b/src/Numeric/SGD/Momentum.hs --- a/src/Numeric/SGD/Momentum.hs +++ b/src/Numeric/SGD/Momentum.hs @@ -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 diff --git a/src/Numeric/SGD/ParamSet.hs b/src/Numeric/SGD/ParamSet.hs --- a/src/Numeric/SGD/ParamSet.hs +++ b/src/Numeric/SGD/ParamSet.hs @@ -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 diff --git a/src/Numeric/SGD/Pipe.hs b/src/Numeric/SGD/Pipe.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/SGD/Pipe.hs @@ -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 diff --git a/src/Numeric/SGD/Sparse/Dataset.hs b/src/Numeric/SGD/Sparse/Dataset.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/SGD/Sparse/Dataset.hs @@ -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 diff --git a/test/Spec.hs b/test/Spec.hs new file mode 100644 --- /dev/null +++ b/test/Spec.hs @@ -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