ad 3.1.4 → 3.2

raw patch · 18 files changed

+505/−496 lines, 18 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

• CHANGELOG.markdown +6/−0
• ad.cabal +4/−4
• src/Numeric/AD/Halley.hs +1/−1
• src/Numeric/AD/Internal/Chain.hs +0/−226
• src/Numeric/AD/Internal/Classes.hs +15/−15
• src/Numeric/AD/Internal/Composition.hs +6/−6
• src/Numeric/AD/Internal/Dense.hs +3/−3
• src/Numeric/AD/Internal/Forward.hs +10/−10
• src/Numeric/AD/Internal/Identity.hs +1/−1
• src/Numeric/AD/Internal/Reverse.hs +6/−6
• src/Numeric/AD/Internal/Sparse.hs +12/−12
• src/Numeric/AD/Internal/Tower.hs +3/−3
• src/Numeric/AD/Internal/Wengert.hs +227/−0
• src/Numeric/AD/Mode/Chain.hs +0/−194
• src/Numeric/AD/Mode/Directed.hs +8/−8
• src/Numeric/AD/Mode/Forward.hs +7/−5
• src/Numeric/AD/Mode/Wengert.hs +194/−0
• src/Numeric/AD/Newton.hs +2/−2

@@ -1,3 +1,9 @@+3.2+---+* Renamed `Chain` to `Wengert` to reflect its use of Wengert lists for reverse mode.+* Renamed `lift` to `auto` to avoid conflict with the more prevalent `transformers` library.+* Fixed a bug in `Numeric.AD.Forward.gradWith'`, which caused it to return the wrong value for the primal.+ 3.1.4 ----- * Added a better "convergence" test for `findZero`

@@ -1,5 +1,5 @@ name:         ad-version:      3.1.4+version:      3.2 license:      BSD3 license-File: LICENSE copyright:    (c) Edward Kmett 2010-2012,@@ -26,7 +26,7 @@     .     * @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.Chain@ uses benign side-effects to compute reverse-mode AD. It is good for computing gradients in one pass. It generates a 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.Sparse@ computes a sparse forward-mode AD tower. It is good for higher derivatives or large numbers of outputs.     .@@ -104,7 +104,7 @@     Numeric.AD.Newton     Numeric.AD.Halley -    Numeric.AD.Mode.Chain+    Numeric.AD.Mode.Wengert     Numeric.AD.Mode.Directed     Numeric.AD.Mode.Forward     Numeric.AD.Mode.Reverse@@ -121,7 +121,7 @@     Numeric.AD.Internal.Tower     Numeric.AD.Internal.Reverse     Numeric.AD.Internal.Var-    Numeric.AD.Internal.Chain+    Numeric.AD.Internal.Wengert     Numeric.AD.Internal.Sparse     Numeric.AD.Internal.Dense     Numeric.AD.Internal.Composition

