ad 3.2.2 → 3.3.0.1

raw patch · 14 files changed

+737/−736 lines, 14 filesdep ~freePVP ok

version bump matches the API change (PVP)

Dependency ranges changed: free

API changes (from Hackage documentation)

Files

• CHANGELOG.markdown +6/−0
• README.markdown +1/−1
• ad.cabal +7/−7
• src/Numeric/AD/Internal/Kahn.hs +264/−0
• src/Numeric/AD/Internal/Reverse.hs +159/−195
• src/Numeric/AD/Internal/Var.hs +1/−7
• src/Numeric/AD/Internal/Wengert.hs +0/−227
• src/Numeric/AD/Mode/Directed.hs +15/−17
• src/Numeric/AD/Mode/Kahn.hs +204/−0
• src/Numeric/AD/Mode/Reverse.hs +52/−60
• src/Numeric/AD/Mode/Wengert.hs +0/−194
• src/Numeric/AD/Variadic.hs +1/−1
• src/Numeric/AD/Variadic/Kahn.hs +27/−0
• src/Numeric/AD/Variadic/Reverse.hs +0/−27

@@ -1,3 +1,9 @@+3.3+---+* Renamed `Reverse` to `Kahn` and `Wengert` to `Reverse`. We use Arthur Kahn's topological sorting algorithm to+  sort the tape after the fact in Kahn mode, while the stock Reverse mode builds a Wengert list as it goes, which+  is more efficient in practice.+ 3.2.2 ----- * Export of the `conjugateGradientDescent` and `gradientDescent` from `Numeric.AD`

@@ -87,7 +87,7 @@  * `Numeric.AD.Mode.Forward` provides basic forward-mode AD. It is good for computing simple derivatives.  * `Numeric.AD.Mode.Sparse` computes a sparse forward-mode AD tower. It is good for higher derivatives or large numbers of outputs.  * `Numeric.AD.Mode.Reverse` computes with reverse-mode AD. It is good for computing a few outputs given many inputs.- * `Numeric.AD.Mode.Chain` computes with reverse-mode AD. It is good for computing a few outputs given many inputs, when not using sparks.+ * `Numeric.AD.Mode.Wengert` computes with reverse-mode AD. It is good for computing a few outputs given many inputs, when not using sparks.  * `Numeric.AD.Mode.Tower` computes a dense forward-mode AD tower useful for higher derivatives of single input functions.   * `Numeric.AD.Newton` provides a number of combinators for root finding using Newton's method with quadratic convergence.

@@ -1,5 +1,5 @@ name:         ad-version:      3.2.2+version:      3.3.0.1 license:      BSD3 license-File: LICENSE copyright:    (c) Edward Kmett 2010-2012,@@ -24,9 +24,9 @@     .     * @Numeric.AD.Mode.Forward@ provides basic forward-mode AD. It is good for computing simple derivatives.     .-    * @Numeric.AD.Mode.Reverse@ uses benign side-effects to compute reverse-mode AD. It is good for computing gradients in one pass. It generates a tree-like tape that needs to be topologically sorted in the end.+    * @Numeric.AD.Mode.Reverse@ uses benign side-effects to compute reverse-mode AD. It is good for computing gradients in one pass. It generates a Wengert list (linear tape) using @Data.Reflection@.     .-    * @Numeric.AD.Mode.Wengert@ uses benign side-effects to compute reverse-mode AD. It is good for computing gradients in one pass. It generates a Wengert list (linear tape) using @Data.Reflection@.+    * @Numeric.AD.Mode.Kahn@ uses benign side-effects to compute reverse-mode AD. It is good for computing gradients in one pass. It generates a tree-like tape that needs to be topologically sorted in the end.     .     * @Numeric.AD.Mode.Sparse@ computes a sparse forward-mode AD tower. It is good for higher derivatives or large numbers of outputs.     .@@ -91,7 +91,7 @@     comonad          == 3.0.*,     containers       >= 0.2 && < 0.6,     data-reify       >= 0.6 && < 0.7,-    free             >= 3.0 && <= 3.3,+    free             >= 3.0 && <= 3.4,     mtl,     reflection       >= 1.1.6 && < 1.2,     tagged           >= 0.4.2.1 && < 0.5,@@ -104,24 +104,24 @@     Numeric.AD.Newton     Numeric.AD.Halley -    Numeric.AD.Mode.Wengert     Numeric.AD.Mode.Directed     Numeric.AD.Mode.Forward+    Numeric.AD.Mode.Kahn     Numeric.AD.Mode.Reverse     Numeric.AD.Mode.Tower     Numeric.AD.Mode.Sparse      Numeric.AD.Variadic-    Numeric.AD.Variadic.Reverse+    Numeric.AD.Variadic.Kahn     Numeric.AD.Variadic.Sparse      Numeric.AD.Internal.Classes     Numeric.AD.Internal.Combinators     Numeric.AD.Internal.Forward     Numeric.AD.Internal.Tower+    Numeric.AD.Internal.Kahn     Numeric.AD.Internal.Reverse     Numeric.AD.Internal.Var-    Numeric.AD.Internal.Wengert     Numeric.AD.Internal.Sparse     Numeric.AD.Internal.Dense     Numeric.AD.Internal.Composition