@@ -58,7 +58,7 @@ -- Note: the @take 10 $ inverse sqrt 1 (sqrt 10)@ example that works for Newton's method -- fails with Halley's method because the preconditions do not hold! inverse :: (Fractional a, Eq a) => (forall s. Mode s => AD s a -> AD s a) -> a -> a -> [a]-inverse f x0 y = findZero (\x -> f x - lift y) x0+inverse f x0 y = findZero (\x -> f x - auto y) x0 {-# INLINE inverse  #-}  -- | The 'fixedPoint' function find a fixedpoint of a scalar

@@ -1,226 +0,0 @@-{-# LANGUAGE Rank2Types, TypeFamilies, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, FlexibleContexts, TemplateHaskell, UndecidableInstances, DeriveDataTypeable, GADTs, ScopedTypeVariables #-}--- {-# OPTIONS_HADDOCK hide, prune #-}--------------------------------------------------------------------------------- |--- Module      :  Numeric.AD.Internal.Chain--- Copyright   :  (c) Edward Kmett 2012--- License     :  BSD3--- Maintainer  :  ekmett@gmail.com--- Stability   :  experimental--- Portability :  GHC only------ Reverse-Mode Automatic Differentiation using a single 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.-----------------------------------------------------------------------------------module Numeric.AD.Internal.Chain-    ( Chain(..)-    , 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 -> Chain s a-unarily f di i b = Chain (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 -> Chain s a-binarily f di dj i b j c = Chain (unsafePerformIO (modifyTape (Proxy :: Proxy s) (bin i j di dj))) $! f b c-{-# INLINE binarily #-}--data Chain s a where-  Zero :: Chain s a-  Lift :: a -> Chain s a-  Chain :: {-# UNPACK #-} !Int -> a -> Chain s a-  deriving (Show, Typeable)--instance (Reifies s Tape, Lifted (Chain s)) => Mode (Chain s) where-  isKnownZero Zero = True-  isKnownZero _    = False--  isKnownConstant Chain{} = False-  isKnownConstant _ = True--  lift = Lift-  zero = Zero-  (<+>)  = binary (+) one one-  a *^ b = lift1 (a *) (\_ -> lift a) b-  a ^* b = lift1 (* b) (\_ -> lift b) a-  a ^/ b = lift1 (/ b) (\_ -> lift (recip b)) a--  Zero <**> y      = lift (0 ** primal y)-  _    <**> Zero   = lift 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 (Chain s) where-    primal Zero = 0-    primal (Lift a) = a-    primal (Chain _ a) = a--instance (Reifies s Tape, Lifted (Chain s)) => Jacobian (Chain s) where-    type D (Chain s) = Id--    unary f _         (Zero)   = Lift (f 0)-    unary f _         (Lift a) = Lift (f a)-    unary f (Id dadi) (Chain 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        (Chain i c) = unarily (f 0) dadc i c-    binary f _         (Id dadc) (Lift b)    (Chain i c) = unarily (f b) dadc i c-    binary f (Id dadb) _         (Chain i b) Zero        = unarily (`f` 0) dadb i b-    binary f (Id dadb) _         (Chain i b) (Lift c)    = unarily (`f` c) dadb i b-    binary f (Id dadb) (Id dadc) (Chain i b) (Chain 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 ''Chain `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 (Chain 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 (Chain 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 (Chain s) Double -> [Double] #-}-partials :: forall s a. (Reifies s Tape, Num a) => AD (Chain s) a -> [a]-partials (AD Zero)        = []-partials (AD (Lift _))    = []-partials (AD (Chain 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 (Chain 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 (Chain 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 (Chain s) where-    var a v = Chain v a-    varId (Chain v _) = v-    varId _ = error "varId: not a Var"

@@ -92,7 +92,7 @@     isKnownZero _ = False      -- | Embed a constant-    lift  :: Num a => a -> t a+    auto  :: Num a => a -> t a      -- | Vector sum     (<+>) :: Num a => t a -> t a -> t a@@ -113,19 +113,19 @@     -- | > 'zero' = 'lift' 0     zero :: Num a => t a -    a *^ b = lift a *! b-    a ^* b = a *! lift b+    a *^ b = auto a *! b+    a ^* b = a *! auto b      a ^/ b = a ^* recip b -    zero = lift 0+    zero = auto 0  one :: (Mode t, Num a) => t a-one = lift 1+one = auto 1 {-# INLINE one #-}  negOne :: (Mode t, Num a) => t a-negOne = lift (-1)+negOne = auto (-1) {-# INLINE negOne #-}  -- | 'Primal' is used by 'deriveMode' but is not exposed@@ -133,7 +133,7 @@ -- via the AD data type. -- -- It provides direct access to the result, stripped of its derivative information,--- but this is unsafe in general as (lift . primal) would discard derivative+-- but this is unsafe in general as (auto . primal) would discard derivative -- information. The end user is protected from accidentally using this function -- by the universal quantification on the various combinators we expose. @@ -200,11 +200,11 @@        instance Lifted $_t where         (==!)         = (==) `on` primal         compare1      = compare `on` primal-        maxBound1     = lift maxBound-        minBound1     = lift minBound+        maxBound1     = auto maxBound+        minBound1     = auto minBound         showsPrec1 d  = showsPrec d . primal         fromInteger1 0 = zero-        fromInteger1 n = lift (fromInteger n)+        fromInteger1 n = auto (fromInteger n)         (+!)          = (<+>) -- binary (+) one one         (-!)          = binary (-) one negOne -- TODO: <-> ? as it is, this might be pretty bad for Tower         (*!)          = lift2 (*) (\x y -> (y, x))@@ -212,14 +212,14 @@         abs1          = lift1 abs signum1         signum1       = lift1 signum (const zero)         fromRational1 0 = zero-        fromRational1 r = lift (fromRational r)+        fromRational1 r = auto (fromRational r)         x /! y        = x *! recip1 y         recip1        = lift1_ recip (const . negate1 . square1)-        pi1       = lift pi+        pi1       = auto pi         exp1      = lift1_ exp const         log1      = lift1 log recip1         logBase1 x y = log1 y /! log1 x-        sqrt1     = lift1_ sqrt (\z _ -> recip1 (lift 2 *! z))+        sqrt1     = lift1_ sqrt (\z _ -> recip1 (auto 2 *! z))         (**!)     = (<**>)         --x **! y         --   | isKnownZero y     = 1@@ -240,7 +240,7 @@          succ1                 = lift1 succ (const one)         pred1                 = lift1 pred (const one)-        toEnum1               = lift . toEnum+        toEnum1               = auto . toEnum         fromEnum1             = discrete1 fromEnum         enumFrom1 a           = withPrimal a <$> discrete1 enumFrom a         enumFromTo1 a b       = withPrimal a <$> discrete2 enumFromTo a b@@ -252,7 +252,7 @@         floatDigits1     = discrete1 floatDigits         floatRange1      = discrete1 floatRange         decodeFloat1     = discrete1 decodeFloat-        encodeFloat1 m e = lift (encodeFloat m e)+        encodeFloat1 m e = auto (encodeFloat m e)         isNaN1           = discrete1 isNaN         isInfinite1      = discrete1 isInfinite         isDenormalized1  = discrete1 isDenormalized

@@ -93,18 +93,18 @@     primal = primal . primal . runComposeMode  instance (Mode f, Mode g) => Mode (ComposeMode f g) where-    lift = ComposeMode . lift . lift+    auto = ComposeMode . auto . auto     ComposeMode a <+> ComposeMode b = ComposeMode (a <+> b)-    a *^ ComposeMode b = ComposeMode (lift a *^ b)-    ComposeMode a ^* b = ComposeMode (a ^* lift b)-    ComposeMode a ^/ b = ComposeMode (a ^/ lift b)+    a *^ ComposeMode b = ComposeMode (auto a *^ b)+    ComposeMode a ^* b = ComposeMode (a ^* auto b)+    ComposeMode a ^/ b = ComposeMode (a ^/ auto b)     ComposeMode a <**> ComposeMode b = ComposeMode (a <**> b)  instance (Mode f, Mode g) => Lifted (ComposeMode f g) where     showsPrec1 n (ComposeMode a) = showsPrec1 n a     ComposeMode a ==! ComposeMode b  = a ==! b     compare1 (ComposeMode a) (ComposeMode b) = compare1 a b-    fromInteger1 = ComposeMode . lift . fromInteger1+    fromInteger1 = ComposeMode . auto . fromInteger1     ComposeMode a +! ComposeMode b = ComposeMode (a +! b)     ComposeMode a -! ComposeMode b = ComposeMode (a -! b)     ComposeMode a *! ComposeMode b = ComposeMode (a *! b)@@ -113,7 +113,7 @@     signum1 (ComposeMode a) = ComposeMode (signum1 a)     ComposeMode a /! ComposeMode b = ComposeMode (a /! b)     recip1 (ComposeMode a) = ComposeMode (recip1 a)-    fromRational1 = ComposeMode . lift . fromRational1+    fromRational1 = ComposeMode . auto . fromRational1     toRational1 (ComposeMode a) = toRational1 a     pi1 = ComposeMode pi1     exp1 (ComposeMode a) = ComposeMode (exp1 a)

@@ -80,7 +80,7 @@     primal (Dense a _) = a  instance (Traversable f, Lifted (Dense f)) => Mode (Dense f) where-    lift = Lift+    auto = Lift     zero = Zero      Zero <+> a = a@@ -90,8 +90,8 @@     Dense a da <+> Lift b     = Dense (a + b) da     Dense a da <+> Dense b db = Dense (a + b) $ zipWithT (+) da db -    Zero <**> y      = lift (0 ** primal y)-    _    <**> Zero   = lift 1+    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 

@@ -37,7 +37,7 @@ import Numeric.AD.Internal.Classes import Numeric.AD.Internal.Identity --- | 'Forward' mode AD.+-- | 'Forward' mode AD data Forward a   = Forward !a a   | Lift !a@@ -70,7 +70,7 @@     primal Zero = 0  instance Lifted Forward => Mode Forward where-    lift = Lift+    auto = Lift     zero = Zero      isKnownZero Zero = True@@ -86,8 +86,8 @@     Lift a <+> Forward b db = Forward (a + b) db     Lift a <+> Lift b = Lift (a + b) -    Zero <**> y      = lift (0 ** primal y)-    _    <**> Zero   = lift 1+    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 @@ -168,28 +168,28 @@ bind f as = snd $ mapAccumL outer (0 :: Int) as     where         outer !i _ = (i + 1, f $ snd $ mapAccumL (inner i) 0 as)-        inner !i !j a = (j + 1, if i == j then bundle a 1 else lift a)+        inner !i !j a = (j + 1, if i == j then bundle a 1 else auto a)  bind' :: (Traversable f, Num a) => (f (AD Forward a) -> b) -> f a -> (b, f b) bind' f as = dropIx $ mapAccumL outer (0 :: Int, b0) as     where         outer (!i, _) _ = let b = f $ snd $ mapAccumL (inner i) (0 :: Int) as in ((i + 1, b), b)-        inner !i !j a = (j + 1, if i == j then bundle a 1 else lift a)-        b0 = f (lift <$> as)+        inner !i !j a = (j + 1, if i == j then bundle a 1 else auto a)+        b0 = f (auto <$> as)         dropIx ((_,b),bs) = (b,bs)  bindWith :: (Traversable f, Num a) => (a -> b -> c) -> (f (AD Forward a) -> b) -> f a -> f c bindWith g f as = snd $ mapAccumL outer (0 :: Int) as     where         outer !i a = (i + 1, g a $ f $ snd $ mapAccumL (inner i) 0 as)-        inner !i !j a = (j + 1, if i == j then bundle a 1 else AD Zero)+        inner !i !j a = (j + 1, if i == j then bundle a 1 else auto a)  bindWith' :: (Traversable f, Num a) => (a -> b -> c) -> (f (AD Forward a) -> b) -> f a -> (b, f c) bindWith' g f as = dropIx $ mapAccumL outer (0 :: Int, b0) as     where         outer (!i, _) a = let b = f $ snd $ mapAccumL (inner i) (0 :: Int) as in ((i + 1, b), g a b)-        inner !i !j a = (j + 1, if i == j then bundle a 1 else AD Zero)-        b0 = f (lift <$> as)+        inner !i !j a = (j + 1, if i == j then bundle a 1 else auto a)+        b0 = f (auto <$> as)         dropIx ((_,b),bs) = (b,bs)  -- we can't transpose arbitrary traversables, since we can't construct one out of whole cloth, and the outer

@@ -129,7 +129,7 @@     maxBound1 = maxBound  instance Mode Id where-    lift = Id+    auto = Id     Id a ^* b = Id (a * b)     a *^ Id b = Id (a * b)     Id a <+> Id b = Id (a + b)

@@ -81,15 +81,15 @@     isKnownConstant (Reverse (Lift _)) = True     isKnownConstant _ = False -    lift a = Reverse (Lift a)+    auto a = Reverse (Lift a)     zero   = Reverse Zero     (<+>)  = binary (+) one one-    a *^ b = lift1 (a *) (\_ -> lift a) b-    a ^* b = lift1 (* b) (\_ -> lift b) a-    a ^/ b = lift1 (/ b) (\_ -> lift (recip b)) a+    a *^ b = lift1 (a *) (\_ -> auto a) b+    a ^* b = lift1 (* b) (\_ -> auto b) a+    a ^/ b = lift1 (/ b) (\_ -> auto (recip b)) a -    Reverse Zero <**> y                = lift (0 ** primal y)-    _            <**> Reverse Zero     = lift 1+    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 

@@ -72,7 +72,7 @@ vars :: (Traversable f, Num a) => f a -> f (AD Sparse a) vars = snd . mapAccumL var 0     where-        var !n a = (n + 1, AD $ Sparse a $ singleton n $ lift 1)+        var !n a = (n + 1, AD $ Sparse a $ singleton n $ auto 1) {-# INLINE vars #-}  apply :: (Traversable f, Num a) => (f (AD Sparse a) -> b) -> f a -> b@@ -105,7 +105,7 @@  {- vvars :: Num a => Vector a -> Vector (AD Sparse a)-vvars = Vector.imap (\n a -> AD $ Sparse a $ singleton n $ lift 1)+vvars = Vector.imap (\n a -> AD $ Sparse a $ singleton n $ auto 1) {-# INLINE vvars #-}  vapply :: Num a => (Vector (AD Sparse a) -> b) -> Vector a -> b@@ -131,7 +131,7 @@  partial :: Num a => [Int] -> Sparse a -> a partial []     (Sparse a _)  = a-partial (n:ns) (Sparse _ da) = partial ns $ findWithDefault (lift 0) n da+partial (n:ns) (Sparse _ da) = partial ns $ findWithDefault (auto 0) n da partial _      Zero          = 0 {-# INLINE partial #-} @@ -148,10 +148,10 @@     primal Zero = 0  instance Lifted Sparse => Mode Sparse where-    lift a = Sparse a IntMap.empty+    auto a = Sparse a IntMap.empty     zero = Zero-    Zero <**> y    = lift (0 ** primal y)-    _    <**> Zero = lift 1+    Zero <**> y    = auto (0 ** primal y)+    _    <**> Zero = auto 1     x    <**> y@(Sparse b bs)       | IntMap.null bs = lift1 (**b) (\z -> (b *^ z <**> Sparse (b-1) IntMap.empty)) x       | otherwise      = lift2_ (**) (\z xi yi -> (yi *! z /! xi, z *! log1 xi)) x y@@ -167,17 +167,17 @@  instance Lifted Sparse => Jacobian Sparse where     type D Sparse = Sparse-    unary f _ Zero = lift (f 0)+    unary f _ Zero = auto (f 0)     unary f dadb (Sparse pb bs) = Sparse (f pb) $ mapWithKey (times dadb) bs -    lift1 f _ Zero = lift (f 0)+    lift1 f _ Zero = auto (f 0)     lift1 f df b@(Sparse pb bs) = Sparse (f pb) $ mapWithKey (times (df b)) bs -    lift1_ f _  Zero = lift (f 0)+    lift1_ f _  Zero = auto (f 0)     lift1_ f df b@(Sparse pb bs) = a where         a = Sparse (f pb) $ mapWithKey (times (df a b)) bs -    binary f _    _    Zero           Zero           = lift (f 0 0)+    binary f _    _    Zero           Zero           = auto (f 0 0)     binary f _    dadc Zero           (Sparse pc dc) = Sparse (f 0  pc) $ mapWithKey (times dadc) dc     binary f dadb _    (Sparse pb db) Zero           = Sparse (f pb 0 ) $ mapWithKey (times dadb) db     binary f dadb dadc (Sparse pb db) (Sparse pc dc) = Sparse (f pb pc) $@@ -185,7 +185,7 @@             (mapWithKey (times dadb) db)             (mapWithKey (times dadc) dc) -    lift2 f _  Zero             Zero = lift (f 0 0)+    lift2 f _  Zero             Zero = auto (f 0 0)     lift2 f df Zero c@(Sparse pc dc) = Sparse (f 0 pc) $ mapWithKey (times dadc) dc where dadc = snd (df zero c)     lift2 f df b@(Sparse pb db) Zero = Sparse (f pb 0) $ mapWithKey (times dadb) db where dadb = fst (df b zero)     lift2 f df b@(Sparse pb db) c@(Sparse pc dc) = Sparse (f pb pc) da where@@ -194,7 +194,7 @@             (mapWithKey (times dadb) db)             (mapWithKey (times dadc) dc) -    lift2_ f _  Zero             Zero = lift (f 0 0)+    lift2_ f _  Zero             Zero = auto (f 0 0)     lift2_ f df b@(Sparse pb db) Zero = a where a = Sparse (f pb 0) (mapWithKey (times (fst (df a b zero))) db)     lift2_ f df Zero c@(Sparse pc dc) = a where a = Sparse (f 0 pc) (mapWithKey (times (snd (df a zero c))) dc)     lift2_ f df b@(Sparse pb db) c@(Sparse pc dc) = a where

@@ -103,10 +103,10 @@     primal _ = 0  instance Lifted Tower => Mode Tower where-    lift a = Tower [a]+    auto a = Tower [a]     zero = Tower []-    Tower [] <**> y         = lift (0 ** primal y)-    _        <**> Tower []  = lift 1+    Tower [] <**> y         = auto (0 ** primal y)+    _        <**> Tower []  = auto 1     x        <**> Tower [y] = lift1 (**y) (\z -> (y *^ z <**> Tower [y-1])) x     x        <**> y         = lift2_ (**) (\z xi yi -> (yi *! z /! xi, z *! log1 xi)) x y 

@@ -0,0 +1,227 @@+{-# 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"

@@ -1,194 +0,0 @@-{-# LANGUAGE Rank2Types, TemplateHaskell, BangPatterns, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, UndecidableInstances, ScopedTypeVariables #-}--------------------------------------------------------------------------------- |--- Module      :  Numeric.AD.Mode.Chain--- Copyright   :  (c) Edward Kmett 2010--- License     :  BSD3--- Maintainer  :  ekmett@gmail.com--- Stability   :  experimental--- Portability :  GHC only------ Reverse Automatic Differentiation using Data.Reflection-----------------------------------------------------------------------------------module Numeric.AD.Mode.Chain-    (-    -- * 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.Chain-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))

@@ -33,7 +33,7 @@ import qualified Numeric.AD.Mode.Reverse as R import qualified Numeric.AD.Mode.Forward as F import qualified Numeric.AD.Mode.Tower as T-import qualified Numeric.AD.Mode.Chain as C+import qualified Numeric.AD.Mode.Wengert as W import qualified Numeric.AD as M import Data.Ix @@ -42,7 +42,7 @@ data Direction     = Forward     | Reverse-    | Chain+    | Wengert     | Tower     | Mixed     deriving (Show, Eq, Ord, Read, Bounded, Enum, Ix)@@ -50,7 +50,7 @@ diff :: Num a => Direction -> (forall s. Mode s => AD s a -> AD s a) -> a -> a diff Forward = F.diff diff Reverse = R.diff-diff Chain = C.diff+diff Wengert = W.diff diff Tower = T.diff diff Mixed = F.diff {-# INLINE diff #-}@@ -58,7 +58,7 @@ diff' :: Num a => Direction -> (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a) diff' Forward = F.diff' diff' Reverse = R.diff'-diff' Chain = C.diff'+diff' Wengert = W.diff' diff' Tower = T.diff' diff' Mixed = F.diff' {-# INLINE diff' #-}@@ -66,7 +66,7 @@ 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 Reverse = R.jacobian-jacobian Chain = C.jacobian+jacobian Wengert = W.jacobian jacobian Tower = F.jacobian -- error "jacobian Tower: unimplemented" jacobian Mixed = M.jacobian {-# INLINE jacobian #-}@@ -74,7 +74,7 @@ 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' Reverse = R.jacobian'-jacobian' Chain = C.jacobian'+jacobian' Wengert = W.jacobian' jacobian' Tower = F.jacobian' -- error "jacobian' Tower: unimplemented" jacobian' Mixed = M.jacobian' {-# INLINE jacobian' #-}@@ -82,7 +82,7 @@ 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 Reverse = R.grad-grad Chain   = C.grad+grad Wengert   = W.grad grad Tower   = F.grad -- error "grad Tower: unimplemented" grad Mixed   = M.grad {-# INLINE grad #-}@@ -90,7 +90,7 @@ 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' Reverse = R.grad'-grad' Chain   = C.grad'+grad' Wengert   = W.grad' grad' Tower   = F.grad' -- error "grad' Tower: unimplemented" grad' Mixed   = M.grad' {-# INLINE grad' #-}

@@ -159,14 +159,14 @@  -- | Compute the gradient of a function using forward mode AD. ----- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Chain.grad' for @n@ inputs, in exchange for better space utilization.+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Wengert.grad' for @n@ inputs, in exchange for better space utilization. grad :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f a grad f = bind (tangent . f) {-# INLINE grad #-}  -- | Compute the gradient and answer to a function using forward mode AD. ----- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Chain.grad'' for @n@ inputs, in exchange for better space utilization.+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Wengert.grad'' for @n@ inputs, in exchange for better space utilization. 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 b, tangent <$> bs)     where@@ -175,7 +175,7 @@  -- | Compute the gradient of a function using forward mode AD and combine the result with the input using a user-specified function. ----- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Chain.gradWith' for @n@ inputs, in exchange for better space utilization.+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Wengert.gradWith' for @n@ inputs, in exchange for better space utilization. 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 = bindWith g (tangent . f) {-# INLINE gradWith #-}@@ -183,9 +183,11 @@ -- | Compute the gradient of a function using forward mode AD and the answer, and combine the result with the input using a -- user-specified function. ----- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Chain.gradWith'' for @n@ inputs, in exchange for better space utilization.+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Wengert.gradWith'' for @n@ inputs, in exchange for better space utilization.+-- >>> gradWith' (,) sum [0..4]+-- (10,[(0,1),(1,1),(2,1),(3,1),(4,1)]) 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 = bindWith' g (tangent . f)+gradWith' g f as = (primal $ f (AD . Lift <$> as), bindWith g (tangent . f) as) {-# INLINE gradWith' #-}  -- | Compute the product of a vector with the Hessian using forward-on-forward-mode AD.

@@ -0,0 +1,194 @@+{-# 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))

@@ -64,7 +64,7 @@ -- >>> last $ take 10 $ inverse sqrt 1 (sqrt 10) -- 10.0 inverse :: (Fractional a, Eq a) => (forall s. Mode s => AD s a -> AD s a) -> a -> a -> [a]-inverse f x0 y = findZero (\x -> f x - lift y) x0+inverse f x0 y = findZero (\x -> f x - auto y) x0 {-# INLINE inverse  #-}  -- | The 'fixedPoint' function find a fixedpoint of a scalar@@ -129,7 +129,7 @@     d0 = negate <$> grad f x0     go xi ri di = xi : go xi1 ri1 di1       where-        ai  = last $ take 20 $ extremum (\a -> f $ zipWithT (\x d -> lift x + a * lift d) xi di) 0+        ai  = last $ take 20 $ extremum (\a -> f $ zipWithT (\x d -> auto x + a * auto d) xi di) 0         xi1 = zipWithT (\x d -> x + ai*d) xi di         ri1 = negate <$> grad f xi1         bi1 = max 0 $ dot ri1 (zipWithT (-) ri1 ri) / dot ri1 ri1