@@ -0,0 +1,264 @@+{-# LANGUAGE Rank2Types, TypeFamilies, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, FlexibleContexts, TemplateHaskell, UndecidableInstances, DeriveDataTypeable #-}+-- {-# OPTIONS_HADDOCK hide, prune #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.AD.Internal.Kahn+-- Copyright   :  (c) Edward Kmett 2010+-- License     :  BSD3+-- Maintainer  :  ekmett@gmail.com+-- Stability   :  experimental+-- Portability :  GHC only+--+-- This module provides reverse-mode Automatic Differentiation implementation using+-- linear time topological sorting after the fact.+--+-- For this form of reverse-mode AD we use 'System.Mem.StableName.StableName' to recover+-- sharing information from the tape to avoid combinatorial explosion, and thus+-- run asymptotically faster than it could without such sharing information, but the use+-- of side-effects contained herein is benign.+--+-----------------------------------------------------------------------------++module Numeric.AD.Internal.Kahn+    ( Kahn(..)+    , Tape(..)+    , partials+    , partialArray+    , partialMap+    , derivative+    , derivative'+    , vgrad, vgrad'+    , Grad(..)+    ) where++import Prelude hiding (mapM)+import Control.Applicative (Applicative(..),(<$>))+import Control.Monad.ST+import Control.Monad (forM_)+import Data.List (foldl')+import Data.Array.ST+import Data.Array+import Data.IntMap (IntMap, fromListWith)+import Data.Graph (Vertex, transposeG, Graph)+import Data.Reify (reifyGraph, MuRef(..))+import qualified Data.Reify.Graph as Reified+import System.IO.Unsafe (unsafePerformIO)+import Language.Haskell.TH+import Data.Data (Data)+import Data.Typeable (Typeable)+import Numeric.AD.Internal.Types+import Numeric.AD.Internal.Classes+import Numeric.AD.Internal.Identity+import Numeric.AD.Internal.Var++-- | A @Tape@ records the information needed back propagate from the output to each input during reverse 'Mode' AD.+data Tape a t+    = Zero+    | Lift !a+    | Var !a {-# UNPACK #-} !Int+    | Binary !a a a t t+    | Unary !a a t+    deriving (Show, Data, Typeable)++-- | @Kahn@ is a 'Mode' using reverse-mode automatic differentiation that provides fast 'diffFU', 'diff2FU', 'grad', 'grad2' and a fast 'jacobian' when you have a significantly smaller number of outputs than inputs.+newtype Kahn a = Kahn (Tape a (Kahn a)) deriving (Show, Typeable)++-- deriving instance (Data (Tape a (Kahn a)) => Data (Kahn a)++instance MuRef (Kahn a) where+    type DeRef (Kahn a) = Tape a++    mapDeRef _ (Kahn Zero) = pure Zero+    mapDeRef _ (Kahn (Lift a)) = pure (Lift a)+    mapDeRef _ (Kahn (Var a v)) = pure (Var a v)+    mapDeRef f (Kahn (Binary a dadb dadc b c)) = Binary a dadb dadc <$> f b <*> f c+    mapDeRef f (Kahn (Unary a dadb b)) = Unary a dadb <$> f b++instance Lifted Kahn => Mode Kahn where+    isKnownZero (Kahn Zero) = True+    isKnownZero _    = False++    isKnownConstant (Kahn Zero) = True+    isKnownConstant (Kahn (Lift _)) = True+    isKnownConstant _ = False++    auto a = Kahn (Lift a)+    zero   = Kahn Zero+    (<+>)  = binary (+) one one+    a *^ b = lift1 (a *) (\_ -> auto a) b+    a ^* b = lift1 (* b) (\_ -> auto b) a+    a ^/ b = lift1 (/ b) (\_ -> auto (recip b)) a++    Kahn Zero <**> y                = auto (0 ** primal y)+    _            <**> Kahn Zero     = auto 1+    x            <**> Kahn (Lift y) = lift1 (**y) (\z -> (y *^ z ** Id (y-1))) x+    x            <**> y                = lift2_ (**) (\z xi yi -> (yi *! z /! xi, z *! log1 xi)) x y++instance Primal Kahn where+    primal (Kahn Zero) = 0+    primal (Kahn (Lift a)) = a+    primal (Kahn (Var a _)) = a+    primal (Kahn (Binary a _ _ _ _)) = a+    primal (Kahn (Unary a _ _)) = a++instance Lifted Kahn => Jacobian Kahn where+    type D Kahn = Id++    unary f _         (Kahn Zero)     = Kahn (Lift (f 0))+    unary f _         (Kahn (Lift a)) = Kahn (Lift (f a))+    unary f (Id dadb) b                  = Kahn (Unary (f (primal b)) dadb b)++    lift1 f df b = unary f (df (Id pb)) b+        where pb = primal b++    lift1_ f df b = unary (const a) (df (Id a) (Id pb)) b+        where pb = primal b+              a = f pb++    binary f _         _         (Kahn Zero)     (Kahn Zero)     = Kahn (Lift (f 0 0))+    binary f _         _         (Kahn Zero)     (Kahn (Lift c)) = Kahn (Lift (f 0 c))+    binary f _         _         (Kahn (Lift b)) (Kahn Zero)     = Kahn (Lift (f b 0))+    binary f _         _         (Kahn (Lift b)) (Kahn (Lift c)) = Kahn (Lift (f b c))+    binary f _         (Id dadc) (Kahn Zero)     c                  = Kahn (Unary (f 0 (primal c)) dadc c)+    binary f _         (Id dadc) (Kahn (Lift b)) c                  = Kahn (Unary (f b (primal c)) dadc c)+    binary f (Id dadb) _         b                  (Kahn Zero)     = Kahn (Unary (f (primal b) 0) dadb b)+    binary f (Id dadb) _         b                  (Kahn (Lift c)) = Kahn (Unary (f (primal b) c) dadb b)+    binary f (Id dadb) (Id dadc) b                  c                  = Kahn (Binary (f (primal b) (primal c)) dadb dadc b c)++    lift2 f df b c = binary f dadb dadc b c+        where (dadb, dadc) = df (Id (primal b)) (Id (primal c))++    lift2_ f df b c = binary (\_ _ -> a) dadb dadc b c+        where+            pb = primal b+            pc = primal c+            a = f pb pc+            (dadb, dadc) = df (Id a) (Id pb) (Id pc)++deriveLifted id (conT ''Kahn)++derivative :: Num a => AD Kahn a -> a+derivative = sum . map snd . partials+{-# INLINE derivative #-}++derivative' :: Num a => AD Kahn a -> (a, a)+derivative' r = (primal r, derivative r)+{-# INLINE derivative' #-}++-- | back propagate sensitivities along a tape.+backPropagate :: Num a => (Vertex -> (Tape a Int, Int, [Int])) -> STArray s Int a -> Vertex -> ST s ()+backPropagate vmap ss v = do+        case node of+            Unary _ g b -> do+                da <- readArray ss i+                db <- readArray ss b+                writeArray ss b (db + g*da)+            Binary _ gb gc b c -> do+                da <- readArray ss i+                db <- readArray ss b+                writeArray ss b (db + gb*da)+                dc <- readArray ss c+                writeArray ss c (dc + gc*da)+            _ -> return ()+    where+        (node, i, _) = vmap v+        -- this isn't _quite_ right, as it should allow negative zeros to multiply through++topSortAcyclic :: Graph -> [Vertex]+topSortAcyclic g = reverse $ runST $ do+    del <- newArray (bounds g) False :: ST s (STUArray s Int Bool)+    let tg = transposeG g+        starters = [ n | (n, []) <- assocs tg ]+        loop [] rs = return rs+        loop (n:ns) rs = do+            writeArray del n True+            let add [] = return ns+                add (m:ms) = do+                    b <- ok (tg!m)+                    ms' <- add ms+                    if b then return (m:ms') else return ms'+                ok [] = return True+                ok (x:xs) = do b <- readArray del x; if b then ok xs else return False+            ns' <- add (g!n)+            loop ns' (n : rs)+    loop starters []++-- | This returns a list of contributions to the partials.+-- The variable ids returned in the list are likely /not/ unique!+{-# SPECIALIZE partials :: AD Kahn Double -> [(Int, Double)] #-}+partials :: forall a . Num a => AD Kahn a -> [(Int, a)]+partials (AD tape) = [ let v = sensitivities ! ix in seq v (ident, v) | (ix, Var _ ident) <- xs ]+    where+        Reified.Graph xs start = unsafePerformIO $ reifyGraph tape+        g = array xsBounds [ (i, successors t) | (i, t) <- xs ]+        vertexMap = array xsBounds xs+        vmap i = (vertexMap ! i, i, [])+        xsBounds = sbounds xs++        sensitivities = runSTArray $ do+            ss <- newArray xsBounds 0+            writeArray ss start 1+            forM_ (topSortAcyclic g) $+                backPropagate vmap ss+            return ss++        sbounds ((a,_):as) = foldl' (\(lo,hi) (b,_) -> let lo' = min lo b; hi' = max hi b in lo' `seq` hi' `seq` (lo', hi')) (a,a) as+        sbounds _ = undefined -- the graph can't be empty, it contains the output node!++        successors :: Tape a t -> [t]+        successors (Unary _ _ b) = [b]+        successors (Binary _ _ _ b c) = [b,c]+        successors _ = []++-- | Return an 'Array' of 'partials' given bounds for the variable IDs.+partialArray :: Num a => (Int, Int) -> AD Kahn a -> Array Int a+partialArray vbounds tape = accumArray (+) 0 vbounds (partials tape)+{-# INLINE partialArray #-}++-- | Return an 'IntMap' of sparse partials+partialMap :: Num a => AD Kahn a -> IntMap a+partialMap = fromListWith (+) . partials+{-# INLINE partialMap #-}++-- A simple fresh variable supply monad+newtype S a = S { runS :: Int -> (a,Int) }+instance Monad S where+    return a = S (\s -> (a,s))+    S g >>= f = S (\s -> let (a,s') = g s in runS (f a) s')++instance Var Kahn where+    var a v = Kahn (Var a v)+    varId (Kahn (Var _ v)) = v+    varId _ = error "varId: not a Var"++class Num a => Grad i o o' a | i -> a o o', o -> a i o', o' -> a i o where+    pack :: i -> [AD Kahn a] -> AD Kahn a+    unpack :: ([a] -> [a]) -> o+    unpack' :: ([a] -> (a, [a])) -> o'++instance Num a => Grad (AD Kahn a) [a] (a, [a]) a where+    pack i _ = i+    unpack f = f []+    unpack' f = f []++instance Grad i o o' a => Grad (AD Kahn a -> i) (a -> o) (a -> o') a where+    pack f (a:as) = pack (f a) as+    pack _ [] = error "Grad.pack: logic error"+    unpack f a = unpack (f . (a:))+    unpack' f a = unpack' (f . (a:))++vgrad :: Grad i o o' a => i -> o+vgrad i = unpack (unsafeGrad (pack i))+    where+        unsafeGrad f as = unbind vs (partialArray bds $ f vs)+            where+                (vs,bds) = bind as++vgrad' :: Grad i o o' a => i -> o'+vgrad' i = unpack' (unsafeGrad' (pack i))+    where+        unsafeGrad' f as = (primal r, unbind vs (partialArray bds r))+            where+                r = f vs+                (vs,bds) = bind as+

@@ -1,112 +1,137 @@-{-# LANGUAGE Rank2Types, TypeFamilies, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, FlexibleContexts, TemplateHaskell, UndecidableInstances, DeriveDataTypeable #-}+{-# LANGUAGE Rank2Types, TypeFamilies, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, FlexibleContexts, TemplateHaskell, UndecidableInstances, DeriveDataTypeable, GADTs, ScopedTypeVariables #-} -- {-# OPTIONS_HADDOCK hide, prune #-} ----------------------------------------------------------------------------- -- | -- Module      :  Numeric.AD.Internal.Reverse--- Copyright   :  (c) Edward Kmett 2010+-- Copyright   :  (c) Edward Kmett 2012 -- License     :  BSD3 -- Maintainer  :  ekmett@gmail.com -- Stability   :  experimental -- Portability :  GHC only ----- Reverse-Mode Automatic Differentiation implementation details+-- Reverse-Mode Automatic Differentiation using a single Wengert list (or \"tape\"). ----- For reverse mode AD we use 'System.Mem.StableName.StableName' to recover sharing information from--- the tape to avoid combinatorial explosion, and thus run asymptotically faster--- than it could without such sharing information, but the use of side-effects--- contained herein is benign.+-- This version uses @Data.Reflection@ to find and update the tape. --+-- This is asymptotically faster than using @Reverse@, which+-- is forced to reify and topologically sort the graph, but it requires+-- a fairly expensive rendezvous during construction when updated using+-- multiple threads.+-- -----------------------------------------------------------------------------  module Numeric.AD.Internal.Reverse     ( Reverse(..)     , Tape(..)+    , Head(..)+    , Cells(..)+    , reifyTape     , partials-    , partialArray-    , partialMap-    , derivative-    , derivative'-    , vgrad, vgrad'-    , Grad(..)+    , partialArrayOf+    , partialMapOf+    , derivativeOf+    , derivativeOf'     ) where -import Prelude hiding (mapM)-import Control.Applicative (Applicative(..),(<$>)) import Control.Monad.ST-import Control.Monad (forM_)-import Data.List (foldl') import Data.Array.ST import Data.Array-import Data.IntMap (IntMap, fromListWith)-import Data.Graph (Vertex, transposeG, Graph)-import Data.Reify (reifyGraph, MuRef(..))-import qualified Data.Reify.Graph as Reified-import System.IO.Unsafe (unsafePerformIO)-import Language.Haskell.TH-import Data.Data (Data)-import Data.Typeable (Typeable)+import Data.Array.Unsafe as Unsafe+import Data.IORef+import Data.IntMap (IntMap, fromDistinctAscList)+import Data.Proxy+import Data.Reflection+import Data.Typeable+import Language.Haskell.TH hiding (reify) import Numeric.AD.Internal.Types import Numeric.AD.Internal.Classes import Numeric.AD.Internal.Identity import Numeric.AD.Internal.Var+import Prelude hiding (mapM)+import System.IO.Unsafe (unsafePerformIO)+import Unsafe.Coerce --- | A @Tape@ records the information needed back propagate from the output to each input during 'Reverse' 'Mode' AD.-data Tape a t-    = Zero-    | Lift !a-    | Var !a {-# UNPACK #-} !Int-    | Binary !a a a t t-    | Unary !a a t-    deriving (Show, Data, Typeable)+-- evil untyped tape+data Cells where+  Nil    :: Cells+  Unary  :: {-# UNPACK #-} !Int -> a -> Cells -> Cells+  Binary :: {-# UNPACK #-} !Int -> {-# UNPACK #-} !Int -> a -> a -> Cells -> Cells --- | @Reverse@ is a 'Mode' using reverse-mode automatic differentiation that provides fast 'diffFU', 'diff2FU', 'grad', 'grad2' and a fast 'jacobian' when you have a significantly smaller number of outputs than inputs.-newtype Reverse a = Reverse (Tape a (Reverse a)) deriving (Show, Typeable)+dropCells :: Int -> Cells -> Cells+dropCells 0 xs = xs+dropCells _ Nil = Nil+dropCells n (Unary _ _ xs)      = (dropCells $! n - 1) xs+dropCells n (Binary _ _ _ _ xs) = (dropCells $! n - 1) xs --- deriving instance (Data (Tape a (Reverse a)) => Data (Reverse a)+data Head = Head {-# UNPACK #-} !Int Cells -instance MuRef (Reverse a) where-    type DeRef (Reverse a) = Tape a+newtype Tape = Tape { getTape :: IORef Head } -    mapDeRef _ (Reverse Zero) = pure Zero-    mapDeRef _ (Reverse (Lift a)) = pure (Lift a)-    mapDeRef _ (Reverse (Var a v)) = pure (Var a v)-    mapDeRef f (Reverse (Binary a dadb dadc b c)) = Binary a dadb dadc <$> f b <*> f c-    mapDeRef f (Reverse (Unary a dadb b)) = Unary a dadb <$> f b+un :: Int -> a -> Head -> (Head, Int)+un i di (Head r t) = h `seq` r' `seq` (h, r') where+  r' = r + 1+  h = Head r' (Unary i di t)+{-# INLINE un #-} -instance Lifted Reverse => Mode Reverse where-    isKnownZero (Reverse Zero) = True-    isKnownZero _    = False+bin :: Int -> Int -> a -> a -> Head -> (Head, Int)+bin i j di dj (Head r t) = h `seq` r' `seq` (h, r') where+  r' = r + 1+  h = Head r' (Binary i j di dj t)+{-# INLINE bin #-} -    isKnownConstant (Reverse Zero) = True-    isKnownConstant (Reverse (Lift _)) = True-    isKnownConstant _ = False+modifyTape :: Reifies s Tape => p s -> (Head -> (Head, r)) -> IO r+modifyTape p = atomicModifyIORef (getTape (reflect p))+{-# INLINE modifyTape #-} -    auto a = Reverse (Lift a)-    zero   = Reverse Zero-    (<+>)  = binary (+) one one-    a *^ b = lift1 (a *) (\_ -> auto a) b-    a ^* b = lift1 (* b) (\_ -> auto b) a-    a ^/ b = lift1 (/ b) (\_ -> auto (recip b)) a+-- | This is used to create a new entry on the chain given a unary function, its derivative with respect to its input,+-- the variable ID of its input, and the value of its input. Used by 'unary' and 'binary' internally.+unarily :: forall s a. Reifies s Tape => (a -> a) -> a -> Int -> a -> Reverse s a+unarily f di i b = Reverse (unsafePerformIO (modifyTape (Proxy :: Proxy s) (un i di))) $! f b+{-# INLINE unarily #-} -    Reverse Zero <**> y                = auto (0 ** primal y)-    _            <**> Reverse Zero     = auto 1-    x            <**> Reverse (Lift y) = lift1 (**y) (\z -> (y *^ z ** Id (y-1))) x-    x            <**> y                = lift2_ (**) (\z xi yi -> (yi *! z /! xi, z *! log1 xi)) x y+-- | This is used to create a new entry on the chain given a binary function, its derivatives with respect to its inputs,+-- their variable IDs and values. Used by 'binary' internally.+binarily :: forall s a. Reifies s Tape => (a -> a -> a) -> a -> a -> Int -> a -> Int -> a -> Reverse s a+binarily f di dj i b j c = Reverse (unsafePerformIO (modifyTape (Proxy :: Proxy s) (bin i j di dj))) $! f b c+{-# INLINE binarily #-} -instance Primal Reverse where-    primal (Reverse Zero) = 0-    primal (Reverse (Lift a)) = a-    primal (Reverse (Var a _)) = a-    primal (Reverse (Binary a _ _ _ _)) = a-    primal (Reverse (Unary a _ _)) = a+data Reverse s a where+  Zero :: Reverse s a+  Lift :: a -> Reverse s a+  Reverse :: {-# UNPACK #-} !Int -> a -> Reverse s a+  deriving (Show, Typeable) -instance Lifted Reverse => Jacobian Reverse where-    type D Reverse = Id+instance (Reifies s Tape, Lifted (Reverse s)) => Mode (Reverse s) where+  isKnownZero Zero = True+  isKnownZero _    = False -    unary f _         (Reverse Zero)     = Reverse (Lift (f 0))-    unary f _         (Reverse (Lift a)) = Reverse (Lift (f a))-    unary f (Id dadb) b                  = Reverse (Unary (f (primal b)) dadb b)+  isKnownConstant Reverse{} = False+  isKnownConstant _ = True +  auto = Lift+  zero = Zero+  (<+>)  = binary (+) one one+  a *^ b = lift1 (a *) (\_ -> auto a) b+  a ^* b = lift1 (* b) (\_ -> auto b) a+  a ^/ b = lift1 (/ b) (\_ -> auto (recip b)) a++  Zero <**> y      = auto (0 ** primal y)+  _    <**> Zero   = auto 1+  x    <**> Lift y = lift1 (**y) (\z -> (y *^ z ** Id (y-1))) x+  x    <**> y      = lift2_ (**) (\z xi yi -> (yi *! z /! xi, z *! log1 xi)) x y++instance Primal (Reverse s) where+    primal Zero = 0+    primal (Lift a) = a+    primal (Reverse _ a) = a++instance (Reifies s Tape, Lifted (Reverse s)) => Jacobian (Reverse s) where+    type D (Reverse s) = Id++    unary f _         (Zero)   = Lift (f 0)+    unary f _         (Lift a) = Lift (f a)+    unary f (Id dadi) (Reverse i b) = unarily f dadi i b+     lift1 f df b = unary f (df (Id pb)) b         where pb = primal b @@ -114,16 +139,17 @@         where pb = primal b               a = f pb -    binary f _         _         (Reverse Zero)     (Reverse Zero)     = Reverse (Lift (f 0 0))-    binary f _         _         (Reverse Zero)     (Reverse (Lift c)) = Reverse (Lift (f 0 c))-    binary f _         _         (Reverse (Lift b)) (Reverse Zero)     = Reverse (Lift (f b 0))-    binary f _         _         (Reverse (Lift b)) (Reverse (Lift c)) = Reverse (Lift (f b c))-    binary f _         (Id dadc) (Reverse Zero)     c                  = Reverse (Unary (f 0 (primal c)) dadc c)-    binary f _         (Id dadc) (Reverse (Lift b)) c                  = Reverse (Unary (f b (primal c)) dadc c)-    binary f (Id dadb) _         b                  (Reverse Zero)     = Reverse (Unary (f (primal b) 0) dadb b)-    binary f (Id dadb) _         b                  (Reverse (Lift c)) = Reverse (Unary (f (primal b) c) dadb b)-    binary f (Id dadb) (Id dadc) b                  c                  = Reverse (Binary (f (primal b) (primal c)) dadb dadc b c)+    binary f _         _         Zero     Zero     = Lift (f 0 0)+    binary f _         _         Zero     (Lift c) = Lift (f 0 c)+    binary f _         _         (Lift b) Zero     = Lift (f b 0)+    binary f _         _         (Lift b) (Lift c) = Lift (f b c) +    binary f _         (Id dadc) Zero        (Reverse i c) = unarily (f 0) dadc i c+    binary f _         (Id dadc) (Lift b)    (Reverse i c) = unarily (f b) dadc i c+    binary f (Id dadb) _         (Reverse i b) Zero        = unarily (`f` 0) dadb i b+    binary f (Id dadb) _         (Reverse i b) (Lift c)    = unarily (`f` c) dadb i b+    binary f (Id dadb) (Id dadc) (Reverse i b) (Reverse j c) = binarily f dadb dadc i b j c+     lift2 f df b c = binary f dadb dadc b c         where (dadb, dadc) = df (Id (primal b)) (Id (primal c)) @@ -134,130 +160,68 @@             a = f pb pc             (dadb, dadc) = df (Id a) (Id pb) (Id pc) -deriveLifted id (conT ''Reverse)--derivative :: Num a => AD Reverse a -> a-derivative = sum . map snd . partials-{-# INLINE derivative #-}--derivative' :: Num a => AD Reverse a -> (a, a)-derivative' r = (primal r, derivative r)-{-# INLINE derivative' #-}---- | back propagate sensitivities along a tape.-backPropagate :: Num a => (Vertex -> (Tape a Int, Int, [Int])) -> STArray s Int a -> Vertex -> ST s ()-backPropagate vmap ss v = do-        case node of-            Unary _ g b -> do-                da <- readArray ss i-                db <- readArray ss b-                writeArray ss b (db + g*da)-            Binary _ gb gc b c -> do-                da <- readArray ss i-                db <- readArray ss b-                writeArray ss b (db + gb*da)-                dc <- readArray ss c-                writeArray ss c (dc + gc*da)-            _ -> return ()-    where-        (node, i, _) = vmap v-        -- this isn't _quite_ right, as it should allow negative zeros to multiply through--topSortAcyclic :: Graph -> [Vertex]-topSortAcyclic g = reverse $ runST $ do-    del <- newArray (bounds g) False :: ST s (STUArray s Int Bool)-    let tg = transposeG g-        starters = [ n | (n, []) <- assocs tg ]-        loop [] rs = return rs-        loop (n:ns) rs = do-            writeArray del n True-            let add [] = return ns-                add (m:ms) = do-                    b <- ok (tg!m)-                    ms' <- add ms-                    if b then return (m:ms') else return ms'-                ok [] = return True-                ok (x:xs) = do b <- readArray del x; if b then ok xs else return False-            ns' <- add (g!n)-            loop ns' (n : rs)-    loop starters []+let s = varT (mkName "s") in+  deriveLifted (classP ''Reifies [s, conT ''Tape] :) (conT ''Reverse `appT` s) --- | This returns a list of contributions to the partials.--- The variable ids returned in the list are likely /not/ unique!-{-# SPECIALIZE partials :: AD Reverse Double -> [(Int, Double)] #-}-partials :: forall a . Num a => AD Reverse a -> [(Int, a)]-partials (AD tape) = [ let v = sensitivities ! ix in seq v (ident, v) | (ix, Var _ ident) <- xs ]-    where-        Reified.Graph xs start = unsafePerformIO $ reifyGraph tape-        g = array xsBounds [ (i, successors t) | (i, t) <- xs ]-        vertexMap = array xsBounds xs-        vmap i = (vertexMap ! i, i, [])-        xsBounds = sbounds xs+-- | Helper that extracts the derivative of a chain when the chain was constructed with one variable.+derivativeOf :: (Reifies s Tape, Num a) => Proxy s -> AD (Reverse s) a -> a+derivativeOf _ = sum . partials+{-# INLINE derivativeOf #-} -        sensitivities = runSTArray $ do-            ss <- newArray xsBounds 0-            writeArray ss start 1-            forM_ (topSortAcyclic g) $-                backPropagate vmap ss-            return ss+-- | Helper that extracts both the primal and derivative of a chain when the chain was constructed with one variable.+derivativeOf' :: (Reifies s Tape, Num a) => Proxy s -> AD (Reverse s) a -> (a, a)+derivativeOf' p r = (primal r, derivativeOf p r)+{-# INLINE derivativeOf' #-} -        sbounds ((a,_):as) = foldl' (\(lo,hi) (b,_) -> let lo' = min lo b; hi' = max hi b in lo' `seq` hi' `seq` (lo', hi')) (a,a) as-        sbounds _ = undefined -- the graph can't be empty, it contains the output node!+-- | Used internally to push sensitivities down the chain.+backPropagate :: Num a => Int -> Cells -> STArray s Int a -> ST s Int+backPropagate k Nil _ = return k+backPropagate k (Unary i g xs) ss = do+  da <- readArray ss k+  db <- readArray ss i+  writeArray ss i $! db + unsafeCoerce g*da+  (backPropagate $! k - 1) xs ss+backPropagate k (Binary i j g h xs) ss = do+  da <- readArray ss k+  db <- readArray ss i+  writeArray ss i $! db + unsafeCoerce g*da+  dc <- readArray ss j+  writeArray ss j $! dc + unsafeCoerce h*da+  (backPropagate $! k - 1) xs ss -        successors :: Tape a t -> [t]-        successors (Unary _ _ b) = [b]-        successors (Binary _ _ _ b c) = [b,c]-        successors _ = []+-- | Extract the partials from the current chain for a given AD variable.+{-# SPECIALIZE partials :: Reifies s Tape => AD (Reverse s) Double -> [Double] #-}+partials :: forall s a. (Reifies s Tape, Num a) => AD (Reverse s) a -> [a]+partials (AD Zero)        = []+partials (AD (Lift _))    = []+partials (AD (Reverse k _)) = map (sensitivities !) [0..vs] where+   Head n t = unsafePerformIO $ readIORef (getTape (reflect (Proxy :: Proxy s)))+   tk = dropCells (n - k) t+   (vs,sensitivities) = runST $ do+     ss <- newArray (0, k) 0+     writeArray ss k 1+     v <- backPropagate k tk ss+     as <- Unsafe.unsafeFreeze ss+     return (v, as)  -- | Return an 'Array' of 'partials' given bounds for the variable IDs.-partialArray :: Num a => (Int, Int) -> AD Reverse a -> Array Int a-partialArray vbounds tape = accumArray (+) 0 vbounds (partials tape)-{-# INLINE partialArray #-}+partialArrayOf :: (Reifies s Tape, Num a) => Proxy s -> (Int, Int) -> AD (Reverse s) a -> Array Int a+partialArrayOf _ vbounds = accumArray (+) 0 vbounds . zip [0..] . partials+{-# INLINE partialArrayOf #-}  -- | Return an 'IntMap' of sparse partials-partialMap :: Num a => AD Reverse a -> IntMap a-partialMap = fromListWith (+) . partials-{-# INLINE partialMap #-}+partialMapOf :: (Reifies s Tape, Num a) => Proxy s -> AD (Reverse s) a -> IntMap a+partialMapOf _ = fromDistinctAscList . zip [0..] . partials+{-# INLINE partialMapOf #-} --- A simple fresh variable supply monad-newtype S a = S { runS :: Int -> (a,Int) }-instance Monad S where-    return a = S (\s -> (a,s))-    S g >>= f = S (\s -> let (a,s') = g s in runS (f a) s')+-- | Construct a tape that starts with @n@ variables.+reifyTape :: Int -> (forall s. Reifies s Tape => Proxy s -> r) -> r+reifyTape vs k = unsafePerformIO $ do+  h <- newIORef (Head vs Nil)+  return (reify (Tape h) k)+{-# NOINLINE reifyTape #-} -instance Var Reverse where-    var a v = Reverse (Var a v)-    varId (Reverse (Var _ v)) = v+instance Var (Reverse s) where+    var a v = Reverse v a+    varId (Reverse v _) = v     varId _ = error "varId: not a Var"--class Num a => Grad i o o' a | i -> a o o', o -> a i o', o' -> a i o where-    pack :: i -> [AD Reverse a] -> AD Reverse a-    unpack :: ([a] -> [a]) -> o-    unpack' :: ([a] -> (a, [a])) -> o'--instance Num a => Grad (AD Reverse a) [a] (a, [a]) a where-    pack i _ = i-    unpack f = f []-    unpack' f = f []--instance Grad i o o' a => Grad (AD Reverse a -> i) (a -> o) (a -> o') a where-    pack f (a:as) = pack (f a) as-    pack _ [] = error "Grad.pack: logic error"-    unpack f a = unpack (f . (a:))-    unpack' f a = unpack' (f . (a:))--vgrad :: Grad i o o' a => i -> o-vgrad i = unpack (unsafeGrad (pack i))-    where-        unsafeGrad f as = unbind vs (partialArray bds $ f vs)-            where-                (vs,bds) = bind as--vgrad' :: Grad i o o' a => i -> o'-vgrad' i = unpack' (unsafeGrad' (pack i))-    where-        unsafeGrad' f as = (primal r, unbind vs (partialArray bds r))-            where-                r = f vs-                (vs,bds) = bind as-

@@ -8,13 +8,7 @@ -- Stability   :  experimental -- Portability :  GHC only ----- Reverse-Mode Automatic Differentiation using a single tape.------ This version uses @Data.Reflection@ to update a single tape.------ This is asymptotically faster than using @Reverse@, which--- is forced to reify and topologically sort the graph, but it is--- less friendly to the use of sparks.+-- Variables used for reverse-mode automatic differentiation. -----------------------------------------------------------------------------  module Numeric.AD.Internal.Var

@@ -1,227 +0,0 @@-{-# LANGUAGE Rank2Types, TypeFamilies, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, FlexibleContexts, TemplateHaskell, UndecidableInstances, DeriveDataTypeable, GADTs, ScopedTypeVariables #-}--- {-# OPTIONS_HADDOCK hide, prune #-}--------------------------------------------------------------------------------- |--- Module      :  Numeric.AD.Internal.Wengert--- Copyright   :  (c) Edward Kmett 2012--- License     :  BSD3--- Maintainer  :  ekmett@gmail.com--- Stability   :  experimental--- Portability :  GHC only------ Reverse-Mode Automatic Differentiation using a single Wengert list (or \"tape\").------ This version uses @Data.Reflection@ to find and update the tape.------ This is asymptotically faster than using @Reverse@, which--- is forced to reify and topologically sort the graph, but it requires--- a fairly expensive rendezvous during construction when updated using--- multiple threads.-----------------------------------------------------------------------------------module Numeric.AD.Internal.Wengert-    ( Wengert(..)-    , Tape(..)-    , Head(..)-    , Cells(..)-    , reifyTape-    , partials-    , partialArrayOf-    , partialMapOf-    , derivativeOf-    , derivativeOf'-    ) where--import Control.Monad.ST-import Data.Array.ST-import Data.Array-import Data.Array.Unsafe as Unsafe-import Data.IORef-import Data.IntMap (IntMap, fromDistinctAscList)-import Data.Proxy-import Data.Reflection-import Data.Typeable-import Language.Haskell.TH hiding (reify)-import Numeric.AD.Internal.Types-import Numeric.AD.Internal.Classes-import Numeric.AD.Internal.Identity-import Numeric.AD.Internal.Var-import Prelude hiding (mapM)-import System.IO.Unsafe (unsafePerformIO)-import Unsafe.Coerce---- evil untyped tape-data Cells where-  Nil    :: Cells-  Unary  :: {-# UNPACK #-} !Int -> a -> Cells -> Cells-  Binary :: {-# UNPACK #-} !Int -> {-# UNPACK #-} !Int -> a -> a -> Cells -> Cells--dropCells :: Int -> Cells -> Cells-dropCells 0 xs = xs-dropCells _ Nil = Nil-dropCells n (Unary _ _ xs)      = (dropCells $! n - 1) xs-dropCells n (Binary _ _ _ _ xs) = (dropCells $! n - 1) xs--data Head = Head {-# UNPACK #-} !Int Cells--newtype Tape = Tape { getTape :: IORef Head }--un :: Int -> a -> Head -> (Head, Int)-un i di (Head r t) = h `seq` r' `seq` (h, r') where-  r' = r + 1-  h = Head r' (Unary i di t)-{-# INLINE un #-}--bin :: Int -> Int -> a -> a -> Head -> (Head, Int)-bin i j di dj (Head r t) = h `seq` r' `seq` (h, r') where-  r' = r + 1-  h = Head r' (Binary i j di dj t)-{-# INLINE bin #-}--modifyTape :: Reifies s Tape => p s -> (Head -> (Head, r)) -> IO r-modifyTape p = atomicModifyIORef (getTape (reflect p))-{-# INLINE modifyTape #-}---- | This is used to create a new entry on the chain given a unary function, its derivative with respect to its input,--- the variable ID of its input, and the value of its input. Used by 'unary' and 'binary' internally.-unarily :: forall s a. Reifies s Tape => (a -> a) -> a -> Int -> a -> Wengert s a-unarily f di i b = Wengert (unsafePerformIO (modifyTape (Proxy :: Proxy s) (un i di))) $! f b-{-# INLINE unarily #-}---- | This is used to create a new entry on the chain given a binary function, its derivatives with respect to its inputs,--- their variable IDs and values. Used by 'binary' internally.-binarily :: forall s a. Reifies s Tape => (a -> a -> a) -> a -> a -> Int -> a -> Int -> a -> Wengert s a-binarily f di dj i b j c = Wengert (unsafePerformIO (modifyTape (Proxy :: Proxy s) (bin i j di dj))) $! f b c-{-# INLINE binarily #-}--data Wengert s a where-  Zero :: Wengert s a-  Lift :: a -> Wengert s a-  Wengert :: {-# UNPACK #-} !Int -> a -> Wengert s a-  deriving (Show, Typeable)--instance (Reifies s Tape, Lifted (Wengert s)) => Mode (Wengert s) where-  isKnownZero Zero = True-  isKnownZero _    = False--  isKnownConstant Wengert{} = False-  isKnownConstant _ = True--  auto = Lift-  zero = Zero-  (<+>)  = binary (+) one one-  a *^ b = lift1 (a *) (\_ -> auto a) b-  a ^* b = lift1 (* b) (\_ -> auto b) a-  a ^/ b = lift1 (/ b) (\_ -> auto (recip b)) a--  Zero <**> y      = auto (0 ** primal y)-  _    <**> Zero   = auto 1-  x    <**> Lift y = lift1 (**y) (\z -> (y *^ z ** Id (y-1))) x-  x    <**> y      = lift2_ (**) (\z xi yi -> (yi *! z /! xi, z *! log1 xi)) x y--instance Primal (Wengert s) where-    primal Zero = 0-    primal (Lift a) = a-    primal (Wengert _ a) = a--instance (Reifies s Tape, Lifted (Wengert s)) => Jacobian (Wengert s) where-    type D (Wengert s) = Id--    unary f _         (Zero)   = Lift (f 0)-    unary f _         (Lift a) = Lift (f a)-    unary f (Id dadi) (Wengert i b) = unarily f dadi i b--    lift1 f df b = unary f (df (Id pb)) b-        where pb = primal b--    lift1_ f df b = unary (const a) (df (Id a) (Id pb)) b-        where pb = primal b-              a = f pb--    binary f _         _         Zero     Zero     = Lift (f 0 0)-    binary f _         _         Zero     (Lift c) = Lift (f 0 c)-    binary f _         _         (Lift b) Zero     = Lift (f b 0)-    binary f _         _         (Lift b) (Lift c) = Lift (f b c)--    binary f _         (Id dadc) Zero        (Wengert i c) = unarily (f 0) dadc i c-    binary f _         (Id dadc) (Lift b)    (Wengert i c) = unarily (f b) dadc i c-    binary f (Id dadb) _         (Wengert i b) Zero        = unarily (`f` 0) dadb i b-    binary f (Id dadb) _         (Wengert i b) (Lift c)    = unarily (`f` c) dadb i b-    binary f (Id dadb) (Id dadc) (Wengert i b) (Wengert j c) = binarily f dadb dadc i b j c--    lift2 f df b c = binary f dadb dadc b c-        where (dadb, dadc) = df (Id (primal b)) (Id (primal c))--    lift2_ f df b c = binary (\_ _ -> a) dadb dadc b c-        where-            pb = primal b-            pc = primal c-            a = f pb pc-            (dadb, dadc) = df (Id a) (Id pb) (Id pc)--let s = varT (mkName "s") in-  deriveLifted (classP ''Reifies [s, conT ''Tape] :) (conT ''Wengert `appT` s)---- | Helper that extracts the derivative of a chain when the chain was constructed with one variable.-derivativeOf :: (Reifies s Tape, Num a) => Proxy s -> AD (Wengert s) a -> a-derivativeOf _ = sum . partials-{-# INLINE derivativeOf #-}---- | Helper that extracts both the primal and derivative of a chain when the chain was constructed with one variable.-derivativeOf' :: (Reifies s Tape, Num a) => Proxy s -> AD (Wengert s) a -> (a, a)-derivativeOf' p r = (primal r, derivativeOf p r)-{-# INLINE derivativeOf' #-}---- | Used internally to push sensitivities down the chain.-backPropagate :: Num a => Int -> Cells -> STArray s Int a -> ST s Int-backPropagate k Nil _ = return k-backPropagate k (Unary i g xs) ss = do-  da <- readArray ss k-  db <- readArray ss i-  writeArray ss i $! db + unsafeCoerce g*da-  (backPropagate $! k - 1) xs ss-backPropagate k (Binary i j g h xs) ss = do-  da <- readArray ss k-  db <- readArray ss i-  writeArray ss i $! db + unsafeCoerce g*da-  dc <- readArray ss j-  writeArray ss j $! dc + unsafeCoerce h*da-  (backPropagate $! k - 1) xs ss---- | Extract the partials from the current chain for a given AD variable.-{-# SPECIALIZE partials :: Reifies s Tape => AD (Wengert s) Double -> [Double] #-}-partials :: forall s a. (Reifies s Tape, Num a) => AD (Wengert s) a -> [a]-partials (AD Zero)        = []-partials (AD (Lift _))    = []-partials (AD (Wengert k _)) = map (sensitivities !) [0..vs] where-   Head n t = unsafePerformIO $ readIORef (getTape (reflect (Proxy :: Proxy s)))-   tk = dropCells (n - k) t-   (vs,sensitivities) = runST $ do-     ss <- newArray (0, k) 0-     writeArray ss k 1-     v <- backPropagate k tk ss-     as <- Unsafe.unsafeFreeze ss-     return (v, as)---- | Return an 'Array' of 'partials' given bounds for the variable IDs.-partialArrayOf :: (Reifies s Tape, Num a) => Proxy s -> (Int, Int) -> AD (Wengert s) a -> Array Int a-partialArrayOf _ vbounds = accumArray (+) 0 vbounds . zip [0..] . partials-{-# INLINE partialArrayOf #-}---- | Return an 'IntMap' of sparse partials-partialMapOf :: (Reifies s Tape, Num a) => Proxy s -> AD (Wengert s) a -> IntMap a-partialMapOf _ = fromDistinctAscList . zip [0..] . partials-{-# INLINE partialMapOf #-}---- | Construct a tape that starts with @n@ variables.-reifyTape :: Int -> (forall s. Reifies s Tape => Proxy s -> r) -> r-reifyTape vs k = unsafePerformIO $ do-  h <- newIORef (Head vs Nil)-  return (reify (Tape h) k)-{-# NOINLINE reifyTape #-}--instance Var (Wengert s) where-    var a v = Wengert v a-    varId (Wengert v _) = v-    varId _ = error "varId: not a Var"

@@ -30,67 +30,65 @@ import Prelude hiding (reverse) import Numeric.AD.Types import Data.Traversable (Traversable)-import qualified Numeric.AD.Mode.Reverse as R+import qualified Numeric.AD.Mode.Kahn as K import qualified Numeric.AD.Mode.Forward as F import qualified Numeric.AD.Mode.Tower as T-import qualified Numeric.AD.Mode.Wengert as W+import qualified Numeric.AD.Mode.Reverse as R import qualified Numeric.AD as M import Data.Ix --- TODO: use a data types a la carte approach, so we can expose more methods here--- rather than just the intersection of all of the functionality data Direction     = Forward+    | Kahn     | Reverse-    | Wengert     | Tower     | Mixed     deriving (Show, Eq, Ord, Read, Bounded, Enum, Ix)  diff :: Num a => Direction -> (forall s. Mode s => AD s a -> AD s a) -> a -> a diff Forward = F.diff+diff Kahn    = K.diff diff Reverse = R.diff-diff Wengert = W.diff-diff Tower = T.diff-diff Mixed = F.diff+diff Tower   = T.diff+diff Mixed   = F.diff {-# INLINE diff #-}  diff' :: Num a => Direction -> (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a) diff' Forward = F.diff'+diff' Kahn = K.diff' diff' Reverse = R.diff'-diff' Wengert = W.diff' diff' Tower = T.diff' diff' Mixed = F.diff' {-# INLINE diff' #-}  jacobian :: (Traversable f, Traversable g, Num a) => Direction -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f a) jacobian Forward = F.jacobian+jacobian Kahn    = K.jacobian jacobian Reverse = R.jacobian-jacobian Wengert = W.jacobian-jacobian Tower = F.jacobian -- error "jacobian Tower: unimplemented"-jacobian Mixed = M.jacobian+jacobian Tower   = F.jacobian -- error "jacobian Tower: unimplemented"+jacobian Mixed   = M.jacobian {-# INLINE jacobian #-}  jacobian' :: (Traversable f, Traversable g, Num a) => Direction -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f a) jacobian' Forward = F.jacobian'+jacobian' Kahn    = K.jacobian' jacobian' Reverse = R.jacobian'-jacobian' Wengert = W.jacobian'-jacobian' Tower = F.jacobian' -- error "jacobian' Tower: unimplemented"-jacobian' Mixed = M.jacobian'+jacobian' Tower   = F.jacobian' -- error "jacobian' Tower: unimplemented"+jacobian' Mixed   = M.jacobian' {-# INLINE jacobian' #-}  grad :: (Traversable f, Num a) => Direction -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f a grad Forward = F.grad+grad Kahn    = K.grad grad Reverse = R.grad-grad Wengert   = W.grad grad Tower   = F.grad -- error "grad Tower: unimplemented" grad Mixed   = M.grad {-# INLINE grad #-}  grad' :: (Traversable f, Num a) => Direction -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f a) grad' Forward = F.grad'+grad' Kahn    = K.grad' grad' Reverse = R.grad'-grad' Wengert   = W.grad' grad' Tower   = F.grad' -- error "grad' Tower: unimplemented" grad' Mixed   = M.grad' {-# INLINE grad' #-}

@@ -0,0 +1,204 @@+{-# LANGUAGE Rank2Types, TemplateHaskell, BangPatterns, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, UndecidableInstances, ScopedTypeVariables #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.AD.Mode.Kahn+-- Copyright   :  (c) Edward Kmett 2010+-- License     :  BSD3+-- Maintainer  :  ekmett@gmail.com+-- Stability   :  experimental+-- Portability :  GHC only+--+-- This module provides reverse-mode Automatic Differentiation using post-hoc linear time+-- topological sorting.+--+-- For reverse mode AD we use 'System.Mem.StableName.StableName' to recover sharing information from+-- the tape to avoid combinatorial explosion, and thus run asymptotically faster+-- than it could without such sharing information, but the use of side-effects+-- contained herein is benign.+--+-----------------------------------------------------------------------------++module Numeric.AD.Mode.Kahn+    (+    -- * Gradient+      grad+    , grad'+    , gradWith+    , gradWith'++    -- * Jacobian+    , jacobian+    , jacobian'+    , jacobianWith+    , jacobianWith'+    -- * Hessian+    , hessian+    , hessianF+    -- * Derivatives+    , diff+    , diff'+    , diffF+    , diffF'+    -- * Unsafe Variadic Gradient+    , vgrad, vgrad'+    , Grad+    ) where++import Control.Applicative ((<$>))+import Data.Traversable (Traversable)++import Numeric.AD.Types+import Numeric.AD.Internal.Classes+import Numeric.AD.Internal.Composition+import Numeric.AD.Internal.Kahn+import Numeric.AD.Internal.Var++-- | The 'grad' function calculates the gradient of a non-scalar-to-scalar function with reverse-mode AD in a single pass.+--+-- >>> grad (\[x,y,z] -> x*y+z) [1,2,3]+-- [2,1,1]++grad :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f a+grad f as = unbind vs (partialArray bds $ f vs)+    where (vs,bds) = bind as+{-# INLINE grad #-}++-- | The 'grad'' function calculates the result and gradient of a non-scalar-to-scalar function with reverse-mode AD in a single pass.+--+-- >>> grad' (\[x,y,z] -> 4*x*exp y+cos z) [1,2,3]+-- (28.566231899122155,[29.5562243957226,29.5562243957226,-0.1411200080598672])+grad' :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f a)+grad' f as = (primal r, unbind vs $ partialArray bds r)+    where (vs, bds) = bind as+          r = f vs+{-# INLINE grad' #-}++-- | @'grad' g f@ function calculates the gradient of a non-scalar-to-scalar function @f@ with reverse-mode AD in a single pass.+-- The gradient is combined element-wise with the argument using the function @g@.+--+-- @+-- 'grad' = 'gradWith' (\_ dx -> dx)+-- 'id' = 'gradWith' const+-- @+--+--+gradWith :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f b+gradWith g f as = unbindWith g vs (partialArray bds $ f vs)+    where (vs,bds) = bind as+{-# INLINE gradWith #-}++-- | @'grad'' g f@ calculates the result and gradient of a non-scalar-to-scalar function @f@ with reverse-mode AD in a single pass+-- the gradient is combined element-wise with the argument using the function @g@.+--+-- @'grad'' == 'gradWith'' (\_ dx -> dx)@+gradWith' :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f b)+gradWith' g f as = (primal r, unbindWith g vs $ partialArray bds r)+    where (vs, bds) = bind as+          r = f vs+{-# INLINE gradWith' #-}++-- | The 'jacobian' function calculates the jacobian of a non-scalar-to-non-scalar function with reverse AD lazily in @m@ passes for @m@ outputs.+--+-- >>> jacobian (\[x,y] -> [y,x,x*y]) [2,1]+-- [[0,1],[1,0],[1,2]]+--+-- >>> jacobian (\[x,y] -> [exp y,cos x,x+y]) [1,2]+-- [[0.0,7.38905609893065],[-0.8414709848078965,0.0],[1.0,1.0]]+jacobian :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f a)+jacobian f as = unbind vs . partialArray bds <$> f vs where+    (vs, bds) = bind as+{-# INLINE jacobian #-}++-- | The 'jacobian'' function calculates both the result and the Jacobian of a nonscalar-to-nonscalar function, using @m@ invocations of reverse AD,+-- where @m@ is the output dimensionality. Applying @fmap snd@ to the result will recover the result of 'jacobian'+-- | An alias for 'gradF''+--+-- ghci> jacobian' (\[x,y] -> [y,x,x*y]) [2,1]+-- [(1,[0,1]),(2,[1,0]),(2,[1,2])]+jacobian' :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f a)+jacobian' f as = row <$> f vs where+    (vs, bds) = bind as+    row a = (primal a, unbind vs (partialArray bds a))+{-# INLINE jacobian' #-}++-- | 'jacobianWith g f' calculates the Jacobian of a non-scalar-to-non-scalar function @f@ with reverse AD lazily in @m@ passes for @m@ outputs.+--+-- Instead of returning the Jacobian matrix, the elements of the matrix are combined with the input using the @g@.+--+-- @+-- 'jacobian' = 'jacobianWith' (\_ dx -> dx)+-- 'jacobianWith' 'const' = (\f x -> 'const' x '<$>' f x)+-- @+jacobianWith :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f b)+jacobianWith g f as = unbindWith g vs . partialArray bds <$> f vs where+    (vs, bds) = bind as+{-# INLINE jacobianWith #-}++-- | 'jacobianWith' g f' calculates both the result and the Jacobian of a nonscalar-to-nonscalar function @f@, using @m@ invocations of reverse AD,+-- where @m@ is the output dimensionality. Applying @fmap snd@ to the result will recover the result of 'jacobianWith'+--+-- Instead of returning the Jacobian matrix, the elements of the matrix are combined with the input using the @g@.+--+-- @'jacobian'' == 'jacobianWith'' (\_ dx -> dx)@+jacobianWith' :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f b)+jacobianWith' g f as = row <$> f vs where+    (vs, bds) = bind as+    row a = (primal a, unbindWith g vs (partialArray bds a))+{-# INLINE jacobianWith' #-}++-- | Compute the derivative of a function.+--+-- >>> diff sin 0+-- 1.0+--+-- >>> cos 0+-- 1.0+diff :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> a+diff f a = derivative $ f (var a 0)+{-# INLINE diff #-}++-- | The 'diff'' function calculates the value and derivative, as a+-- pair, of a scalar-to-scalar function.+--+--+-- >>> diff' sin 0+-- (0.0,1.0)+diff' :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a)+diff' f a = derivative' $ f (var a 0)+{-# INLINE diff' #-}++-- | Compute the derivatives of a function that returns a vector with regards to its single input.+--+-- >>> diffF (\a -> [sin a, cos a]) 0+-- [1.0,0.0]+diffF :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f a+diffF f a = derivative <$> f (var a 0)+{-# INLINE diffF #-}++-- | Compute the derivatives of a function that returns a vector with regards to its single input+-- as well as the primal answer.+--+-- >>> diffF' (\a -> [sin a, cos a]) 0+-- [(0.0,1.0),(1.0,0.0)]+diffF' :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f (a, a)+diffF' f a = derivative' <$> f (var a 0)+{-# INLINE diffF' #-}++-- | Compute the 'hessian' via the 'jacobian' of the gradient. gradient is computed in reverse mode and then the 'jacobian' is computed in reverse mode.+--+-- However, since the @'grad' f :: f a -> f a@ is square this is not as fast as using the forward-mode 'jacobian' of a reverse mode gradient provided by 'Numeric.AD.hessian'.+--+-- >>> hessian (\[x,y] -> x*y) [1,2]+-- [[0,1],[1,0]]+hessian :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f (f a)+hessian f = jacobian (grad (decomposeMode . f . fmap composeMode))++-- | Compute the order 3 Hessian tensor on a non-scalar-to-non-scalar function via the reverse-mode Jacobian of the reverse-mode Jacobian of the function.+--+-- Less efficient than 'Numeric.AD.Mode.Mixed.hessianF'.+--+-- >>> hessianF (\[x,y] -> [x*y,x+y,exp x*cos y]) [1,2]+-- [[[0.0,1.0],[1.0,0.0]],[[0.0,0.0],[0.0,0.0]],[[-1.1312043837568135,-2.4717266720048188],[-2.4717266720048188,1.1312043837568135]]]+hessianF :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f (f a))+hessianF f = decomposeFunctor . jacobian (ComposeFunctor . jacobian (fmap decomposeMode . f . fmap composeMode))+

@@ -8,12 +8,8 @@ -- Stability   :  experimental -- Portability :  GHC only ----- Mixed-Mode Automatic Differentiation.------ For reverse mode AD we use 'System.Mem.StableName.StableName' to recover sharing information from--- the tape to avoid combinatorial explosion, and thus run asymptotically faster--- than it could without such sharing information, but the use of side-effects--- contained herein is benign.+-- Reverse-mode automatic differentiation using Wengert lists and+-- Data.Reflection -- ----------------------------------------------------------------------------- @@ -30,17 +26,16 @@     , jacobian'     , jacobianWith     , jacobianWith'+     -- * Hessian     , hessian     , hessianF+     -- * Derivatives     , diff     , diff'     , diffF     , diffF'-    -- * Unsafe Variadic Gradient-    , vgrad, vgrad'-    , Grad     ) where  import Control.Applicative ((<$>))@@ -52,72 +47,70 @@ import Numeric.AD.Internal.Reverse import Numeric.AD.Internal.Var --- | The 'grad' function calculates the gradient of a non-scalar-to-scalar function with 'Reverse' AD in a single pass.+-- | The 'grad' function calculates the gradient of a non-scalar-to-scalar function with reverse-mode AD in a single pass. --+-- -- >>> grad (\[x,y,z] -> x*y+z) [1,2,3] -- [2,1,1]- grad :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f a-grad f as = unbind vs (partialArray bds $ f vs)-    where (vs,bds) = bind as+grad f as = reifyTape (snd bds) $ \p -> unbind vs $! partialArrayOf p bds $! f $ vary <$> vs+  where (vs, bds) = bind as {-# INLINE grad #-} --- | The 'grad'' function calculates the result and gradient of a non-scalar-to-scalar function with 'Reverse' AD in a single pass.+-- | The 'grad'' function calculates the result and gradient of a non-scalar-to-scalar function with reverse-mode AD in a single pass. ----- >>> grad' (\[x,y,z] -> 4*x*exp y+cos z) [1,2,3]--- (28.566231899122155,[29.5562243957226,29.5562243957226,-0.1411200080598672])+-- >>> grad' (\[x,y,z] -> x*y+z) [1,2,3]+-- (5,[2,1,1]) grad' :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f a)-grad' f as = (primal r, unbind vs $ partialArray bds r)-    where (vs, bds) = bind as-          r = f vs+grad' f as = reifyTape (snd bds) $ \p ->+  let r = f (fmap vary vs) in (primal r, unbind vs $! partialArrayOf p bds $! r)+  where (vs, bds) = bind as {-# INLINE grad' #-}  -- | @'grad' g f@ function calculates the gradient of a non-scalar-to-scalar function @f@ with reverse-mode AD in a single pass. -- The gradient is combined element-wise with the argument using the function @g@. -- -- @--- 'grad' = 'gradWith' (\_ dx -> dx)--- 'id' = 'gradWith' const+-- 'grad' == 'gradWith' (\_ dx -> dx)+-- 'id' == 'gradWith' 'const' -- @------ gradWith :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f b-gradWith g f as = unbindWith g vs (partialArray bds $ f vs)-    where (vs,bds) = bind as+gradWith g f as = reifyTape (snd bds) $ \p -> unbindWith g vs $! partialArrayOf p bds $! f $ vary <$> vs+  where (vs,bds) = bind as {-# INLINE gradWith #-} --- | @'grad'' g f@ calculates the result and gradient of a non-scalar-to-scalar function @f@ with 'Reverse' AD in a single pass+-- | @'grad'' g f@ calculates the result and gradient of a non-scalar-to-scalar function @f@ with reverse-mode AD in a single pass -- the gradient is combined element-wise with the argument using the function @g@. ----- @'grad'' == 'gradWith'' (\_ dx -> dx)@+-- @+-- 'grad'' == 'gradWith'' (\_ dx -> dx)+-- @ gradWith' :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f b)-gradWith' g f as = (primal r, unbindWith g vs $ partialArray bds r)+gradWith' g f as = reifyTape (snd bds) $ \p ->+   let r = f (fmap vary vs) in (primal r, unbindWith g vs $! partialArrayOf p bds $! r)     where (vs, bds) = bind as-          r = f vs {-# INLINE gradWith' #-}  -- | The 'jacobian' function calculates the jacobian of a non-scalar-to-non-scalar function with reverse AD lazily in @m@ passes for @m@ outputs. -- -- >>> jacobian (\[x,y] -> [y,x,x*y]) [2,1] -- [[0,1],[1,0],[1,2]]------ >>> jacobian (\[x,y] -> [exp y,cos x,x+y]) [1,2]--- [[0.0,7.38905609893065],[-0.8414709848078965,0.0],[1.0,1.0]] jacobian :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f a)-jacobian f as = unbind vs . partialArray bds <$> f vs where-    (vs, bds) = bind as+jacobian f as = reifyTape (snd bds) $ \p -> unbind vs . partialArrayOf p bds <$> f (fmap vary vs)+  where (vs, bds) = bind as {-# INLINE jacobian #-}  -- | The 'jacobian'' function calculates both the result and the Jacobian of a nonscalar-to-nonscalar function, using @m@ invocations of reverse AD, -- where @m@ is the output dimensionality. Applying @fmap snd@ to the result will recover the result of 'jacobian' -- | An alias for 'gradF'' ----- ghci> jacobian' (\[x,y] -> [y,x,x*y]) [2,1]+-- >>> jacobian' (\[x,y] -> [y,x,x*y]) [2,1] -- [(1,[0,1]),(2,[1,0]),(2,[1,2])] jacobian' :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f a)-jacobian' f as = row <$> f vs where-    (vs, bds) = bind as-    row a = (primal a, unbind vs (partialArray bds a))+jacobian' f as = reifyTape (snd bds) $ \p ->+  let row a = (primal a, unbind vs $! partialArrayOf p bds $! a)+  in row <$> f (vary <$> vs)+  where (vs, bds) = bind as {-# INLINE jacobian' #-}  -- | 'jacobianWith g f' calculates the Jacobian of a non-scalar-to-non-scalar function @f@ with reverse AD lazily in @m@ passes for @m@ outputs.@@ -125,11 +118,11 @@ -- Instead of returning the Jacobian matrix, the elements of the matrix are combined with the input using the @g@. -- -- @--- 'jacobian' = 'jacobianWith' (\_ dx -> dx)--- 'jacobianWith' 'const' = (\f x -> 'const' x '<$>' f x)+-- 'jacobian' == 'jacobianWith' (\_ dx -> dx)+-- 'jacobianWith' 'const' == (\f x -> 'const' x '<$>' f x) -- @ jacobianWith :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f b)-jacobianWith g f as = unbindWith g vs . partialArray bds <$> f vs where+jacobianWith g f as = reifyTape (snd bds) $ \p -> unbindWith g vs . partialArrayOf p bds <$> f (fmap vary vs) where     (vs, bds) = bind as {-# INLINE jacobianWith #-} @@ -139,53 +132,53 @@ -- Instead of returning the Jacobian matrix, the elements of the matrix are combined with the input using the @g@. -- -- @'jacobian'' == 'jacobianWith'' (\_ dx -> dx)@+-- jacobianWith' :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f b)-jacobianWith' g f as = row <$> f vs where-    (vs, bds) = bind as-    row a = (primal a, unbindWith g vs (partialArray bds a))+jacobianWith' g f as = reifyTape (snd bds) $ \p ->+  let row a = (primal a, unbindWith g vs $! partialArrayOf p bds $! a)+  in row <$> f (vary <$> vs)+  where (vs, bds) = bind as {-# INLINE jacobianWith' #-}  -- | Compute the derivative of a function. -- -- >>> diff sin 0 -- 1.0------ >>> cos 0--- 1.0 diff :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> a-diff f a = derivative $ f (var a 0)+diff f a = reifyTape 1 $ \p -> derivativeOf p $! f (var a 0) {-# INLINE diff #-} --- | The 'diff'' function calculates the value and derivative, as a--- pair, of a scalar-to-scalar function.---+-- | The 'diff'' function calculates the result and derivative, as a pair, of a scalar-to-scalar function. -- -- >>> diff' sin 0 -- (0.0,1.0)+--+-- >>> diff' exp 0+-- (1.0,1.0) diff' :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a)-diff' f a = derivative' $ f (var a 0)+diff' f a = reifyTape 1 $ \p -> derivativeOf' p $! f (var a 0) {-# INLINE diff' #-} --- | Compute the derivatives of a function that returns a vector with regards to its single input.+-- | Compute the derivatives of each result of a scalar-to-vector function with regards to its input. -- -- >>> diffF (\a -> [sin a, cos a]) 0 -- [1.0,0.0]+-- diffF :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f a-diffF f a = derivative <$> f (var a 0)+diffF f a = reifyTape 1 $ \p -> derivativeOf p <$> f (var a 0) {-# INLINE diffF #-} --- | Compute the derivatives of a function that returns a vector with regards to its single input--- as well as the primal answer.+-- | Compute the derivatives of each result of a scalar-to-vector function with regards to its input along with the answer. -- -- >>> diffF' (\a -> [sin a, cos a]) 0 -- [(0.0,1.0),(1.0,0.0)] diffF' :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f (a, a)-diffF' f a = derivative' <$> f (var a 0)+diffF' f a = reifyTape 1 $ \p -> derivativeOf' p <$> f (var a 0) {-# INLINE diffF' #-} --- | Compute the 'hessian' via the 'jacobian' of the gradient. gradient is computed in reverse mode and then the 'jacobian' is computed in reverse mode.+-- | Compute the hessian via the jacobian of the gradient. gradient is computed in reverse mode and then the jacobian is computed in reverse mode. ----- However, since the @'grad' f :: f a -> f a@ is square this is not as fast as using the forward-mode 'jacobian' of a reverse mode gradient provided by 'Numeric.AD.hessian'.+-- However, since the @'grad' f :: f a -> f a@ is square this is not as fast as using the forward-mode Jacobian of a reverse mode gradient provided by 'Numeric.AD.hessian'. -- -- >>> hessian (\[x,y] -> x*y) [1,2] -- [[0,1],[1,0]]@@ -200,4 +193,3 @@ -- [[[0.0,1.0],[1.0,0.0]],[[0.0,0.0],[0.0,0.0]],[[-1.1312043837568135,-2.4717266720048188],[-2.4717266720048188,1.1312043837568135]]] hessianF :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f (f a)) hessianF f = decomposeFunctor . jacobian (ComposeFunctor . jacobian (fmap decomposeMode . f . fmap composeMode))-

@@ -1,194 +0,0 @@-{-# LANGUAGE Rank2Types, TemplateHaskell, BangPatterns, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, UndecidableInstances, ScopedTypeVariables #-}--------------------------------------------------------------------------------- |--- Module      :  Numeric.AD.Mode.Wengert--- Copyright   :  (c) Edward Kmett 2010--- License     :  BSD3--- Maintainer  :  ekmett@gmail.com--- Stability   :  experimental--- Portability :  GHC only------ Reverse-mode automatic differentiation using Wengert lists and Data.Reflection-----------------------------------------------------------------------------------module Numeric.AD.Mode.Wengert-    (-    -- * Gradient-      grad-    , grad'-    , gradWith-    , gradWith'--    -- * Jacobian-    , jacobian-    , jacobian'-    , jacobianWith-    , jacobianWith'--    -- * Hessian-    , hessian-    , hessianF--    -- * Derivatives-    , diff-    , diff'-    , diffF-    , diffF'-    ) where--import Control.Applicative ((<$>))-import Data.Traversable (Traversable)--import Numeric.AD.Types-import Numeric.AD.Internal.Classes-import Numeric.AD.Internal.Composition-import Numeric.AD.Internal.Wengert-import Numeric.AD.Internal.Var---- | The 'grad' function calculates the gradient of a non-scalar-to-scalar function with reverse-mode AD in a single pass.--------- >>> grad (\[x,y,z] -> x*y+z) [1,2,3]--- [2,1,1]-grad :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f a-grad f as = reifyTape (snd bds) $ \p -> unbind vs $! partialArrayOf p bds $! f $ vary <$> vs-  where (vs, bds) = bind as-{-# INLINE grad #-}---- | The 'grad'' function calculates the result and gradient of a non-scalar-to-scalar function with reverse-mode AD in a single pass.------ >>> grad' (\[x,y,z] -> x*y+z) [1,2,3]--- (5,[2,1,1])-grad' :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f a)-grad' f as = reifyTape (snd bds) $ \p ->-  let r = f (fmap vary vs) in (primal r, unbind vs $! partialArrayOf p bds $! r)-  where (vs, bds) = bind as-{-# INLINE grad' #-}---- | @'grad' g f@ function calculates the gradient of a non-scalar-to-scalar function @f@ with reverse-mode AD in a single pass.--- The gradient is combined element-wise with the argument using the function @g@.------ @--- 'grad' == 'gradWith' (\_ dx -> dx)--- 'id' == 'gradWith' 'const'--- @-gradWith :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f b-gradWith g f as = reifyTape (snd bds) $ \p -> unbindWith g vs $! partialArrayOf p bds $! f $ vary <$> vs-  where (vs,bds) = bind as-{-# INLINE gradWith #-}---- | @'grad'' g f@ calculates the result and gradient of a non-scalar-to-scalar function @f@ with reverse-mode AD in a single pass--- the gradient is combined element-wise with the argument using the function @g@.------ @--- 'grad'' == 'gradWith'' (\_ dx -> dx)--- @-gradWith' :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f b)-gradWith' g f as = reifyTape (snd bds) $ \p ->-   let r = f (fmap vary vs) in (primal r, unbindWith g vs $! partialArrayOf p bds $! r)-    where (vs, bds) = bind as-{-# INLINE gradWith' #-}---- | The 'jacobian' function calculates the jacobian of a non-scalar-to-non-scalar function with reverse AD lazily in @m@ passes for @m@ outputs.------ >>> jacobian (\[x,y] -> [y,x,x*y]) [2,1]--- [[0,1],[1,0],[1,2]]-jacobian :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f a)-jacobian f as = reifyTape (snd bds) $ \p -> unbind vs . partialArrayOf p bds <$> f (fmap vary vs)-  where (vs, bds) = bind as-{-# INLINE jacobian #-}---- | The 'jacobian'' function calculates both the result and the Jacobian of a nonscalar-to-nonscalar function, using @m@ invocations of reverse AD,--- where @m@ is the output dimensionality. Applying @fmap snd@ to the result will recover the result of 'jacobian'--- | An alias for 'gradF''------ >>> jacobian' (\[x,y] -> [y,x,x*y]) [2,1]--- [(1,[0,1]),(2,[1,0]),(2,[1,2])]-jacobian' :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f a)-jacobian' f as = reifyTape (snd bds) $ \p ->-  let row a = (primal a, unbind vs $! partialArrayOf p bds $! a)-  in row <$> f (vary <$> vs)-  where (vs, bds) = bind as-{-# INLINE jacobian' #-}---- | 'jacobianWith g f' calculates the Jacobian of a non-scalar-to-non-scalar function @f@ with reverse AD lazily in @m@ passes for @m@ outputs.------ Instead of returning the Jacobian matrix, the elements of the matrix are combined with the input using the @g@.------ @--- 'jacobian' == 'jacobianWith' (\_ dx -> dx)--- 'jacobianWith' 'const' == (\f x -> 'const' x '<$>' f x)--- @-jacobianWith :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f b)-jacobianWith g f as = reifyTape (snd bds) $ \p -> unbindWith g vs . partialArrayOf p bds <$> f (fmap vary vs) where-    (vs, bds) = bind as-{-# INLINE jacobianWith #-}---- | 'jacobianWith' g f' calculates both the result and the Jacobian of a nonscalar-to-nonscalar function @f@, using @m@ invocations of reverse AD,--- where @m@ is the output dimensionality. Applying @fmap snd@ to the result will recover the result of 'jacobianWith'------ Instead of returning the Jacobian matrix, the elements of the matrix are combined with the input using the @g@.------ @'jacobian'' == 'jacobianWith'' (\_ dx -> dx)@----jacobianWith' :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f b)-jacobianWith' g f as = reifyTape (snd bds) $ \p ->-  let row a = (primal a, unbindWith g vs $! partialArrayOf p bds $! a)-  in row <$> f (vary <$> vs)-  where (vs, bds) = bind as-{-# INLINE jacobianWith' #-}---- | Compute the derivative of a function.------ >>> diff sin 0--- 1.0-diff :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> a-diff f a = reifyTape 1 $ \p -> derivativeOf p $! f (var a 0)-{-# INLINE diff #-}---- | The 'diff'' function calculates the result and derivative, as a pair, of a scalar-to-scalar function.------ >>> diff' sin 0--- (0.0,1.0)------ >>> diff' exp 0--- (1.0,1.0)-diff' :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a)-diff' f a = reifyTape 1 $ \p -> derivativeOf' p $! f (var a 0)-{-# INLINE diff' #-}---- | Compute the derivatives of each result of a scalar-to-vector function with regards to its input.------ >>> diffF (\a -> [sin a, cos a]) 0--- [1.0,0.0]----diffF :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f a-diffF f a = reifyTape 1 $ \p -> derivativeOf p <$> f (var a 0)-{-# INLINE diffF #-}---- | Compute the derivatives of each result of a scalar-to-vector function with regards to its input along with the answer.------ >>> diffF' (\a -> [sin a, cos a]) 0--- [(0.0,1.0),(1.0,0.0)]-diffF' :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f (a, a)-diffF' f a = reifyTape 1 $ \p -> derivativeOf' p <$> f (var a 0)-{-# INLINE diffF' #-}---- | Compute the hessian via the jacobian of the gradient. gradient is computed in reverse mode and then the jacobian is computed in reverse mode.------ However, since the @'grad' f :: f a -> f a@ is square this is not as fast as using the forward-mode Jacobian of a reverse mode gradient provided by 'Numeric.AD.hessian'.------ >>> hessian (\[x,y] -> x*y) [1,2]--- [[0,1],[1,0]]-hessian :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f (f a)-hessian f = jacobian (grad (decomposeMode . f . fmap composeMode))---- | Compute the order 3 Hessian tensor on a non-scalar-to-non-scalar function via the reverse-mode Jacobian of the reverse-mode Jacobian of the function.------ Less efficient than 'Numeric.AD.Mode.Mixed.hessianF'.------ >>> hessianF (\[x,y] -> [x*y,x+y,exp x*cos y]) [1,2]--- [[[0.0,1.0],[1.0,0.0]],[[0.0,0.0],[0.0,0.0]],[[-1.1312043837568135,-2.4717266720048188],[-2.4717266720048188,1.1312043837568135]]]-hessianF :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f (f a))-hessianF f = decomposeFunctor . jacobian (ComposeFunctor . jacobian (fmap decomposeMode . f . fmap composeMode))

@@ -25,5 +25,5 @@     , Grads, vgrads     ) where -import Numeric.AD.Variadic.Reverse+import Numeric.AD.Variadic.Kahn import Numeric.AD.Variadic.Sparse (Grads, vgrads)

@@ -0,0 +1,27 @@+{-# LANGUAGE Rank2Types, TemplateHaskell, BangPatterns, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, UndecidableInstances, ScopedTypeVariables #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.AD.Variadic.Kahn+-- Copyright   :  (c) Edward Kmett 2010-2012+-- License     :  BSD3+-- Maintainer  :  ekmett@gmail.com+-- Stability   :  experimental+-- Portability :  non-portable+--+-- Variadic combinators for reverse-mode automatic differentiation.+--+-- Unfortunately, variadicity comes at the expense of being able to use+-- quantification to avoid sensitivity confusion, so be careful when+-- counting the number of @lift@ you use when taking the gradient of a+-- function that takes gradients!+--+-----------------------------------------------------------------------------++module Numeric.AD.Variadic.Kahn+    (+    -- * Unsafe Variadic Gradient+      vgrad, vgrad'+    , Grad+    ) where++import Numeric.AD.Internal.Kahn

@@ -1,27 +0,0 @@-{-# LANGUAGE Rank2Types, TemplateHaskell, BangPatterns, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, UndecidableInstances, ScopedTypeVariables #-}--------------------------------------------------------------------------------- |--- Module      :  Numeric.AD.Variadic.Reverse--- Copyright   :  (c) Edward Kmett 2010-2012--- License     :  BSD3--- Maintainer  :  ekmett@gmail.com--- Stability   :  experimental--- Portability :  non-portable------ Variadic combinators for reverse-mode automatic differentiation.------ Unfortunately, variadicity comes at the expense of being able to use--- quantification to avoid sensitivity confusion, so be careful when--- counting the number of @lift@ you use when taking the gradient of a--- function that takes gradients!-----------------------------------------------------------------------------------module Numeric.AD.Variadic.Reverse-    (-    -- * Unsafe Variadic Gradient-      vgrad, vgrad'-    , Grad-    ) where--import Numeric.AD.Internal.Reverse