ad 1.4 → 1.5

raw patch · 46 files changed

+3096/−3128 lines, 46 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

• Numeric/AD.hs +0/−20
• Numeric/AD/Classes.hs +0/−16
• Numeric/AD/Halley.hs +0/−88
• Numeric/AD/Internal/Classes.hs +0/−326
• Numeric/AD/Internal/Combinators.hs +0/−28
• Numeric/AD/Internal/Composition.hs +0/−183
• Numeric/AD/Internal/Dense.hs +0/−185
• Numeric/AD/Internal/Forward.hs +0/−199
• Numeric/AD/Internal/Identity.hs +0/−139
• Numeric/AD/Internal/Reverse.hs +0/−280
• Numeric/AD/Internal/Sparse.hs +0/−255
• Numeric/AD/Internal/Tensors.hs +0/−85
• Numeric/AD/Internal/Tower.hs +0/−139
• Numeric/AD/Internal/Types.hs +0/−69
• Numeric/AD/Mode/Directed.hs +0/−93
• Numeric/AD/Mode/Forward.hs +0/−165
• Numeric/AD/Mode/Mixed.hs +0/−226
• Numeric/AD/Mode/Reverse.hs +0/−161
• Numeric/AD/Mode/Sparse.hs +0/−125
• Numeric/AD/Mode/Tower.hs +0/−128
• Numeric/AD/Newton.hs +0/−113
• Numeric/AD/Types.hs +0/−51
• ad.cabal +25/−54
• src/Numeric/AD.hs +215/−0
• src/Numeric/AD/Halley.hs +80/−0
• src/Numeric/AD/Internal/Classes.hs +328/−0
• src/Numeric/AD/Internal/Combinators.hs +28/−0
• src/Numeric/AD/Internal/Composition.hs +183/−0
• src/Numeric/AD/Internal/Dense.hs +186/−0
• src/Numeric/AD/Internal/Forward.hs +200/−0
• src/Numeric/AD/Internal/Identity.hs +139/−0
• src/Numeric/AD/Internal/Jet.hs +86/−0
• src/Numeric/AD/Internal/Reverse.hs +281/−0
• src/Numeric/AD/Internal/Sparse.hs +256/−0
• src/Numeric/AD/Internal/Tower.hs +140/−0
• src/Numeric/AD/Internal/Types.hs +59/−0
• src/Numeric/AD/Mode/Directed.hs +89/−0
• src/Numeric/AD/Mode/Forward.hs +161/−0
• src/Numeric/AD/Mode/Reverse.hs +156/−0
• src/Numeric/AD/Mode/Sparse.hs +121/−0
• src/Numeric/AD/Mode/Tower.hs +123/−0
• src/Numeric/AD/Newton.hs +108/−0
• src/Numeric/AD/Types.hs +49/−0
• src/Numeric/AD/Variadic.hs +29/−0
• src/Numeric/AD/Variadic/Reverse.hs +27/−0
• src/Numeric/AD/Variadic/Sparse.hs +27/−0

@@ -1,20 +0,0 @@--------------------------------------------------------------------------------- |--- Module      :  Numeric.AD--- Copyright   :  (c) Edward Kmett 2010--- License     :  BSD3--- Maintainer  :  ekmett@gmail.com--- Stability   :  experimental--- Portability :  GHC only------ Mixed-mode automatic differentiation combinators.-----------------------------------------------------------------------------------module Numeric.AD-    ( module Numeric.AD.Mode.Mixed-    , module Numeric.AD.Newton-    ) where--import Numeric.AD.Mode.Mixed-import Numeric.AD.Newton hiding (Mode(..), AD(..), UU, UF, FU, FF)

@@ -1,16 +0,0 @@--------------------------------------------------------------------------------- |--- Module      :  Numeric.AD.Classes--- Copyright   :  (c) Edward Kmett 2010--- License     :  BSD3--- Maintainer  :  ekmett@gmail.com--- Stability   :  experimental--- Portability :  GHC only-----------------------------------------------------------------------------------module Numeric.AD.Classes-    ( Mode(..)-    ) where--import Numeric.AD.Internal.Classes

@@ -1,88 +0,0 @@-{-# LANGUAGE Rank2Types, BangPatterns, ScopedTypeVariables #-}--------------------------------------------------------------------------------- |--- Module      :  Numeric.AD.Halley--- Copyright   :  (c) Edward Kmett 2010--- License     :  BSD3--- Maintainer  :  ekmett@gmail.com--- Stability   :  experimental--- Portability :  GHC only------ Root finding using Halley's rational method (the second in --- the class of Householder methods). Assumes the function is three --- times continuously differentiable and converges cubically when --- progress can be made.--- --------------------------------------------------------------------------------module Numeric.AD.Halley-    (-    -- * Halley's Method (Tower AD)-      findZero-    , inverse-    , fixedPoint-    , extremum-    -- * Exposed Types-    , UU, UF, FU, FF-    , AD(..)-    , Mode(..)-    ) where--import Prelude hiding (all)--- import Data.Foldable (all)--- import Data.Traversable (Traversable)-import Numeric.AD.Types-import Numeric.AD.Classes-import Numeric.AD.Mode.Tower (diffs0)-import Numeric.AD.Mode.Forward (diff) -- , diff')--- import Numeric.AD.Mode.Reverse (gradWith')-import Numeric.AD.Internal.Composition---- | The 'findZero' function finds a zero of a scalar function using--- Halley's method; its output is a stream of increasingly accurate--- results.  (Modulo the usual caveats.) ------ Examples:------  > take 10 $ findZero (\\x->x^2-4) 1  -- converge to 2.0------  > module Data.Complex---  > take 10 $ findZero ((+1).(^2)) (1 :+ 1)  -- converge to (0 :+ 1)@----findZero :: (Fractional a, Eq a) => UU a -> a -> [a]-findZero f = go-    where-        go x = x : if y == 0 then [] else go (x - 2*y*y'/(2*y'*y'-y*y''))-            where-                (y:y':y'':_) = diffs0 f x-{-# INLINE findZero #-}---- | The 'inverse' function inverts a scalar function using--- Halley's method; its output is a stream of increasingly accurate--- results.  (Modulo the usual caveats.)------ 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) => UU a -> a -> a -> [a]-inverse f x0 y = findZero (\x -> f x - lift y) x0-{-# INLINE inverse  #-}---- | The 'fixedPoint' function find a fixedpoint of a scalar--- function using Halley's method; its output is a stream of--- increasingly accurate results.  (Modulo the usual caveats.)--- --- > take 10 $ fixedPoint cos 1 -- converges to 0.7390851332151607-fixedPoint :: (Fractional a, Eq a) => UU a -> a -> [a]-fixedPoint f = findZero (\x -> f x - x)-{-# INLINE fixedPoint #-}---- | The 'extremum' function finds an extremum of a scalar--- function using Halley's method; produces a stream of increasingly--- accurate results.  (Modulo the usual caveats.)------ > take 10 $ extremum cos 1 -- convert to 0 -extremum :: (Fractional a, Eq a) => UU a -> a -> [a]-extremum f = findZero (diff (decomposeMode . f . composeMode))-{-# INLINE extremum #-}-

@@ -1,326 +0,0 @@-{-# LANGUAGE Rank2Types, TypeFamilies, FlexibleInstances, MultiParamTypeClasses, PatternGuards, CPP #-}-{-# LANGUAGE FlexibleContexts, FunctionalDependencies, UndecidableInstances, GeneralizedNewtypeDeriving, TemplateHaskell #-}--- {-# OPTIONS_HADDOCK hide #-}--------------------------------------------------------------------------------- |--- Module      :  Numeric.AD.Internal.Classes--- Copyright   :  (c) Edward Kmett 2010--- License     :  BSD3--- Maintainer  :  ekmett@gmail.com--- Stability   :  experimental--- Portability :  GHC only-----------------------------------------------------------------------------------module Numeric.AD.Internal.Classes-    (-    -- * AD modes-      Mode(..)-    , one-    -- * Automatically Deriving AD-    , Jacobian(..)-    , Primal(..)-    , deriveLifted-    , deriveNumeric-    , Lifted(..)-    , Iso(..)-    ) where--import Control.Applicative hiding ((<**>))-import Data.Char-import Language.Haskell.TH-import Numeric.AD.Internal.Combinators (on)--infixr 8 **!, <**>-infixl 7 *!, /!, ^*, *^, ^/-infixl 6 +!, -!, <+>-infix 4 ==!--class Iso a b where-    iso :: f a -> f b-    osi :: f b -> f a--instance Iso a a where-    iso = id-    osi = id--class Lifted t where-    showsPrec1          :: (Num a, Show a) => Int -> t a -> ShowS-    (==!)               :: (Num a, Eq a) => t a -> t a -> Bool-    compare1            :: (Num a, Ord a) => t a -> t a -> Ordering-    fromInteger1        :: Num a => Integer -> t a-    (+!),(-!),(*!)      :: Num a => t a -> t a -> t a-    negate1, abs1, signum1 :: Num a => t a -> t a-    (/!)                :: Fractional a => t a -> t a -> t a-    recip1              :: Fractional a => t a -> t a-    fromRational1       :: Fractional a => Rational -> t a-    toRational1         :: Real a => t a -> Rational -- unsafe-    pi1                 :: Floating a => t a-    exp1, log1, sqrt1   :: Floating a => t a -> t a-    (**!), logBase1     :: Floating a => t a -> t a -> t a-    sin1, cos1, tan1, asin1, acos1, atan1 :: Floating a => t a -> t a-    sinh1, cosh1, tanh1, asinh1, acosh1, atanh1 :: Floating a => t a -> t a-    properFraction1 :: (RealFrac a, Integral b) => t a -> (b, t a)-    truncate1, round1, ceiling1, floor1 :: (RealFrac a, Integral b) => t a -> b-    floatRadix1     :: RealFloat a => t a -> Integer-    floatDigits1    :: RealFloat a => t a -> Int-    floatRange1     :: RealFloat a => t a -> (Int, Int)-    decodeFloat1    :: RealFloat a => t a -> (Integer, Int)-    encodeFloat1    :: RealFloat a => Integer -> Int -> t a-    exponent1       :: RealFloat a => t a -> Int-    significand1    :: RealFloat a => t a -> t a-    scaleFloat1     :: RealFloat a => Int -> t a -> t a-    isNaN1, isInfinite1, isDenormalized1, isNegativeZero1, isIEEE1 :: RealFloat a => t a -> Bool-    atan21          :: RealFloat a => t a -> t a -> t a-    succ1, pred1    :: (Num a, Enum a) => t a -> t a-    toEnum1         :: (Num a, Enum a) => Int -> t a-    fromEnum1       :: (Num a, Enum a) => t a -> Int-    enumFrom1       :: (Num a, Enum a) => t a -> [t a]-    enumFromThen1   :: (Num a, Enum a) => t a -> t a -> [t a]-    enumFromTo1     :: (Num a, Enum a) => t a -> t a -> [t a]-    enumFromThenTo1 :: (Num a, Enum a) => t a -> t a -> t a -> [t a]-    minBound1       :: (Num a, Bounded a) => t a-    maxBound1       :: (Num a, Bounded a) => t a--class Lifted t => Mode t where-    -- | allowed to return False for items with a zero derivative, but we'll give more NaNs than strictly necessary-    isKnownConstant :: t a -> Bool-    isKnownConstant _ = False--    -- | allowed to return False for zero, but we give more NaN's than strictly necessary then-    isKnownZero :: Num a => t a -> Bool-    isKnownZero _ = False--    -- | Embed a constant-    lift  :: Num a => a -> t a--    -- | Vector sum-    (<+>) :: Num a => t a -> t a -> t a--    -- | Scalar-vector multiplication-    (*^) :: Num a => a -> t a -> t a--    -- | Vector-scalar multiplication-    (^*) :: Num a => t a -> a -> t a--    -- | Scalar division-    (^/) :: Fractional a => t a -> a -> t a--    -- | Exponentiation, this should be overloaded if you can figure out anything about what is constant!-    (<**>) :: Floating a => t a -> t a -> t a---  x <**> y = lift2_ (**) (\z xi yi -> (yi *! z /! xi, z *! log1 xi)) x y--    -- | > 'zero' = 'lift' 0-    zero :: Num a => t a--    a *^ b = lift a *! b-    a ^* b = a *! lift b--    a ^/ b = a ^* recip b--    zero = lift 0--one :: (Mode t, Num a) => t a-one = lift 1-{-# INLINE one #-}--negOne :: (Mode t, Num a) => t a-negOne = lift (-1)-{-# INLINE negOne #-}---- | 'Primal' is used by 'deriveMode' but is not exposed--- via the 'Mode' class to prevent its abuse by end users--- 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--- information. The end user is protected from accidentally using this function--- by the universal quantification on the various combinators we expose.--class Primal t where-    primal :: Num a => t a -> a---- | 'Jacobian' is used by 'deriveMode' but is not exposed--- via 'Mode' to prevent its abuse by end users--- via the 'AD' data type.-class (Mode t, Mode (D t)) => Jacobian t where-    type D t :: * -> *--    unary  :: Num a => (a -> a) -> D t a -> t a -> t a-    lift1  :: Num a => (a -> a) -> (D t a -> D t a) -> t a -> t a-    lift1_ :: Num a => (a -> a) -> (D t a -> D t a -> D t a) -> t a -> t a--    binary :: Num a => (a -> a -> a) -> D t a -> D t a -> t a -> t a -> t a-    lift2  :: Num a => (a -> a -> a) -> (D t a -> D t a -> (D t a, D t a)) -> t a -> t a -> t a-    lift2_ :: Num a => (a -> a -> a) -> (D t a -> D t a -> D t a -> (D t a, D t a)) -> t a -> t a -> t a--withPrimal :: (Jacobian t, Num a) => t a -> a -> t a-withPrimal t a = unary (const a) one t-{-# INLINE withPrimal #-}--fromBy :: (Jacobian t, Num a) => t a -> t a -> Int -> a -> t a-fromBy a delta n x = binary (\_ _ -> x) one (fromIntegral1 n) a delta--fromIntegral1 :: (Integral n, Lifted t, Num a) => n -> t a-fromIntegral1 = fromInteger1 . fromIntegral-{-# INLINE fromIntegral1 #-}--square1 :: (Lifted t, Num a) => t a -> t a-square1 x = x *! x-{-# INLINE square1 #-}--discrete1 :: (Primal t, Num a) => (a -> c) -> t a -> c-discrete1 f x = f (primal x)-{-# INLINE discrete1 #-}--discrete2 :: (Primal t, Num a) => (a -> a -> c) -> t a -> t a -> c-discrete2 f x y = f (primal x) (primal y)-{-# INLINE discrete2 #-}--discrete3 :: (Primal t, Num a) => (a -> a -> a -> d) -> t a -> t a -> t a -> d-discrete3 f x y z = f (primal x) (primal y) (primal z)-{-# INLINE discrete3 #-}---- | @'deriveLifted' t@ provides------ > instance Lifted $t------ given supplied instances for------ > instance Lifted $t => Primal $t where ...--- > instance Lifted $t => Jacobian $t where ...------ The seemingly redundant @'Lifted' $t@ constraints are caused by Template Haskell staging restrictions.-deriveLifted :: ([Q Pred] -> [Q Pred]) -> Q Type -> Q [Dec]-deriveLifted f _t = do-        [InstanceD cxt0 type0 dec0] <- lifted-        return <$> instanceD (cxt (f (return <$> cxt0))) (return type0) (return <$> dec0)-    where-      lifted = [d|-       instance Lifted $_t where-        (==!)         = (==) `on` primal-        compare1      = compare `on` primal-        maxBound1     = lift maxBound-        minBound1     = lift minBound-        showsPrec1 d  = showsPrec d . primal-        fromInteger1  = lift . fromInteger-        (+!)          = (<+>) -- binary (+) one one-        (-!)          = binary (-) one negOne -- TODO: <-> ? as it is, this might be pretty bad for Tower-        (*!)          = lift2 (*) (\x y -> (y, x))-        negate1       = lift1 negate (const negOne)-        abs1          = lift1 abs signum1-        signum1       = lift1 signum (const zero)-        fromRational1 = lift . fromRational-        x /! y        = x *! recip1 y-        recip1        = lift1_ recip (const . negate1 . square1)-        pi1       = lift pi-        exp1      = lift1_ exp const-        log1      = lift1 log recip1-        logBase1 x y = log1 y /! log1 x-        sqrt1     = lift1_ sqrt (\z _ -> recip1 (lift 2 *! z))-        (**!)     = (<**>)-        --x **! y-        --   | isKnownZero y     = 1-        --   | isKnownConstant y, y' <- primal y = lift1 (** y') ((y'*) . (**(y'-1))) x-        --   | otherwise         = lift2_ (**) (\z xi yi -> (yi *! z /! xi, z *! log1 xi)) x y-        sin1      = lift1 sin cos1-        cos1      = lift1 cos $ negate1 . sin1-        tan1 x    = sin1 x /! cos1 x-        asin1     = lift1 asin $ \x -> recip1 (sqrt1 (one -! square1 x))-        acos1     = lift1 acos $ \x -> negate1 (recip1 (sqrt1 (one -! square1 x)))-        atan1     = lift1 atan $ \x -> recip1 (one +! square1 x)-        sinh1     = lift1 sinh cosh1-        cosh1     = lift1 cosh sinh1-        tanh1 x   = sinh1 x /! cosh1 x-        asinh1    = lift1 asinh $ \x -> recip1 (sqrt1 (one +! square1 x))-        acosh1    = lift1 acosh $ \x -> recip1 (sqrt1 (square1 x -! one))-        atanh1    = lift1 atanh $ \x -> recip1 (one -! square1 x)--        succ1                 = lift1 succ (const one)-        pred1                 = lift1 pred (const one)-        toEnum1               = lift . toEnum-        fromEnum1             = discrete1 fromEnum-        enumFrom1 a           = withPrimal a <$> discrete1 enumFrom a-        enumFromTo1 a b       = withPrimal a <$> discrete2 enumFromTo a b-        enumFromThen1 a b     = zipWith (fromBy a delta) [0..] $ discrete2 enumFromThen a b where delta = b -! a-        enumFromThenTo1 a b c = zipWith (fromBy a delta) [0..] $ discrete3 enumFromThenTo a b c where delta = b -! a--        toRational1      = discrete1 toRational-        floatRadix1      = discrete1 floatRadix-        floatDigits1     = discrete1 floatDigits-        floatRange1      = discrete1 floatRange-        decodeFloat1     = discrete1 decodeFloat-        encodeFloat1 m e = lift (encodeFloat m e)-        isNaN1           = discrete1 isNaN-        isInfinite1      = discrete1 isInfinite-        isDenormalized1  = discrete1 isDenormalized-        isNegativeZero1  = discrete1 isNegativeZero-        isIEEE1          = discrete1 isIEEE-        exponent1 = exponent . primal-        scaleFloat1 n = unary (scaleFloat n) (scaleFloat1 n one)-        significand1 x =  unary significand (scaleFloat1 (- floatDigits1 x) one) x-        atan21 = lift2 atan2 $ \vx vy -> let r = recip1 (square1 vx +! square1 vy) in (vy *! r, negate1 vx *! r)-        properFraction1 a = (w, a `withPrimal` pb) where-             pa = primal a-             (w, pb) = properFraction pa-        truncate1 = discrete1 truncate-        round1    = discrete1 round-        ceiling1  = discrete1 ceiling-        floor1    = discrete1 floor |]--varA :: Q Type-varA = varT (mkName "a")---- | Find all the members defined in the 'Lifted' data type-liftedMembers :: Q [String]-liftedMembers = do-#ifdef OldClassI-    ClassI (ClassD _ _ _ _ ds) <- reify ''Lifted-#else-    ClassI (ClassD _ _ _ _ ds) _ <- reify ''Lifted-#endif-    return [ nameBase n | SigD n _ <- ds]---- | @'deriveNumeric' f g@ provides the following instances:------ > instance ('Lifted' $f, 'Num' a, 'Enum' a) => 'Enum' ($g a)--- > instance ('Lifted' $f, 'Num' a, 'Eq' a) => 'Eq' ($g a)--- > instance ('Lifted' $f, 'Num' a, 'Ord' a) => 'Ord' ($g a)--- > instance ('Lifted' $f, 'Num' a, 'Bounded' a) => 'Bounded' ($g a)------ > instance ('Lifted' $f, 'Show' a) => 'Show' ($g a)--- > instance ('Lifted' $f, 'Num' a) => 'Num' ($g a)--- > instance ('Lifted' $f, 'Fractional' a) => 'Fractional' ($g a)--- > instance ('Lifted' $f, 'Floating' a) => 'Floating' ($g a)--- > instance ('Lifted' $f, 'RealFloat' a) => 'RealFloat' ($g a)--- > instance ('Lifted' $f, 'RealFrac' a) => 'RealFrac' ($g a)--- > instance ('Lifted' $f, 'Real' a) => 'Real' ($g a)-deriveNumeric :: ([Q Pred] -> [Q Pred]) -> Q Type -> Q [Dec]-deriveNumeric f t = do-    members <- liftedMembers-    let keep n = nameBase n `elem` members-    xs <- lowerInstance keep ((classP ''Num [varA]:) . f) t `mapM` [''Enum, ''Eq, ''Ord, ''Bounded, ''Show]-    ys <- lowerInstance keep f                            t `mapM` [''Num, ''Fractional, ''Floating, ''RealFloat,''RealFrac, ''Real]-    return (xs ++ ys)--lowerInstance :: (Name -> Bool) -> ([Q Pred] -> [Q Pred]) -> Q Type -> Name -> Q Dec-lowerInstance p f t n = do-#ifdef OldClassI-    ClassI (ClassD _ _ _ _ ds) <- reify n-#else-    ClassI (ClassD _ _ _ _ ds) _ <- reify n-#endif-    instanceD (cxt (f [classP n [varA]]))-              (conT n `appT` (t `appT` varA))-              (concatMap lower1 ds)-    where-        lower1 :: Dec -> [Q Dec]-        lower1 (SigD n' _) | p n'' = [valD (varP n') (normalB (varE n'')) []] where n'' = primed n'-        lower1 _          = []--        primed n' = mkName $ base ++ [prime]-            where-                base = nameBase n'-                h = head base-                prime | isSymbol h || h `elem` "/*-<>" = '!'-                      | otherwise = '1'

@@ -1,28 +0,0 @@-{-# OPTIONS_GHC -fno-warn-incomplete-patterns #-}--------------------------------------------------------------------------------- |--- Module      :  Numeric.AD.Internal.Combinators--- Copyright   :  (c) Edward Kmett 2010--- License     :  BSD3--- Maintainer  :  ekmett@gmail.com--- Stability   :  experimental--- Portability :  GHC only----------------------------------------------------------------------------------module Numeric.AD.Internal.Combinators-    ( zipWithT-    , zipWithDefaultT-    , on-    ) where--import Data.Traversable (Traversable, mapAccumL)-import Data.Foldable (Foldable, toList)--on :: (a -> a -> b) -> (c -> a) -> c -> c -> b-on f g a b = f (g a) (g b)--zipWithT :: (Foldable f, Traversable g) => (a -> b -> c) -> f a -> g b -> g c-zipWithT f as = snd . mapAccumL (\(a:as') b -> (as', f a b)) (toList as)--zipWithDefaultT :: (Foldable f, Traversable g) => a -> (a -> b -> c) -> f a -> g b -> g c-zipWithDefaultT z f as = zipWithT f (toList as ++ repeat z)

@@ -1,183 +0,0 @@-{-# LANGUAGE Rank2Types, TypeFamilies, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, FlexibleContexts, TemplateHaskell, UndecidableInstances, TypeOperators #-}--- {-# OPTIONS_HADDOCK hide, prune #-}--------------------------------------------------------------------------------- |--- Module      :  Numeric.AD.Internal.Composition--- Copyright   :  (c) Edward Kmett 2010--- License     :  BSD3--- Maintainer  :  ekmett@gmail.com--- Stability   :  experimental--- Portability :  GHC only-----------------------------------------------------------------------------------module Numeric.AD.Internal.Composition-    ( ComposeFunctor(..)-    , ComposeMode(..)-    , composeMode-    , decomposeMode-    ) where--import Control.Applicative hiding ((<**>))-import Data.Data (Data(..), mkDataType, DataType, mkConstr, Constr, constrIndex, Fixity(..))-import Data.Typeable (Typeable1(..), Typeable(..), TyCon, mkTyCon3, mkTyConApp, typeOfDefault, gcast1)-import Data.Foldable (Foldable(foldMap))-import Data.Traversable (Traversable(traverse))-import Numeric.AD.Internal.Classes-import Numeric.AD.Internal.Types---- | Functor composition, used to nest the use of jacobian and grad-newtype ComposeFunctor f g a = ComposeFunctor { decomposeFunctor :: f (g a) }--instance (Functor f, Functor g) => Functor (ComposeFunctor f g) where-    fmap f (ComposeFunctor a) = ComposeFunctor (fmap (fmap f) a)--instance (Foldable f, Foldable g) => Foldable (ComposeFunctor f g) where-    foldMap f (ComposeFunctor a) = foldMap (foldMap f) a--instance (Traversable f, Traversable g) => Traversable (ComposeFunctor f g) where-    traverse f (ComposeFunctor a) = ComposeFunctor <$> traverse (traverse f) a--instance (Typeable1 f, Typeable1 g) => Typeable1 (ComposeFunctor f g) where-    typeOf1 tfga = mkTyConApp composeFunctorTyCon [typeOf1 (fa tfga), typeOf1 (ga tfga)]-        where fa :: t f (g :: * -> *) a -> f a-              fa = undefined-              ga :: t (f :: * -> *) g a -> g a-              ga = undefined--composeFunctorTyCon :: TyCon-composeFunctorTyCon = mkTyCon3 "ad" "Numeric.AD.Internal.Composition" "ComposeFunctor"-{-# NOINLINE composeFunctorTyCon #-}--composeFunctorConstr :: Constr-composeFunctorConstr = mkConstr composeFunctorDataType "ComposeFunctor" [] Prefix-{-# NOINLINE composeFunctorConstr #-}--composeFunctorDataType :: DataType-composeFunctorDataType = mkDataType "Numeric.AD.Internal.Composition.ComposeFunctor" [composeFunctorConstr]-{-# NOINLINE composeFunctorDataType #-}--instance (Typeable1 f, Typeable1 g, Data (f (g a)), Data a) => Data (ComposeFunctor f g a) where-    gfoldl f z (ComposeFunctor a) = z ComposeFunctor `f` a-    toConstr _ = composeFunctorConstr-    gunfold k z c = case constrIndex c of-        1 -> k (z ComposeFunctor)-        _ -> error "gunfold"-    dataTypeOf _ = composeFunctorDataType-    dataCast1 f = gcast1 f---- | The composition of two AD modes is an AD mode in its own right-newtype ComposeMode f g a = ComposeMode { runComposeMode :: f (AD g a) }--composeMode :: AD f (AD g a) -> AD (ComposeMode f g) a-composeMode (AD a) = AD (ComposeMode a)--decomposeMode :: AD (ComposeMode f g) a -> AD f (AD g a)-decomposeMode (AD (ComposeMode a)) = AD a--instance (Primal f, Mode g, Primal g) => Primal (ComposeMode f g) where-    primal = primal . primal . runComposeMode--instance (Mode f, Mode g) => Mode (ComposeMode f g) where-    lift = ComposeMode . lift . lift-    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)-    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-    ComposeMode a +! ComposeMode b = ComposeMode (a +! b)-    ComposeMode a -! ComposeMode b = ComposeMode (a -! b)-    ComposeMode a *! ComposeMode b = ComposeMode (a *! b)-    negate1 (ComposeMode a) = ComposeMode (negate1 a)-    abs1 (ComposeMode a) = ComposeMode (abs1 a)-    signum1 (ComposeMode a) = ComposeMode (signum1 a)-    ComposeMode a /! ComposeMode b = ComposeMode (a /! b)-    recip1 (ComposeMode a) = ComposeMode (recip1 a)-    fromRational1 = ComposeMode . lift . fromRational1-    toRational1 (ComposeMode a) = toRational1 a-    pi1 = ComposeMode pi1-    exp1 (ComposeMode a) = ComposeMode (exp1 a)-    log1 (ComposeMode a) = ComposeMode (log1 a)-    sqrt1 (ComposeMode a) = ComposeMode (sqrt1 a)-    ComposeMode a **! ComposeMode b = ComposeMode (a **! b)-    logBase1 (ComposeMode a) (ComposeMode b) = ComposeMode (logBase1 a b)-    sin1 (ComposeMode a) = ComposeMode (sin1 a)-    cos1 (ComposeMode a) = ComposeMode (cos1 a)-    tan1 (ComposeMode a) = ComposeMode (tan1 a)-    asin1 (ComposeMode a) = ComposeMode (asin1 a)-    acos1 (ComposeMode a) = ComposeMode (acos1 a)-    atan1 (ComposeMode a) = ComposeMode (atan1 a)-    sinh1 (ComposeMode a) = ComposeMode (sinh1 a)-    cosh1 (ComposeMode a) = ComposeMode (cosh1 a)-    tanh1 (ComposeMode a) = ComposeMode (tanh1 a)-    asinh1 (ComposeMode a) = ComposeMode (asinh1 a)-    acosh1 (ComposeMode a) = ComposeMode (acosh1 a)-    atanh1 (ComposeMode a) = ComposeMode (atanh1 a)-    properFraction1 (ComposeMode a) = (b, ComposeMode c) where-        (b, c) = properFraction1 a-    truncate1 (ComposeMode a) = truncate1 a-    round1 (ComposeMode a) = round1 a-    ceiling1 (ComposeMode a) = ceiling1 a-    floor1 (ComposeMode a) = floor1 a-    floatRadix1 (ComposeMode a) = floatRadix1 a-    floatDigits1 (ComposeMode a) = floatDigits1 a-    floatRange1 (ComposeMode a) = floatRange1 a-    decodeFloat1 (ComposeMode a) = decodeFloat1 a-    encodeFloat1 m e = ComposeMode (encodeFloat1 m e)-    exponent1 (ComposeMode a) = exponent1 a-    significand1 (ComposeMode a) = ComposeMode (significand1 a)-    scaleFloat1 n (ComposeMode a) = ComposeMode (scaleFloat1 n a)-    isNaN1 (ComposeMode a) = isNaN1 a-    isInfinite1 (ComposeMode a) = isInfinite1 a-    isDenormalized1 (ComposeMode a) = isDenormalized1 a-    isNegativeZero1 (ComposeMode a) = isNegativeZero1 a-    isIEEE1 (ComposeMode a) = isIEEE1 a-    atan21 (ComposeMode a) (ComposeMode b) = ComposeMode (atan21 a b)-    succ1 (ComposeMode a) = ComposeMode (succ1 a)-    pred1 (ComposeMode a) = ComposeMode (pred1 a)-    toEnum1 n = ComposeMode (toEnum1 n)-    fromEnum1 (ComposeMode a) = fromEnum1 a-    enumFrom1 (ComposeMode a) = map ComposeMode $ enumFrom1 a-    enumFromThen1 (ComposeMode a) (ComposeMode b) = map ComposeMode $ enumFromThen1 a b-    enumFromTo1 (ComposeMode a) (ComposeMode b) = map ComposeMode $ enumFromTo1 a b-    enumFromThenTo1 (ComposeMode a) (ComposeMode b) (ComposeMode c) = map ComposeMode $ enumFromThenTo1 a b c-    minBound1 = ComposeMode minBound1-    maxBound1 = ComposeMode maxBound1--instance (Typeable1 f, Typeable1 g) => Typeable1 (ComposeMode f g) where-    typeOf1 tfga = mkTyConApp composeModeTyCon [typeOf1 (fa tfga), typeOf1 (ga tfga)]-        where fa :: t f (g :: * -> *) a -> f a-              fa = undefined-              ga :: t (f :: * -> *) g a -> g a-              ga = undefined--instance (Typeable1 f, Typeable1 g, Typeable a) => Typeable (ComposeMode f g a) where-    typeOf = typeOfDefault-    -composeModeTyCon :: TyCon-composeModeTyCon = mkTyCon3 "ad" "Numeric.AD.Internal.Composition" "ComposeMode"-{-# NOINLINE composeModeTyCon #-}--composeModeConstr :: Constr-composeModeConstr = mkConstr composeModeDataType "ComposeMode" [] Prefix-{-# NOINLINE composeModeConstr #-}--composeModeDataType :: DataType-composeModeDataType = mkDataType "Numeric.AD.Internal.Composition.ComposeMode" [composeModeConstr]-{-# NOINLINE composeModeDataType #-}--instance (Typeable1 f, Typeable1 g, Data (f (AD g a)), Data a) => Data (ComposeMode f g a) where-    gfoldl f z (ComposeMode a) = z ComposeMode `f` a-    toConstr _ = composeModeConstr-    gunfold k z c = case constrIndex c of-        1 -> k (z ComposeMode)-        _ -> error "gunfold"-    dataTypeOf _ = composeModeDataType-    dataCast1 f = gcast1 f-

@@ -1,185 +0,0 @@-{-# LANGUAGE Rank2Types, TypeFamilies, FlexibleContexts, UndecidableInstances, TemplateHaskell, DeriveDataTypeable, BangPatterns #-}--- {-# OPTIONS_HADDOCK hide, prune #-}--------------------------------------------------------------------------------- |--- Module      : Numeric.AD.Internal.Dense--- Copyright   : (c) Edward Kmett 2010--- License     : BSD3--- Maintainer  : ekmett@gmail.com--- Stability   : experimental--- Portability : GHC only------ Dense Forward AD. Useful when the result involves the majority of the input--- elements. Do not use for 'Numeric.AD.Mode.Mixed.hessian' and beyond, since--- they only contain a small number of unique @n@th derivatives ----- @(n + k - 1) `choose` k@ for functions of @k@ inputs rather than the--- @k^n@ that would be generated by using 'Dense', not to mention the redundant--- intermediate derivatives that would be--- calculated over and over during that process!------ Assumes all instances of 'f' have the same number of elements.------ NB: We don't need the full power of 'Traversable' here, we could get--- by with a notion of zippable that can plug in 0's for the missing--- entries. This might allow for gradients where @f@ has exponentials like @((->) a)@--------------------------------------------------------------------------------module Numeric.AD.Internal.Dense-    ( Dense(..)-    , ds-    , ds'-    , vars-    , apply-    ) where--import Language.Haskell.TH-import Data.Typeable ()-import Data.Traversable (Traversable, mapAccumL)-import Data.Data ()-import Numeric.AD.Internal.Types-import Numeric.AD.Internal.Combinators-import Numeric.AD.Internal.Classes-import Numeric.AD.Internal.Identity--data Dense f a-    = Lift !a-    | Dense !a (f a)-    | Zero--instance Show a => Show (Dense f a) where-    showsPrec d (Lift a)    = showsPrec d a-    showsPrec d (Dense a _) = showsPrec d a-    showsPrec _ Zero        = showString "0"--ds :: f a -> AD (Dense f) a -> f a-ds _ (AD (Dense _ da)) = da-ds z _ = z-{-# INLINE ds #-}--ds' :: Num a => f a -> AD (Dense f) a -> (a, f a)-ds' _ (AD (Dense a da)) = (a, da)-ds' z (AD (Lift a)) = (a, z)-ds' z (AD Zero) = (0, z)-{-# INLINE ds' #-}---- Bind variables and count inputs-vars :: (Traversable f, Num a) => f a -> f (AD (Dense f) a)-vars as = snd $ mapAccumL outer (0 :: Int) as-    where-        outer !i a = (i + 1, AD $ Dense a $ snd $ mapAccumL (inner i) 0 as)-        inner !i !j _ = (j + 1, if i == j then 1 else 0)-{-# INLINE vars #-}--apply :: (Traversable f, Num a) => (f (AD (Dense f) a) -> b) -> f a -> b-apply f as = f (vars as)-{-# INLINE apply #-}--instance Primal (Dense f) where-    primal Zero = 0-    primal (Lift a) = a-    primal (Dense a _) = a--instance (Traversable f, Lifted (Dense f)) => Mode (Dense f) where-    lift = Lift-    zero = Zero--    Zero <+> a = a-    a <+> Zero = a-    Lift a     <+> Lift b     = Lift (a + b)-    Lift a     <+> Dense b db = Dense (a + b) db-    Dense a da <+> Lift b     = Dense (a + b) da-    Dense a da <+> Dense b db = Dense (a + b) $ zipWithT (+) da db--    _ <**> 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--    _ *^ Zero       = Zero-    a *^ Lift b     = Lift (a * b)-    a *^ Dense b db = Dense (a * b) $ fmap (a*) db-    Zero       ^* _ = Zero-    Lift a     ^* b = Lift (a * b)-    Dense a da ^* b = Dense (a * b) $ fmap (*b) da-    Zero       ^/ _ = Zero-    Lift a     ^/ b = Lift (a / b)-    Dense a da ^/ b = Dense (a / b) $ fmap (/b) da--instance (Traversable f, Lifted (Dense f)) => Jacobian (Dense f) where-    type D (Dense f) = Id-    unary f _         Zero        = Lift (f 0)-    unary f _         (Lift b)    = Lift (f b)-    unary f (Id dadb) (Dense b db) = Dense (f b) (fmap (dadb *) db)--    lift1 f _  Zero        = Lift (f 0)-    lift1 f _  (Lift b)    = Lift (f b)-    lift1 f df (Dense b db) = Dense (f b) (fmap (dadb *) db)-        where-            Id dadb = df (Id b)--    lift1_ f _  Zero         = Lift (f 0)-    lift1_ f _  (Lift b)     = Lift (f b)-    lift1_ f df (Dense b db) = Dense a (fmap (dadb *) db)-        where-            a = f b-            Id dadb = df (Id a) (Id b)--    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         (Dense c dc) = Dense (f 0 c) $ fmap (* dadc) dc-    binary f _         (Id dadc) (Lift b)     (Dense c dc) = Dense (f b c) $ fmap (* dadc) dc-    binary f (Id dadb) _         (Dense b db) Zero         = Dense (f b 0) $ fmap (dadb *) db-    binary f (Id dadb) _         (Dense b db) (Lift c)     = Dense (f b c) $ fmap (dadb *) db-    binary f (Id dadb) (Id dadc) (Dense b db) (Dense c dc) = Dense (f b c) $ zipWithT productRule db dc-        where productRule dbi dci = dadb * dbi + dci * dadc--    lift2 f _  Zero         Zero         = Lift (f 0 0)-    lift2 f _  Zero         (Lift c)     = Lift (f 0 c)-    lift2 f _  (Lift b)     Zero         = Lift (f b 0)-    lift2 f _  (Lift b)     (Lift c)     = Lift (f b c)-    lift2 f df Zero         (Dense c dc) = Dense (f 0 c) $ fmap (*dadc) dc where dadc = runId (snd (df (Id 0) (Id c)))-    lift2 f df (Lift b)     (Dense c dc) = Dense (f b c) $ fmap (*dadc) dc where dadc = runId (snd (df (Id b) (Id c)))-    lift2 f df (Dense b db) Zero         = Dense (f b 0) $ fmap (dadb*) db where dadb = runId (fst (df (Id b) (Id 0)))-    lift2 f df (Dense b db) (Lift c)     = Dense (f b c) $ fmap (dadb*) db where dadb = runId (fst (df (Id b) (Id c)))-    lift2 f df (Dense b db) (Dense c dc) = Dense (f b c) da-        where-            (Id dadb, Id dadc) = df (Id b) (Id c)-            da = zipWithT productRule db dc-            productRule dbi dci = dadb * dbi + dci * dadc--    lift2_ f _  Zero     Zero     = Lift (f 0 0)-    lift2_ f _  Zero     (Lift c) = Lift (f 0 c)-    lift2_ f _  (Lift b) Zero     = Lift (f b 0)-    lift2_ f _  (Lift b) (Lift c) = Lift (f b c)-    lift2_ f df Zero     (Dense c dc)-        = Dense a $ fmap (*dadc) dc-        where-            a = f 0 c-            (_, Id dadc) = df (Id a) (Id 0) (Id c)-    lift2_ f df (Lift b) (Dense c dc)-        = Dense a $ fmap (*dadc) dc-        where-            a = f b c-            (_, Id dadc) = df (Id a) (Id b) (Id c)-    lift2_ f df (Dense b db) Zero-        = Dense a $ fmap (dadb*) db-        where-            a = f b 0-            (Id dadb, _) = df (Id a) (Id b) (Id 0)-    lift2_ f df (Dense b db) (Lift c)-        = Dense a $ fmap (dadb*) db-        where-            a = f b c-            (Id dadb, _) = df (Id a) (Id b) (Id c)-    lift2_ f df (Dense b db) (Dense c dc)-        = Dense a $ zipWithT productRule db dc-        where-            a = f b c-            (Id dadb, Id dadc) = df (Id a) (Id b) (Id c)-            productRule dbi dci = dadb * dbi + dci * dadc--let f = varT (mkName "f") in-    deriveLifted-        (classP ''Traversable [f]:)-        (conT ''Dense `appT` f)

@@ -1,199 +0,0 @@-{-# LANGUAGE Rank2Types, TypeFamilies, DeriveDataTypeable, TemplateHaskell, UndecidableInstances, BangPatterns #-}--- {-# OPTIONS_HADDOCK hide, prune #-}--------------------------------------------------------------------------------- |--- Module      :  Numeric.AD.Internal.Forward--- Copyright   :  (c) Edward Kmett 2010--- License     :  BSD3--- Maintainer  :  ekmett@gmail.com--- Stability   :  experimental--- Portability :  GHC only------ Unsafe and often partial combinators intended for internal usage.------ Handle with care.--------------------------------------------------------------------------------module Numeric.AD.Internal.Forward-    ( Forward(..)-    , tangent-    , bundle-    , unbundle-    , apply-    , bind-    , bind'-    , bindWith-    , bindWith'-    , transposeWith-    ) where--import Language.Haskell.TH-import Data.Typeable-import Data.Traversable (Traversable, mapAccumL)-import Data.Foldable (Foldable, toList)-import Data.Data-import Control.Applicative-import Numeric.AD.Internal.Types-import Numeric.AD.Internal.Classes-import Numeric.AD.Internal.Identity--data Forward a-  = Forward !a a-  | Lift !a-  | Zero-  deriving (Show, Data, Typeable)--tangent :: Num a => AD Forward a -> a-tangent (AD (Forward _ da)) = da-tangent _ = 0-{-# INLINE tangent #-}--unbundle :: Num a => AD Forward a -> (a, a)-unbundle (AD (Forward a da)) = (a, da)-unbundle (AD Zero) = (0,0)-unbundle (AD (Lift a)) = (a, 0)-{-# INLINE unbundle #-}--bundle :: a -> a -> AD Forward a-bundle a da = AD (Forward a da)-{-# INLINE bundle #-}--apply :: Num a => (AD Forward a -> b) -> a -> b-apply f a = f (bundle a 1)-{-# INLINE apply #-}--instance Primal Forward where-    primal (Forward a _) = a-    primal (Lift a) = a-    primal Zero = 0--instance Lifted Forward => Mode Forward where-    lift = Lift-    zero = Zero--    isKnownZero Zero = True-    isKnownZero _    = False--    isKnownConstant Forward{} = False-    isKnownConstant _ = True--    Zero <+> a = a-    a <+> Zero = a-    Forward a da <+> Forward b db = Forward (a + b) (da + db)-    Forward a da <+> Lift b = Forward (a + b) da-    Lift a <+> Forward b db = Forward (a + b) db-    Lift a <+> Lift b = Lift (a + b)--    _ <**> 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--    a *^ Forward b db = Forward (a * b) (a * db)-    a *^ Lift b = Lift (a * b)-    _ *^ Zero = Zero--    Forward a da ^* b = Forward (a * b) (da * b)-    Lift a ^* b = Lift (a * b)-    Zero ^* _ = Zero--    Forward a da ^/ b = Forward (a / b) (da / b)-    Lift a ^/ b = Lift (a / b)-    Zero ^/ _ = Zero--instance Lifted Forward => Jacobian Forward where-    type D Forward = Id---    unary f (Id dadb) (Forward b db) = Forward (f b) (dadb * db)-    unary f _         (Lift b)       = Lift (f b)-    unary f _         Zero           = Lift (f 0)--    lift1 f _ Zero            = Lift (f 0)-    lift1 f _  (Lift b)       = Lift (f b)-    lift1 f df (Forward b db) = Forward (f b) (dadb * db)-        where-            Id dadb = df (Id b)--    lift1_ f _  Zero           = Lift (f 0)-    lift1_ f _  (Lift b)       = Lift (f b)-    lift1_ f df (Forward b db) = Forward a da-        where-            a = f b-            Id da = df (Id a) (Id b) ^* db--    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           (Forward c dc) = Forward (f 0 c) $ dc * dadc-    binary f _         (Id dadc) (Lift b)       (Forward c dc) = Forward (f b c) $ dc * dadc-    binary f (Id dadb) _         (Forward b db) Zero           = Forward (f b 0) $ dadb * db-    binary f (Id dadb) _         (Forward b db) (Lift c)       = Forward (f b c) $ dadb * db-    binary f (Id dadb) (Id dadc) (Forward b db) (Forward c dc) = Forward (f b c) $ dadb * db + dc * dadc--    lift2 f _  Zero           Zero           = Lift (f 0 0)-    lift2 f _  Zero           (Lift c)       = Lift (f 0 c)-    lift2 f _  (Lift b)       Zero           = Lift (f b 0)-    lift2 f _  (Lift b)       (Lift c)       = Lift (f b c)-    lift2 f df Zero           (Forward c dc) = Forward (f 0 c) $ dc * runId (snd (df (Id 0) (Id c)))-    lift2 f df (Lift b)       (Forward c dc) = Forward (f b c) $ dc * runId (snd (df (Id b) (Id c)))-    lift2 f df (Forward b db) Zero           = Forward (f b 0) $ runId (fst (df (Id b) (Id 0))) * db-    lift2 f df (Forward b db) (Lift c)       = Forward (f b c) $ runId (fst (df (Id b) (Id c))) * db-    lift2 f df (Forward b db) (Forward c dc) = Forward a da-        where-            a = f b c-            (Id dadb, Id dadc) = df (Id b) (Id c)-            da = dadb * db + dc * dadc--    lift2_ f _  Zero           Zero           = Lift (f 0 0)-    lift2_ f _  Zero           (Lift c)       = Lift (f 0 c)-    lift2_ f _  (Lift b)       Zero           = Lift (f b 0)-    lift2_ f _  (Lift b)       (Lift c)       = Lift (f b c)-    lift2_ f df Zero           (Forward c dc) = Forward a $ dc * runId (snd (df (Id a) (Id 0) (Id c))) where a = f 0 c-    lift2_ f df (Lift b)       (Forward c dc) = Forward a $ dc * runId (snd (df (Id a) (Id b) (Id c))) where a = f b c-    lift2_ f df (Forward b db) Zero           = Forward a $ runId (fst (df (Id a) (Id b) (Id 0))) * db where a = f b 0-    lift2_ f df (Forward b db) (Lift c)       = Forward a $ runId (fst (df (Id a) (Id b) (Id c))) * db where a = f b c-    lift2_ f df (Forward b db) (Forward c dc) = Forward a da-        where-            a = f b c-            (Id dadb, Id dadc) = df (Id a) (Id b) (Id c)-            da = dadb * db + dc * dadc--deriveLifted id $ conT ''Forward--bind :: (Traversable f, Num a) => (f (AD Forward a) -> b) -> f a -> f b-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 AD Zero)--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 AD Zero)-        b0 = f (lift <$> 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)--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)-        dropIx ((_,b),bs) = (b,bs)---- we can't transpose arbitrary traversables, since we can't construct one out of whole cloth, and the outer--- traversable could be empty. So instead we use one as a 'skeleton'-transposeWith :: (Functor f, Foldable f, Traversable g) => (b -> f a -> c) -> f (g a) -> g b -> g c-transposeWith f as = snd . mapAccumL go xss0-    where-        go xss b = (tail <$> xss, f b (head <$> xss))-        xss0 = toList <$> as-

@@ -1,139 +0,0 @@-{-# LANGUAGE GeneralizedNewtypeDeriving, FlexibleContexts, FlexibleInstances, MultiParamTypeClasses, DeriveDataTypeable #-}-{-# OPTIONS_HADDOCK hide #-}--------------------------------------------------------------------------------- |--- Module      :  Numeric.AD.Internal.Identity--- Copyright   :  (c) Edward Kmett 2010--- License     :  BSD3--- Maintainer  :  ekmett@gmail.com--- Stability   :  experimental--- Portability :  GHC only----------------------------------------------------------------------------------module Numeric.AD.Internal.Identity-    ( Id(..)-    , probe-    , unprobe-    , probed-    , unprobed-    ) where--import Control.Applicative-import Numeric.AD.Internal.Classes-import Numeric.AD.Internal.Types-import Data.Monoid-import Data.Data (Data)-import Data.Typeable (Typeable)-import Data.Traversable (Traversable, traverse)-import Data.Foldable (Foldable, foldMap)--newtype Id a = Id { runId :: a } deriving-    (Iso a, Eq, Ord, Show, Enum, Bounded, Num, Real, Fractional, Floating, RealFrac, RealFloat, Monoid, Data, Typeable)--probe :: a -> AD Id a-probe a = AD (Id a)--unprobe :: AD Id a -> a-unprobe (AD (Id a)) = a--pid :: f a -> f (Id a)-pid = iso--unpid :: f (Id a) -> f a-unpid = osi--probed :: f a -> f (AD Id a)-probed = iso . pid--unprobed :: f (AD Id a) -> f a-unprobed = unpid . osi--instance Functor Id where-    fmap f (Id a) = Id (f a)--instance Foldable Id where-    foldMap f (Id a) = f a--instance Traversable Id where-    traverse f (Id a) = Id <$> f a--instance Applicative Id where-    pure = Id-    Id f <*> Id a = Id (f a)--instance Monad Id where-    return = Id-    Id a >>= f = f a--instance Lifted Id where-    (==!) = (==)-    compare1 = compare-    showsPrec1 = showsPrec-    fromInteger1 = fromInteger-    (+!) = (+)-    (-!) = (-)-    (*!) = (*)-    negate1 = negate-    abs1 = abs-    signum1 = signum-    (/!) = (/)-    recip1 = recip-    fromRational1 = fromRational-    toRational1 = toRational-    pi1 = pi-    exp1 = exp-    log1 = log-    sqrt1 = sqrt-    (**!) = (**)-    logBase1 = logBase-    sin1 = sin-    cos1 = cos-    tan1 = tan-    asin1 = asin-    acos1 = acos-    atan1 = atan-    sinh1 = sinh-    cosh1 = cosh-    tanh1 = tanh-    asinh1 = asinh-    acosh1 = acosh-    atanh1 = atanh-    properFraction1 = properFraction-    truncate1 = truncate-    round1 = round-    ceiling1 = ceiling-    floor1 = floor-    floatRadix1 = floatRadix-    floatDigits1 = floatDigits-    floatRange1 = floatRange-    decodeFloat1 = decodeFloat-    encodeFloat1 = encodeFloat-    exponent1 = exponent-    significand1 = significand-    scaleFloat1 = scaleFloat-    isNaN1 = isNaN-    isInfinite1 = isInfinite-    isDenormalized1 = isDenormalized-    isNegativeZero1 = isNegativeZero-    isIEEE1 = isIEEE-    atan21 = atan2-    succ1 = succ-    pred1 = pred-    toEnum1 = toEnum-    fromEnum1 = fromEnum-    enumFrom1 = enumFrom-    enumFromThen1 = enumFromThen-    enumFromTo1 = enumFromTo-    enumFromThenTo1 = enumFromThenTo-    minBound1 = minBound-    maxBound1 = maxBound--instance Mode Id where-    lift = Id-    Id a ^* b = Id (a * b)-    a *^ Id b = Id (a * b)-    Id a <+> Id b = Id (a + b)-    Id a <**> Id b = Id (a ** b)--instance Primal Id where-    primal (Id a) = a

@@ -1,280 +0,0 @@-{-# LANGUAGE Rank2Types, TypeFamilies, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, FlexibleContexts, TemplateHaskell, UndecidableInstances, DeriveDataTypeable #-}--- {-# OPTIONS_HADDOCK hide, prune #-}--------------------------------------------------------------------------------- |--- Module      :  Numeric.AD.Internal.Reverse--- Copyright   :  (c) Edward Kmett 2010--- License     :  BSD3--- Maintainer  :  ekmett@gmail.com--- Stability   :  experimental--- Portability :  GHC only------ Reverse-Mode Automatic Differentiation implementation details------ 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.Internal.Reverse-    ( Reverse(..)-    , Tape(..)-    , partials-    , partialArray-    , partialMap-    , derivative-    , derivative'-    , Var(..)-    , bind-    , unbind-    , unbindMap-    , unbindWith-    , unbindMapWithDefault-    , vgrad, vgrad'-    , Grad(..)-    ) where--import Prelude hiding (mapM)-import Control.Applicative (Applicative(..),(<$>))-import Control.Monad.ST-import Control.Monad (forM_)-import Data.List (foldl', delete)-import Data.Array.ST-import Data.Array-import Data.IntMap (IntMap, fromListWith, findWithDefault, fromAscList, -                    updateLookupWithKey)-import qualified Data.IntSet as IS-import Data.Graph (graphFromEdges', Vertex, vertices, edges, transposeG, Graph)-import Data.Reify (reifyGraph, MuRef(..))-import qualified Data.Reify.Graph as Reified-import Data.Traversable (Traversable, mapM)-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---- | 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)---- | @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)---- deriving instance (Data (Tape a (Reverse a)) => Data (Reverse a)--instance MuRef (Reverse a) where-    type DeRef (Reverse a) = Tape a--    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--instance Lifted Reverse => Mode Reverse where-    lift 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--    _ <**> Reverse Zero     = lift 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--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--instance Lifted Reverse => Jacobian Reverse where-    type D Reverse = Id--    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)--    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 _         _         (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)--    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 ''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 = go (fromAscList . assocs $ transposeG g) starters-  where starters = IS.toList $ foldl' (flip IS.delete)-                                      (IS.fromList $ vertices g)-                                      (map snd $ edges g)-        go _ [] = []-        go g' (n:ns) = let (g'',ns') = foldl' (uncurry (prune n)) (g',[]) (g!n)-                       in n : go g'' (ns'++ns)-        prune n g' acc m = let f _ = Just . delete n-                               (Just ns, g'') = updateLookupWithKey f m g'-                           in g'' `seq` (g'', if null (tail ns) then m:acc else acc)----- | This returns a list of contributions to the partials.--- The variable ids returned in the list are likely /not/ unique!-partials :: Num a => AD Reverse a -> [(Int, a)]-partials (AD tape) = [ (ident, sensitivities ! ix) | (ix, Var _ ident) <- xs ]-    where-        Reified.Graph xs start = unsafePerformIO $ reifyGraph tape-        (g, vmap) = graphFromEdges' (edgeSet <$> filter nonConst xs)-        sensitivities = runSTArray $ do-            ss <- newArray (sbounds xs) 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!-        edgeSet (i, t) = (t, i, successors t)-        nonConst (_, Lift{}) = False-        nonConst _ = True-        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 Reverse a -> Array Int a-partialArray vbounds tape = accumArray (+) 0 vbounds (partials tape)-{-# INLINE partialArray #-}---- | Return an 'IntMap' of sparse partials-partialMap :: Num a => AD Reverse 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')---- | Used to mark variables for inspection during the reverse pass-class Primal v => Var v where-    var   :: a -> Int -> v a-    varId :: v a -> Int--instance Var Reverse where-    var a v = Reverse (Var a v)-    varId (Reverse (Var _ v)) = v-    varId _ = error "varId: not a Var"--instance Var (AD Reverse) where-    var a v = AD (var a v)-    varId (AD v) = varId v--bind :: (Traversable f, Var v) => f a -> (f (v a), (Int,Int))-bind xs = (r,(0,hi))-    where-        (r,hi) = runS (mapM freshVar xs) 0-        freshVar a = S (\s -> let s' = s + 1 in s' `seq` (var a s, s'))--unbind :: (Functor f, Var v)  => f (v a) -> Array Int a -> f a-unbind xs ys = fmap (\v -> ys ! varId v) xs--unbindWith :: (Functor f, Var v, Num a) => (a -> b -> c) -> f (v a) -> Array Int b -> f c-unbindWith f xs ys = fmap (\v -> f (primal v) (ys ! varId v)) xs--unbindMap :: (Functor f, Var v, Num a) => f (v a) -> IntMap a -> f a-unbindMap xs ys = fmap (\v -> findWithDefault 0 (varId v) ys) xs--unbindMapWithDefault :: (Functor f, Var v, Num a) => b -> (a -> b -> c) -> f (v a) -> IntMap b -> f c-unbindMapWithDefault z f xs ys = fmap (\v -> f (primal v) $ findWithDefault z (varId v) ys) xs--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-

@@ -1,255 +0,0 @@-{-# LANGUAGE BangPatterns, TemplateHaskell, TypeFamilies, TypeOperators, FlexibleContexts, UndecidableInstances, DeriveDataTypeable, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances #-}-{-# OPTIONS_GHC -fno-warn-name-shadowing #-}-module Numeric.AD.Internal.Sparse-    ( Index(..)-    , emptyIndex-    , addToIndex-    , indices-    , Sparse(..)-    , apply-    , vars-    , d, d', ds-    , skeleton-    , spartial-    , partial-    , vgrad-    , vgrad'-    , vgrads-    , Grad(..)-    , Grads(..)-    ) where--import Prelude hiding (lookup)-import Control.Applicative hiding ((<**>))-import Numeric.AD.Internal.Classes-import Control.Comonad.Cofree-import Numeric.AD.Internal.Types-import Data.Data-import Data.Typeable ()-import qualified Data.IntMap as IntMap-import Data.IntMap (IntMap, mapWithKey, unionWith, findWithDefault, toAscList, singleton, insertWith, lookup)-import Data.Traversable-import Language.Haskell.TH--newtype Index = Index (IntMap Int)--emptyIndex :: Index-emptyIndex = Index IntMap.empty-{-# INLINE emptyIndex #-}--addToIndex :: Int -> Index -> Index-addToIndex k (Index m) = Index (insertWith (+) k 1 m)-{-# INLINE addToIndex #-}--indices :: Index -> [Int]-indices (Index as) = uncurry (flip replicate) `concatMap` toAscList as-{-# INLINE indices #-}---- | We only store partials in sorted order, so the map contained in a partial--- will only contain partials with equal or greater keys to that of the map in--- which it was found. This should be key for efficiently computing sparse hessians.--- there are only (n + k - 1) choose k distinct nth partial derivatives of a--- function with k inputs.-data Sparse a-  = Sparse !a (IntMap (Sparse a))-  | Zero-  deriving (Show, Data, Typeable)---- | drop keys below a given value-dropMap :: Int -> IntMap a -> IntMap a-dropMap n = snd . IntMap.split (n - 1)-{-# INLINE dropMap #-}--times :: Num a => Sparse a -> Int -> Sparse a -> Sparse a-times Zero _ _ = Zero-times _ _ Zero = Zero-times (Sparse a as) n (Sparse b bs) = Sparse (a * b) $-    unionWith (<+>)-        (fmap (^* b) (dropMap n as))-        (fmap (a *^) (dropMap n bs))-{-# INLINE times #-}--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)-{-# INLINE vars #-}--apply :: (Traversable f, Num a) => (f (AD Sparse a) -> b) -> f a -> b-apply f = f . vars-{-# INLINE apply #-}--skeleton :: Traversable f => f a -> f Int-skeleton = snd . mapAccumL (\ !n _ -> (n + 1, n)) 0-{-# INLINE skeleton #-}--d :: (Traversable f, Num a) => f b -> AD Sparse a -> f a-d fs (AD Zero) = 0 <$ fs-d fs (AD (Sparse _ da)) = snd $ mapAccumL (\ !n _ -> (n + 1, maybe 0 primal $ lookup n da)) 0 fs-{-# INLINE d #-}--d' :: (Traversable f, Num a) => f a -> AD Sparse a -> (a, f a)-d' fs (AD Zero) = (0, 0 <$ fs)-d' fs (AD (Sparse a da)) = (a, snd $ mapAccumL (\ !n _ -> (n + 1, maybe 0 primal $ lookup n da)) 0 fs)-{-# INLINE d' #-}--ds :: (Traversable f, Num a) => f b -> AD Sparse a -> Cofree f a-ds fs (AD Zero) = r where r = 0 :< (r <$ fs)-ds fs (AD as@(Sparse a _)) = a :< (go emptyIndex <$> fns)-    where-        fns = skeleton fs-        -- go :: Index -> Int -> Cofree f a-        go ix i = partial (indices ix') as :< (go ix' <$> fns)-            where ix' = addToIndex i ix-{-# INLINE ds #-}--{--vvars :: Num a => Vector a -> Vector (AD Sparse a)-vvars = Vector.imap (\n a -> AD $ Sparse a $ singleton n $ lift 1)-{-# INLINE vvars #-}--vapply :: Num a => (Vector (AD Sparse a) -> b) -> Vector a -> b-vapply f = f . vvars-{-# INLINE vapply #-}---vd :: Num a => Int -> AD Sparse a -> Vector a-vd n (AD (Sparse _ da)) = Vector.generate n $ \i -> maybe 0 primal $ lookup i da-{-# INLINE vd #-}--vd' :: Num a => Int -> AD Sparse a -> (a, Vector a)-vd' n (AD (Sparse a da)) = (a , Vector.generate n $ \i -> maybe 0 primal $ lookup i da)-{-# INLINE vd' #-}--vds :: Num a => Int -> AD Sparse a -> Cofree Vector a-vds n (AD as@(Sparse a _)) = a :< Vector.generate n (go emptyIndex)-    where-        go ix i = partial (indices ix') as :< Vector.generate n (go ix')-            where ix' = addToIndex i ix-{-# INLINE vds #-}--}--partial :: Num a => [Int] -> Sparse a -> a-partial []     (Sparse a _)  = a-partial (n:ns) (Sparse _ da) = partial ns $ findWithDefault (lift 0) n da-partial _      Zero          = 0-{-# INLINE partial #-}--spartial :: Num a => [Int] -> Sparse a -> Maybe a-spartial [] (Sparse a _) = Just a-spartial (n:ns) (Sparse _ da) = do-    a' <- lookup n da-    spartial ns a'-spartial _  Zero         = Nothing-{-# INLINE spartial #-}--instance Primal Sparse where-    primal (Sparse a _) = a-    primal Zero = 0--instance Lifted Sparse => Mode Sparse where-    lift a = Sparse a IntMap.empty-    zero = Zero-    _ <**> Zero = lift 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-    Zero <+> a = a-    a <+> Zero = a-    Sparse a as <+> Sparse b bs = Sparse (a + b) $ unionWith (<+>) as bs-    Zero        ^* _ = Zero-    Sparse a as ^* b = Sparse (a * b) $ fmap (^* b) as-    _ *^ Zero        = Zero-    a *^ Sparse b bs = Sparse (a * b) $ fmap (a *^) bs-    Zero        ^/ _ = Zero-    Sparse a as ^/ b = Sparse (a / b) $ fmap (^/ b) as--instance Lifted Sparse => Jacobian Sparse where-    type D Sparse = Sparse-    unary f _ Zero = lift (f 0)-    unary f dadb (Sparse pb bs) = Sparse (f pb) $ mapWithKey (times dadb) bs--    lift1 f _ Zero = lift (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 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 _    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) $-        unionWith (<+>)-            (mapWithKey (times dadb) db)-            (mapWithKey (times dadc) dc)--    lift2 f _  Zero             Zero = lift (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-        (dadb, dadc) = df b c-        da = unionWith (<+>)-            (mapWithKey (times dadb) db)-            (mapWithKey (times dadc) dc)--    lift2_ f _  Zero             Zero = lift (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-        (dadb, dadc) = df a b c-        a = Sparse (f pb pc) da-        da = unionWith (<+>)-            (mapWithKey (times dadb) db)-            (mapWithKey (times dadc) dc)--deriveLifted id $ conT ''Sparse---class Num a => Grad i o o' a | i -> a o o', o -> a i o', o' -> a i o where-    pack :: i -> [AD Sparse a] -> AD Sparse a-    unpack :: ([a] -> [a]) -> o-    unpack' :: ([a] -> (a, [a])) -> o'--instance Num a => Grad (AD Sparse a) [a] (a, [a]) a where-    pack i _ = i-    unpack f = f []-    unpack' f = f []--instance Grad i o o' a => Grad (AD Sparse 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 = d as $ apply f as-{-# INLINE vgrad #-}--vgrad' :: Grad i o o' a => i -> o'-vgrad' i = unpack' (unsafeGrad' (pack i))-    where-        unsafeGrad' f as = d' as $ apply f as-{-# INLINE vgrad' #-}--class Num a => Grads i o a | i -> a o, o -> a i where-    packs :: i -> [AD Sparse a] -> AD Sparse a-    unpacks :: ([a] -> Cofree [] a) -> o--instance Num a => Grads (AD Sparse a) (Cofree [] a) a where-    packs i _ = i-    unpacks f = f []--instance Grads i o a => Grads (AD Sparse a -> i) (a -> o) a where-    packs f (a:as) = packs (f a) as-    packs _ [] = error "Grad.pack: logic error"-    unpacks f a = unpacks (f . (a:))--vgrads :: Grads i o a => i -> o-vgrads i = unpacks (unsafeGrads (packs i))-    where-        unsafeGrads f as = ds as $ apply f as-{-# INLINE vgrads #-}-

@@ -1,85 +0,0 @@-{-# LANGUAGE TypeOperators, TemplateHaskell, ScopedTypeVariables, FlexibleContexts #-}-{-# OPTIONS_HADDOCK hide #-}--------------------------------------------------------------------------------- |--- Module      :  Numeric.AD.Internal.Tensors--- Copyright   :  (c) Edward Kmett 2010--- License     :  BSD3--- Maintainer  :  ekmett@gmail.com--- Stability   :  experimental--- Portability :  GHC only-----------------------------------------------------------------------------------module Numeric.AD.Internal.Tensors-    ( Tensors(..)-    , headT-    , tailT-    , tensors-    ) where--import Control.Applicative-import Data.Foldable-import Data.Traversable-import Data.Monoid-#if __GLASGOW_HASKELL__ < 704-import Data.Typeable (Typeable1(..), TyCon, mkTyCon, mkTyConApp)-#else-import Data.Typeable (Typeable1(..), TyCon, mkTyCon3, mkTyConApp)-#endif-import Control.Comonad.Cofree--infixl 3 :---data Tensors f a = a :- Tensors f (f a)--newtype Showable = Showable (Int -> String -> String)--instance Show Showable where-  showsPrec d (Showable f) = f d--showable :: Show a => a -> Showable-showable a = Showable (\d -> showsPrec d a)---- Polymorphic recursion precludes 'Data' in its current form, as no Data1 class exists--- Polymorphic recursion also breaks 'show' for 'Tensors'!--- factor Show1 out of Lifted?-instance (Functor f, Show (f Showable), Show a) => Show (Tensors f a) where-  showsPrec d (a :- as) = showParen (d > 3) $ -    showsPrec 4 a . showString " :- " . showsPrec 3 (fmap showable <$> as)--instance Functor f => Functor (Tensors f) where-    fmap f (a :- as) = f a :- fmap (fmap f) as--instance Foldable f => Foldable (Tensors f) where-    foldMap f (a :- as) = f a `mappend` foldMap (foldMap f) as--instance Traversable f => Traversable (Tensors f) where-    traverse f (a :- as) = (:-) <$> f a <*> traverse (traverse f) as--tailT :: Tensors f a -> Tensors f (f a)-tailT (_ :- as) = as-{-# INLINE tailT #-}--headT :: Tensors f a -> a-headT (a :- _) = a-{-# INLINE headT #-}--tensors :: Functor f => Cofree f a -> Tensors f a-tensors (a :< as) = a :- dist (tensors <$> as)-    where-        dist :: Functor f => f (Tensors f a) -> Tensors f (f a)-        dist x = (headT <$> x) :- dist (tailT <$> x)--instance Typeable1 f => Typeable1 (Tensors f) where-    typeOf1 tfa = mkTyConApp tensorsTyCon [typeOf1 (undefined `asArgsType` tfa)]-        where asArgsType :: f a -> t f a -> f a-              asArgsType = const--tensorsTyCon :: TyCon-#if __GLASGOW_HASKELL__ < 704-tensorsTyCon = mkTyCon "Numeric.AD.Internal.Tensors.Tensors"-#else-tensorsTyCon = mkTyCon3 "ad" "Numeric.AD.Internal.Tensors" "Tensors"-#endif-{-# NOINLINE tensorsTyCon #-}

@@ -1,139 +0,0 @@-{-# LANGUAGE Rank2Types, TypeFamilies, FlexibleContexts, UndecidableInstances, TemplateHaskell, DeriveDataTypeable #-}-{-# OPTIONS_GHC -fno-warn-name-shadowing #-}--- {-# OPTIONS_HADDOCK hide, prune #-}--------------------------------------------------------------------------------- |--- Module      : Numeric.AD.Tower.Internal--- Copyright   : (c) Edward Kmett 2010--- License     : BSD3--- Maintainer  : ekmett@gmail.com--- Stability   : experimental--- Portability : GHC only-----------------------------------------------------------------------------------module Numeric.AD.Internal.Tower-    ( Tower(..)-    , zeroPad-    , zeroPadF-    , transposePadF-    , d-    , d'-    , withD-    , tangents-    , bundle-    , apply-    , getADTower-    , tower-    ) where--import Prelude hiding (all)-import Control.Applicative hiding ((<**>))-import Data.Foldable-import Data.Data (Data)-import Data.Typeable (Typeable)-import Language.Haskell.TH-import Numeric.AD.Internal.Types-import Numeric.AD.Internal.Classes---- | @Tower@ is an AD 'Mode' that calculates a tangent tower by forward AD, and provides fast 'diffsUU', 'diffsUF'-newtype Tower a = Tower { getTower :: [a] } deriving (Data, Typeable)--instance Show a => Show (Tower a) where-    showsPrec n (Tower as) = showParen (n > 10) $ showString "Tower " . showList as---- Local combinators--zeroPad :: Num a => [a] -> [a]-zeroPad xs = xs ++ repeat 0-{-# INLINE zeroPad #-}--zeroPadF :: (Functor f, Num a) => [f a] -> [f a]-zeroPadF fxs@(fx:_) = fxs ++ repeat (const 0 <$> fx)-zeroPadF _ = error "zeroPadF :: empty list"-{-# INLINE zeroPadF #-}--transposePadF :: (Foldable f, Functor f) => a -> f [a] -> [f a]-transposePadF pad fx-    | all null fx = []-    | otherwise = fmap headPad fx : transposePadF pad (drop1 <$> fx)-    where-        headPad [] = pad-        headPad (x:_) = x-        drop1 (_:xs) = xs-        drop1 xs = xs--d :: Num a => [a] -> a-d (_:da:_) = da-d _ = 0-{-# INLINE d #-}--d' :: Num a => [a] -> (a, a)-d' (a:da:_) = (a, da)-d' (a:_)    = (a, 0)-d' _        = (0, 0)-{-# INLINE d' #-}--tangents :: Tower a -> Tower a-tangents (Tower []) = Tower []-tangents (Tower (_:xs)) = Tower xs-{-# INLINE tangents #-}--bundle :: a -> Tower a -> Tower a-bundle a (Tower as) = Tower (a:as)-{-# INLINE bundle #-}--withD :: (a, a) -> AD Tower a-withD (a, da) = AD (Tower [a,da])-{-# INLINE withD #-}--apply :: Num a => (AD Tower a -> b) -> a -> b-apply f a = f (AD (Tower [a,1]))-{-# INLINE apply #-}--getADTower :: AD Tower a -> [a]-getADTower (AD t) = getTower t-{-# INLINE getADTower #-}--tower :: [a] -> AD Tower a-tower as = AD (Tower as)--instance Primal Tower where-    primal (Tower (x:_)) = x-    primal _ = 0--instance Lifted Tower => Mode Tower where-    lift a = Tower [a]-    zero = Tower []-    _ <**> Tower []  = lift 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--    Tower [] <+> bs = bs-    as <+> Tower [] = as-    Tower (a:as) <+> Tower (b:bs) = Tower (c:cs)-        where-            c = a + b-            Tower cs = Tower as <+> Tower bs--    a *^ Tower bs = Tower (map (a*) bs)-    Tower as ^* b = Tower (map (*b) as)-    Tower as ^/ b = Tower (map (/b) as)--instance Lifted Tower => Jacobian Tower where-    type D Tower = Tower-    unary f dadb b = bundle (f (primal b)) (tangents b *! dadb)-    lift1 f df b   = bundle (f (primal b)) (tangents b *! df b)-    lift1_ f df b = a where-        a = bundle (f (primal b)) (tangents b *! df a b)--    binary f dadb dadc b c = bundle (f (primal b) (primal c)) (tangents b *! dadb +! tangents c *! dadc)-    lift2 f df b c = bundle (f (primal b) (primal c)) (tangents b *! dadb +! tangents c *! dadc) where-        (dadb, dadc) = df b c-    lift2_ f df b c = a where-        a0 = f (primal b) (primal c)-        da = tangents b *! dadb +! tangents c *! dadc-        a = bundle a0 da-        (dadb, dadc) = df a b c--deriveLifted id (conT ''Tower)

@@ -1,69 +0,0 @@-{-# LANGUAGE Rank2Types, GeneralizedNewtypeDeriving, TemplateHaskell, FlexibleContexts, FlexibleInstances, MultiParamTypeClasses, UndecidableInstances #-}-{-# OPTIONS_HADDOCK hide #-}--------------------------------------------------------------------------------- |--- Module      :  Numeric.AD.Internal.Types--- Copyright   :  (c) Edward Kmett 2010--- License     :  BSD3--- Maintainer  :  ekmett@gmail.com--- Stability   :  experimental--- Portability :  GHC only----------------------------------------------------------------------------------module Numeric.AD.Internal.Types-    ( AD(..)-    , UU, UF, FU, FF-    ) where--import Data.Data (Data(..), mkDataType, DataType, mkConstr, Constr, constrIndex, Fixity(..))-import Data.Typeable (Typeable1(..), Typeable(..), TyCon, mkTyCon3, mkTyConApp, gcast1)-import Language.Haskell.TH-import Numeric.AD.Internal.Classes---- | 'AD' serves as a common wrapper for different 'Mode' instances, exposing a traditional--- numerical tower. Universal quantification is used to limit the actions in user code to--- machinery that will return the same answers under all AD modes, allowing us to use modes--- interchangeably as both the type level \"brand\" and dictionary, providing a common API.-newtype AD f a = AD { runAD :: f a } deriving (Iso (f a), Lifted, Mode, Primal)---- > instance (Lifted f, Num a) => Num (AD f a)--- etc.-let f = varT (mkName "f") in -    deriveNumeric -        (classP ''Lifted [f]:) -        (conT ''AD `appT` f)---- | A scalar-to-scalar automatically-differentiable function.-type UU a = forall s. Mode s => AD s a -> AD s a--- | A scalar-to-non-scalar automatically-differentiable function.-type UF f a = forall s. Mode s => AD s a -> f (AD s a)--- | A non-scalar-to-scalar automatically-differentiable function.-type FU f a = forall s. Mode s => f (AD s a) -> AD s a--- | A non-scalar-to-non-scalar automatically-differentiable function.-type FF f g a = forall s. Mode s => f (AD s a) -> g (AD s a)--instance Typeable1 f => Typeable1 (AD f) where-    typeOf1 tfa = mkTyConApp adTyCon [typeOf1 (undefined `asArgsType` tfa)]-        where asArgsType :: f a -> t f a -> f a-              asArgsType = const--adTyCon :: TyCon-adTyCon = mkTyCon3 "ad" "Numeric.AD.Internal.Types" "AD"-{-# NOINLINE adTyCon #-}--adConstr :: Constr-adConstr = mkConstr adDataType "AD" [] Prefix-{-# NOINLINE adConstr #-}--adDataType :: DataType-adDataType = mkDataType "Numeric.AD.Internal.Types.AD" [adConstr]-{-# NOINLINE adDataType #-}--instance (Typeable1 f, Typeable a, Data (f a), Data a) => Data (AD f a) where-    gfoldl f z (AD a) = z AD `f` a-    toConstr _ = adConstr-    gunfold k z c = case constrIndex c of-        1 -> k (z AD)-        _ -> error "gunfold"-    dataTypeOf _ = adDataType-    dataCast1 f = gcast1 f

@@ -1,93 +0,0 @@-{-# LANGUAGE Rank2Types #-}--------------------------------------------------------------------------------- |--- Module      :  Numeric.AD.Mode.Directed--- Copyright   :  (c) Edward Kmett 2010--- License     :  BSD3--- Maintainer  :  ekmett@gmail.com--- Stability   :  experimental--- Portability :  GHC only------ Allows the choice of AD 'Mode' to be specified at the term level for--- benchmarking or more complicated usage patterns.--------------------------------------------------------------------------------module Numeric.AD.Mode.Directed-    (-    -- * Gradients-      grad-    , grad'-    -- * Jacobians-    , jacobian-    , jacobian'-    -- * Derivatives-    , diff-    , diff'-    -- * Exposed Types-    , UU, UF, FU, FF-    , Direction(..)-    , Mode(..)-    , AD(..)-    ) where--import Prelude hiding (reverse)-import Numeric.AD.Types-import Numeric.AD.Classes-import Data.Traversable (Traversable)-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.Mixed 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-    | Reverse-    | Tower-    | Mixed-    deriving (Show, Eq, Ord, Read, Bounded, Enum, Ix)--diff :: Num a => Direction -> UU a -> a -> a-diff Forward = F.diff-diff Reverse = R.diff-diff Tower = T.diff-diff Mixed = F.diff-{-# INLINE diff #-}--diff' :: Num a => Direction -> UU a -> a -> (a, a)-diff' Forward = F.diff'-diff' Reverse = R.diff'-diff' Tower = T.diff'-diff' Mixed = F.diff'-{-# INLINE diff' #-}--jacobian :: (Traversable f, Traversable g, Num a) => Direction -> FF f g a -> f a -> g (f a)-jacobian Forward = F.jacobian-jacobian Reverse = R.jacobian-jacobian Tower = F.jacobian -- error "jacobian Tower: unimplemented"-jacobian Mixed = M.jacobian-{-# INLINE jacobian #-}--jacobian' :: (Traversable f, Traversable g, Num a) => Direction -> FF f g a -> f a -> g (a, f a)-jacobian' Forward = F.jacobian'-jacobian' Reverse = R.jacobian'-jacobian' Tower = F.jacobian' -- error "jacobian' Tower: unimplemented"-jacobian' Mixed = M.jacobian'-{-# INLINE jacobian' #-}--grad :: (Traversable f, Num a) => Direction -> FU f a -> f a -> f a-grad Forward = F.grad-grad Reverse = R.grad-grad Tower   = F.grad -- error "grad Tower: unimplemented"-grad Mixed   = M.grad-{-# INLINE grad #-}--grad' :: (Traversable f, Num a) => Direction -> FU f a -> f a -> (a, f a)-grad' Forward = F.grad'-grad' Reverse = R.grad'-grad' Tower   = F.grad' -- error "grad' Tower: unimplemented"-grad' Mixed   = M.grad'-{-# INLINE grad' #-}-

@@ -1,165 +0,0 @@-{-# LANGUAGE Rank2Types #-}--------------------------------------------------------------------------------- |--- Module      :  Numeric.AD.Mode.Forward--- Copyright   :  (c) Edward Kmett 2010--- License     :  BSD3--- Maintainer  :  ekmett@gmail.com--- Stability   :  experimental--- Portability :  GHC only------ Forward mode automatic differentiation-----------------------------------------------------------------------------------module Numeric.AD.Mode.Forward-    (-    -- * Gradient-      grad-    , grad'-    , gradWith-    , gradWith'-    -- * Jacobian-    , jacobian-    , jacobian'-    , jacobianWith-    , jacobianWith'-    -- * Transposed Jacobian-    , jacobianT-    , jacobianWithT-    -- * Hessian Product-    , hessianProduct-    , hessianProduct'-    -- * Derivatives-    , diff-    , diff'-    , diffF-    , diffF'-    -- * Directional Derivatives-    , du-    , du'-    , duF-    , duF'-    -- * Exposed Types-    , UU, UF, FU, FF-    , AD(..)-    , Mode(..)-    ) where--import Data.Traversable (Traversable)-import Control.Applicative-import Numeric.AD.Types-import Numeric.AD.Internal.Classes-import Numeric.AD.Internal.Composition-import Numeric.AD.Internal.Forward--du :: (Functor f, Num a) => FU f a -> f (a, a) -> a-du f = tangent . f . fmap (uncurry bundle)-{-# INLINE du #-}--du' :: (Functor f, Num a) => FU f a -> f (a, a) -> (a, a)-du' f = unbundle . f . fmap (uncurry bundle)-{-# INLINE du' #-}--duF :: (Functor f, Functor g, Num a) => FF f g a -> f (a, a) -> g a-duF f = fmap tangent . f . fmap (uncurry bundle)-{-# INLINE duF #-}--duF' :: (Functor f, Functor g, Num a) => FF f g a -> f (a, a) -> g (a, a)-duF' f = fmap unbundle . f . fmap (uncurry bundle)-{-# INLINE duF' #-}---- | The 'diff' function calculates the first derivative of a scalar-to-scalar function by forward-mode 'AD'------ > diff sin == cos-diff :: Num a => UU a -> a -> a-diff f a = tangent $ apply f a-{-# INLINE diff #-}---- | The 'd'UU' function calculates the result and first derivative of scalar-to-scalar function by F'orward' 'AD'--- --- > d' sin == sin &&& cos--- > d' f = f &&& d f-diff' :: Num a => UU a -> a -> (a, a)-diff' f a = unbundle $ apply f a-{-# INLINE diff' #-}---- | The 'diffF' function calculates the first derivative of scalar-to-nonscalar function by F'orward' 'AD'-diffF :: (Functor f, Num a) => UF f a -> a -> f a-diffF f a = tangent <$> apply f a-{-# INLINE diffF #-}---- | The 'diffF'' function calculates the result and first derivative of a scalar-to-non-scalar function by F'orward' 'AD'-diffF' :: (Functor f, Num a) => UF f a -> a -> f (a, a)-diffF' f a = unbundle <$> apply f a-{-# INLINE diffF' #-}---- | A fast, simple transposed Jacobian computed with forward-mode AD.-jacobianT :: (Traversable f, Functor g, Num a) => FF f g a -> f a -> f (g a)-jacobianT f = bind (fmap tangent . f)-{-# INLINE jacobianT #-}---- | A fast, simple transposed Jacobian computed with forward-mode AD.-jacobianWithT :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> FF f g a -> f a -> f (g b)-jacobianWithT g f = bindWith g' f-    where g' a ga = g a . tangent <$> ga-{-# INLINE jacobianWithT #-}--jacobian :: (Traversable f, Traversable g, Num a) => FF f g a -> f a -> g (f a)-jacobian f as = transposeWith (const id) t p-    where-        (p, t) = bind' (fmap tangent . f) as-{-# INLINE jacobian #-}--jacobianWith :: (Traversable f, Traversable g, Num a) => (a -> a -> b) -> FF f g a -> f a -> g (f b)-jacobianWith g f as = transposeWith (const id) t p-    where-        (p, t) = bindWith' g' f as-        g' a ga = g a . tangent <$> ga-{-# INLINE jacobianWith #-}--jacobian' :: (Traversable f, Traversable g, Num a) => FF f g a -> f a -> g (a, f a)-jacobian' f as = transposeWith row t p-    where-        (p, t) = bind' f as-        row x as' = (primal x, tangent <$> as')-{-# INLINE jacobian' #-}--jacobianWith' :: (Traversable f, Traversable g, Num a) => (a -> a -> b) -> FF f g a -> f a -> g (a, f b)-jacobianWith' g f as = transposeWith row t p-    where-        (p, t) = bindWith' g' f as-        row x as' = (primal x, as')-        g' a ga = g a . tangent <$> ga-{-# INLINE jacobianWith' #-}--grad :: (Traversable f, Num a) => FU f a -> f a -> f a-grad f = bind (tangent . f)-{-# INLINE grad #-}--grad' :: (Traversable f, Num a) => FU f a -> f a -> (a, f a)-grad' f as = (primal b, tangent <$> bs)-    where-        (b, bs) = bind' f as-{-# INLINE grad' #-}--gradWith :: (Traversable f, Num a) => (a -> a -> b) -> FU f a -> f a -> f b-gradWith g f = bindWith g (tangent . f)-{-# INLINE gradWith #-}--gradWith' :: (Traversable f, Num a) => (a -> a -> b) -> FU f a -> f a -> (a, f b)-gradWith' g f = bindWith' g (tangent . f)-{-# INLINE gradWith' #-}---- | Compute the product of a vector with the Hessian using forward-on-forward-mode AD. -hessianProduct :: (Traversable f, Num a) => FU f a -> f (a, a) -> f a-hessianProduct f = duF $ grad $ decomposeMode . f . fmap composeMode---- | Compute the gradient and hessian product using forward-on-forward-mode AD. -hessianProduct' :: (Traversable f, Num a) => FU f a -> f (a, a) -> f (a, a)-hessianProduct' f = duF' $ grad $ decomposeMode . f . fmap composeMode---- * Experimental---- data f :> a = a :< f (f :> a)--- gradients :: (Traversable f, Num a) => FU f a -> f a -> (f :> a)

@@ -1,226 +0,0 @@-{-# LANGUAGE Rank2Types, TypeFamilies, PatternGuards #-}--------------------------------------------------------------------------------- |--- Module      :  Numeric.AD.Mode.Mixed--- Copyright   :  (c) Edward Kmett 2010--- License     :  BSD3--- Maintainer  :  ekmett@gmail.com--- Stability   :  experimental--- Portability :  GHC only------ Mixed-Mode Automatic Differentiation.------ Each combinator exported from this module chooses an appropriate AD mode.--- The following basic operations are supported, modified as appropriate by the suffixes below:--- --- * 'grad' computes the gradient (partial derivatives) of a function at a point------ * 'jacobian' computes the Jacobian matrix of a function at a point------ * 'diff' computes the derivative of a function at a point------ * 'du' computes a directional derivative of a function at a point--- --- * 'hessian' compute the Hessian matrix (matrix of second partial derivatives) of a function at a point--- --- The suffixes have the following meanings:--- --- * @\'@ -- also return the answer------ * @With@ lets the user supply a function to blend the input with the output------ * @F@ is a version of the base function lifted to return a 'Traversable' (or 'Functor') result------ * @s@ means the function returns all higher derivatives in a list or f-branching 'Stream'------ * @T@ means the result is transposed with respect to the traditional formulation.------ * @0@ means that the resulting derivative list is padded with 0s at the end.--------------------------------------------------------------------------------module Numeric.AD.Mode.Mixed-    (-    -- * Gradients (Reverse Mode)-      grad-    , grad'-    , gradWith-    , gradWith'--    -- * Higher Order Gradients (Sparse-on-Reverse)-    , grads--    -- * Jacobians (Sparse or Reverse)-    , jacobian-    , jacobian'-    , jacobianWith-    , jacobianWith'--    -- * Higher Order Jacobian (Sparse-on-Reverse)-    , jacobians--    -- * Transposed Jacobians (Forward Mode)-    , jacobianT-    , jacobianWithT--    -- * Hessian (Sparse-On-Reverse)-    , hessian-    , hessian'--    -- * Hessian Tensors (Sparse or Sparse-On-Reverse)-    , hessianF-    -- * Hessian Tensors (Sparse)-    , hessianF'--    -- * Hessian Vector Products (Forward-On-Reverse)-    , hessianProduct-    , hessianProduct'--    -- * Derivatives (Forward Mode)-    , diff-    , diffF--    , diff'-    , diffF'--    -- * Derivatives (Tower)-    , diffs-    , diffsF--    , diffs0-    , diffs0F--    -- * Directional Derivatives (Forward Mode)-    , du-    , du'-    , duF-    , duF'--    -- * Directional Derivatives (Tower)-    , dus-    , dus0-    , dusF-    , dus0F--    -- * Taylor Series (Tower)-    , taylor-    , taylor0--    -- * Maclaurin Series (Tower)-    , maclaurin-    , maclaurin0--    -- * Unsafe Variadic Grad-    , vgrad-    , vgrad'-    , vgrads--    -- * Exposed Types-    , module Numeric.AD.Types-    , Mode(..)-    , Grad-    , Grads-    ) where--import Data.Traversable (Traversable)-import Data.Foldable (Foldable, foldr')-import Control.Applicative--import Numeric.AD.Types-import Numeric.AD.Internal.Composition-import Numeric.AD.Classes (Mode(..))--import Numeric.AD.Mode.Forward -    ( diff, diff', diffF, diffF'-    , du, du', duF, duF'-    , jacobianT, jacobianWithT ) --import Numeric.AD.Mode.Tower -    ( diffsF, diffs0F, diffs, diffs0-    , taylor, taylor0, maclaurin, maclaurin0-    , dus, dus0, dusF, dus0F )--import qualified Numeric.AD.Mode.Reverse as Reverse-import Numeric.AD.Mode.Reverse -    ( grad, grad', gradWith, gradWith', vgrad, vgrad', Grad)---- temporary until we make a full sparse mode-import qualified Numeric.AD.Mode.Sparse as Sparse-import Numeric.AD.Mode.Sparse-    ( grads, jacobians, hessian', hessianF', vgrads, Grads)-    --- | Calculate the Jacobian of a non-scalar-to-non-scalar function, automatically choosing between forward and reverse mode AD based on the number of inputs and outputs.------ If you know the relative number of inputs and outputs, consider 'Numeric.AD.Reverse.jacobian' or 'Nuneric.AD.Sparse.jacobian'.-jacobian :: (Traversable f, Functor g, Num a) => FF f g a -> f a -> g (f a)-jacobian f bs = snd <$> jacobian' f bs-{-# INLINE jacobian #-}--data Nat = Z | S Nat deriving (Eq, Ord)--size :: Foldable f => f a -> Nat-size = foldr' (\_ b -> S b) Z --big :: Nat -> Bool-big (S (S (S (S (S (S (S (S (S (S _)))))))))) = True-big _ = False---- | Calculate both the answer and Jacobian of a non-scalar-to-non-scalar function, automatically choosing between forward- and reverse- mode AD based on the relative, based on the number of inputs------ If you know the relative number of inputs and outputs, consider 'Numeric.AD.Reverse.jacobian'' or 'Nuneric.AD.Sparse.jacobian''.-jacobian' :: (Traversable f, Functor g, Num a) => FF f g a -> f a -> g (a, f a)-jacobian' f bs | Z <- n = fmap (\x -> (unprobe x, bs)) (f (probed bs))-               | big n  = Reverse.jacobian' f bs-               | otherwise = Sparse.jacobian' f bs-    where-        n = size bs-{-# INLINE jacobian' #-}---- | @'jacobianWith' g f@ calculates the Jacobian of a non-scalar-to-non-scalar function, automatically choosing between forward and reverse mode AD based on the number of inputs and outputs.------ The resulting Jacobian matrix is then recombined element-wise with the input using @g@.------ If you know the relative number of inputs and outputs, consider 'Numeric.AD.Reverse.jacobianWith' or 'Nuneric.AD.Sparse.jacobianWith'.-jacobianWith :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> FF f g a -> f a -> g (f b)-jacobianWith g f bs = snd <$> jacobianWith' g f bs-{-# INLINE jacobianWith #-}---- | @'jacobianWith'' g f@ calculates the answer and Jacobian of a non-scalar-to-non-scalar function, automatically choosing between sparse and reverse mode AD based on the number of inputs and outputs.------ The resulting Jacobian matrix is then recombined element-wise with the input using @g@.------ If you know the relative number of inputs and outputs, consider 'Numeric.AD.Reverse.jacobianWith'' or 'Nuneric.AD.Sparse.jacobianWith''.-jacobianWith' :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> FF f g a -> f a -> g (a, f b)-jacobianWith' g f bs-    | Z <- n = fmap (\x -> (unprobe x, undefined <$> bs)) (f (probed bs))-    | big n  = Reverse.jacobianWith' g f bs-    | otherwise = Sparse.jacobianWith' g f bs-    where-        n = size bs-{-# INLINE jacobianWith' #-}---- | @'hessianProduct' f wv@ computes the product of the hessian @H@ of a non-scalar-to-scalar function @f@ at @w = 'fst' <$> wv@ with a vector @v = snd <$> wv@ using \"Pearlmutter\'s method\" from <http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.29.6143>, which states:------ > H v = (d/dr) grad_w (w + r v) | r = 0--- --- Or in other words, we take the directional derivative of the gradient. The gradient is calculated in reverse mode, then the directional derivative is calculated in forward mode.----hessianProduct :: (Traversable f, Num a) => FU f a -> f (a, a) -> f a-hessianProduct f = duF (grad (decomposeMode . f . fmap composeMode))---- | @'hessianProduct'' f wv@ computes both the gradient of a non-scalar-to-scalar @f@ at @w = 'fst' <$> wv@ and the product of the hessian @H@ at @w@ with a vector @v = snd <$> wv@ using \"Pearlmutter's method\". The outputs are returned wrapped in the same functor.------ > H v = (d/dr) grad_w (w + r v) | r = 0--- --- Or in other words, we return the gradient and the directional derivative of the gradient. The gradient is calculated in reverse mode, then the directional derivative is calculated in forward mode.-hessianProduct' :: (Traversable f, Num a) => FU f a -> f (a, a) -> f (a, a)-hessianProduct' f = duF' (grad (decomposeMode . f . fmap composeMode))---- | Compute the Hessian via the Jacobian of the gradient. gradient is computed in reverse mode and then the Jacobian is computed in sparse (forward) mode.-hessian :: (Traversable f, Num a) => FU f a -> f a -> f (f a)-hessian f = Sparse.jacobian (grad (decomposeMode . f . fmap composeMode))---- | Compute the order 3 Hessian tensor on a non-scalar-to-non-scalar function using Sparse or Sparse-on-Reverse -hessianF :: (Traversable f, Functor g, Num a) => FF f g a -> f a -> g (f (f a))-hessianF f as -    | big (size as) = decomposeFunctor $ Sparse.jacobian (ComposeFunctor . Reverse.jacobian (fmap decomposeMode . f . fmap composeMode)) as-    | otherwise = Sparse.hessianF f as

@@ -1,161 +0,0 @@--- {-# LANGUAGE Rank2Types, TemplateHaskell, BangPatterns #-}-{-# LANGUAGE Rank2Types, TemplateHaskell, BangPatterns, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, UndecidableInstances, ScopedTypeVariables #-}--------------------------------------------------------------------------------- |--- Module      :  Numeric.AD.Mode.Reverse--- Copyright   :  (c) Edward Kmett 2010--- License     :  BSD3--- Maintainer  :  ekmett@gmail.com--- 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.-----------------------------------------------------------------------------------module Numeric.AD.Mode.Reverse-    (-    -- * Gradient-      grad-    , grad'-    , gradWith-    , gradWith'--    -- * Jacobian-    , jacobian-    , jacobian'-    , jacobianWith-    , jacobianWith'-    -- * Hessian-    , hessian-    , hessianF-    -- * Derivatives-    , diff-    , diff'-    , diffF-    , diffF'-    -- * Unsafe Variadic Gradient-    , vgrad, vgrad'-    -- * Exposed Types-    , UU, UF, FU, FF-    , AD(..)-    , Mode(..)-    , 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.Reverse---- | The 'grad' function calculates the gradient of a non-scalar-to-scalar function with 'Reverse' AD in a single pass.-grad :: (Traversable f, Num a) => FU f 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' AD in a single pass.-grad' :: (Traversable f, Num a) => FU f 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) -> FU f 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' 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) -> FU f 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 :: (Traversable f, Functor g, Num a) => FF f g 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''-jacobian' :: (Traversable f, Functor g, Num a) => FF f g 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) -> FF f g 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) -> FF f g 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' #-}--diff :: Num a => UU a -> a -> a-diff f a = derivative $ f (var a 0)-{-# INLINE diff #-}---- | The 'd'' function calculates the value and derivative, as a--- pair, of a scalar-to-scalar function.-diff' :: Num a => UU a -> a -> (a, a)-diff' f a = derivative' $ f (var a 0)-{-# INLINE diff' #-}--diffF :: (Functor f, Num a) => UF f a -> a -> f a-diffF f a = derivative <$> f (var a 0)-{-# INLINE diffF #-}--diffF' :: (Functor f, Num a) => UF f 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 :: (Traversable f, Num a) => FU f 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 :: (Traversable f, Functor g, Num a) => FF f g a -> f a -> g (f (f a))-hessianF f = decomposeFunctor . jacobian (ComposeFunctor . jacobian (fmap decomposeMode . f . fmap composeMode))-

@@ -1,125 +0,0 @@-{-# LANGUAGE Rank2Types, BangPatterns #-}--------------------------------------------------------------------------------- |--- Module      : Numeric.AD.Mode.Sparse--- Copyright   : (c) Edward Kmett 2010--- License     : BSD3--- Maintainer  : ekmett@gmail.com--- Stability   : experimental--- Portability : GHC only------ Higher order derivatives via a \"dual number tower\".-----------------------------------------------------------------------------------module Numeric.AD.Mode.Sparse-    (-    -- * Sparse Gradients-      grad-    , grad'-    , gradWith-    , gradWith'-    , grads-    -    -- * Sparse Jacobians (synonyms)-    , jacobian-    , jacobian'-    , jacobianWith-    , jacobianWith'-    , jacobians--    -- * Sparse Hessians-    , hessian-    , hessian'--    , hessianF-    , hessianF'--    -- * Unsafe gradients-    , vgrad-    , vgrads--    -- * Exposed Types-    , module Numeric.AD.Types-    , Mode(..)-    , Grad-    , Grads-    ) where--import Control.Comonad-import Control.Applicative ((<$>))-import Data.Traversable-import Control.Comonad.Cofree-import Numeric.AD.Types-import Numeric.AD.Classes-import Numeric.AD.Internal.Sparse-import Numeric.AD.Internal.Combinators--second :: (a -> b) -> (c, a) -> (c, b)-second g (a,b) = (a, g b)-{-# INLINE second #-}--grad :: (Traversable f, Num a) => FU f a -> f a -> f a-grad f as = d as $ apply f as-{-# INLINE grad #-}--grad' :: (Traversable f, Num a) => FU f a -> f a -> (a, f a)-grad' f as = d' as $ apply f as-{-# INLINE grad' #-}--gradWith :: (Traversable f, Num a) => (a -> a -> b) -> FU f a -> f a -> f b-gradWith g f as = zipWithT g as $ grad f as-{-# INLINE gradWith #-}--gradWith' :: (Traversable f, Num a) => (a -> a -> b) -> FU f a -> f a -> (a, f b)-gradWith' g f as = second (zipWithT g as) $ grad' f as-{-# INLINE gradWith' #-}--jacobian :: (Traversable f, Functor g, Num a) => FF f g a -> f a -> g (f a)-jacobian f as = d as <$> apply f as-{-# INLINE jacobian #-}--jacobian' :: (Traversable f, Functor g, Num a) => FF f g a -> f a -> g (a, f a)-jacobian' f as = d' as <$> apply f as-{-# INLINE jacobian' #-}--jacobianWith :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> FF f g a -> f a -> g (f b)-jacobianWith g f as = zipWithT g as <$> jacobian f as-{-# INLINE jacobianWith #-}--jacobianWith' :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> FF f g a -> f a -> g (a, f b)-jacobianWith' g f as = second (zipWithT g as) <$> jacobian' f as-{-# INLINE jacobianWith' #-}--grads :: (Traversable f, Num a) => FU f a -> f a -> Cofree f a-grads f as = ds as $ apply f as-{-# INLINE grads #-}--jacobians :: (Traversable f, Functor g, Num a) => FF f g a -> f a -> g (Cofree f a)-jacobians f as = ds as <$> apply f as-{-# INLINE jacobians #-}--d2 :: Functor f => Cofree f a -> f (f a)-d2 = headT . tailT . tailT . tensors -{-# INLINE d2 #-}--d2' :: Functor f => Cofree f a -> (a, f (a, f a))-d2' (a :< as) = (a, fmap (\(da :< das) -> (da, extract <$> das)) as)-{-# INLINE d2' #-}--hessian :: (Traversable f, Num a) => FU f a -> f a -> f (f a)-hessian f as = d2 $ grads f as-{-# INLINE hessian #-}--hessian' :: (Traversable f, Num a) => FU f a -> f a -> (a, f (a, f a))-hessian' f as = d2' $ grads f as-{-# INLINE hessian' #-}--hessianF :: (Traversable f, Functor g, Num a) => FF f g a -> f a -> g (f (f a))-hessianF f as = d2 <$> jacobians f as-{-# INLINE hessianF #-}--hessianF' :: (Traversable f, Functor g, Num a) => FF f g a -> f a -> g (a, f (a, f a))-hessianF' f as = d2' <$> jacobians f as-{-# INLINE hessianF' #-}-

@@ -1,128 +0,0 @@-{-# LANGUAGE Rank2Types, BangPatterns #-}--------------------------------------------------------------------------------- |--- Module      : Numeric.AD.Mode.Tower--- Copyright   : (c) Edward Kmett 2010--- License     : BSD3--- Maintainer  : ekmett@gmail.com--- Stability   : experimental--- Portability : GHC only------ Higher order derivatives via a \"dual number tower\".-----------------------------------------------------------------------------------module Numeric.AD.Mode.Tower-    (-    -- * Taylor Series-      taylor-    , taylor0-    -- * Maclaurin Series-    , maclaurin-    , maclaurin0-    -- * Derivatives-    , diff    -- first derivative of (a -> a) -    , diff'   -- answer and first derivative of (a -> a) -    , diffs   -- answer and all derivatives of (a -> a) -    , diffs0  -- zero padded derivatives of (a -> a)-    , diffsF  -- answer and all derivatives of (a -> f a)-    , diffs0F -- zero padded derivatives of (a -> f a)-    -- * Directional Derivatives-    , du      -- directional derivative of (a -> a)-    , du'     -- answer and directional derivative of (a -> a)-    , dus     -- answer and all directional derivatives of (a -> a) -    , dus0    -- answer and all zero padded directional derivatives of (a -> a)-    , duF     -- directional derivative of (a -> f a)-    , duF'    -- answer and directional derivative of (a -> f a)-    , dusF    -- answer and all directional derivatives of (a -> f a)-    , dus0F   -- answer and all zero padded directional derivatives of (a -> a)-    -- * Exposed Types-    , UU, UF, FU, FF-    , Mode(..)-    , AD(..)-    ) where--import Control.Applicative ((<$>))-import Numeric.AD.Types-import Numeric.AD.Classes-import Numeric.AD.Internal.Tower--diffs :: Num a => UU a -> a -> [a]-diffs f a = getADTower $ apply f a-{-# INLINE diffs #-}--diffs0 :: Num a => UU a -> a -> [a]-diffs0 f a = zeroPad (diffs f a)-{-# INLINE diffs0 #-}--diffsF :: (Functor f, Num a) => UF f a -> a -> f [a]-diffsF f a = getADTower <$> apply f a-{-# INLINE diffsF #-}--diffs0F :: (Functor f, Num a) => UF f a -> a -> f [a]-diffs0F f a = (zeroPad . getADTower) <$> apply f a-{-# INLINE diffs0F #-}--taylor :: Fractional a => UU a -> a -> a -> [a]-taylor f x dx = go 1 1 (diffs f x)-    where-        go !n !acc (a:as) = a * acc : go (n + 1) (acc * dx / n) as-        go _ _ [] = []--taylor0 :: Fractional a => UU a -> a -> a -> [a]-taylor0 f x dx = zeroPad (taylor f x dx)-{-# INLINE taylor0 #-}--maclaurin :: Fractional a => UU a -> a -> [a]-maclaurin f = taylor f 0-{-# INLINE maclaurin #-}--maclaurin0 :: Fractional a => UU a -> a -> [a]-maclaurin0 f = taylor0 f 0-{-# INLINE maclaurin0 #-}--diff :: Num a => UU a -> a -> a-diff f = d . diffs f-{-# INLINE diff #-}--diff' :: Num a => UU a -> a -> (a, a)-diff' f = d' . diffs f-{-# INLINE diff' #-}--du :: (Functor f, Num a) => FU f a -> f (a, a) -> a-du f = d . getADTower . f . fmap withD-{-# INLINE du #-}--du' :: (Functor f, Num a) => FU f a -> f (a, a) -> (a, a)-du' f = d' . getADTower . f . fmap withD-{-# INLINE du' #-}--duF :: (Functor f, Functor g, Num a) => FF f g a -> f (a, a) -> g a-duF f = fmap (d . getADTower) . f . fmap withD-{-# INLINE duF #-}--duF' :: (Functor f, Functor g, Num a) => FF f g a -> f (a, a) -> g (a, a)-duF' f = fmap (d' . getADTower) . f . fmap withD-{-# INLINE duF' #-}--dus :: (Functor f, Num a) => FU f a -> f [a] -> [a]-dus f = getADTower . f . fmap tower-{-# INLINE dus #-}--dus0 :: (Functor f, Num a) => FU f a -> f [a] -> [a]-dus0 f = zeroPad . getADTower . f . fmap tower-{-# INLINE dus0 #-}--dusF :: (Functor f, Functor g, Num a) => FF f g a -> f [a] -> g [a]-dusF f = fmap getADTower . f . fmap tower-{-# INLINE dusF #-}--dus0F :: (Functor f, Functor g, Num a) => FF f g a -> f [a] -> g [a]-dus0F f = fmap getADTower . f . fmap tower-{-# INLINE dus0F #-}---- TODO: higher order gradients--- data f :> a = a :< f (f :> a) --- gradients  :: (Traversable f, Num a) => FU f a -> f a -> f :> a--- gradientsF, jacobians :: (Traversable f, Functor g, Num a) => FF f g a -> f a -> g (f :> a)--- gradientsM :: (Traversable f, Monad m, Num a) => FF f m a -> f a -> m (f :> a)

@@ -1,113 +0,0 @@-{-# LANGUAGE Rank2Types, BangPatterns, ScopedTypeVariables #-}--------------------------------------------------------------------------------- |--- Module      :  Numeric.AD.Newton--- Copyright   :  (c) Edward Kmett 2010--- License     :  BSD3--- Maintainer  :  ekmett@gmail.com--- Stability   :  experimental--- Portability :  GHC only-----------------------------------------------------------------------------------module Numeric.AD.Newton-    (-    -- * Newton's Method (Forward AD)-      findZero-    , inverse-    , fixedPoint-    , extremum-    -- * Gradient Ascent/Descent (Reverse AD)-    , gradientDescent-    , gradientAscent-    -- * Exposed Types-    , UU, UF, FU, FF-    , AD(..)-    , Mode(..)-    ) where--import Prelude hiding (all)-import Data.Foldable (all)-import Data.Traversable (Traversable)-import Numeric.AD.Types-import Numeric.AD.Classes-import Numeric.AD.Mode.Forward (diff, diff')-import Numeric.AD.Mode.Reverse (gradWith')-import Numeric.AD.Internal.Composition---- | The 'findZero' function finds a zero of a scalar function using--- Newton's method; its output is a stream of increasingly accurate--- results.  (Modulo the usual caveats.)------ Examples:------  > take 10 $ findZero (\\x->x^2-4) 1  -- converge to 2.0------  > module Data.Complex---  > take 10 $ findZero ((+1).(^2)) (1 :+ 1)  -- converge to (0 :+ 1)@----findZero :: (Fractional a, Eq a) => UU a -> a -> [a]-findZero f = go-    where-        go x = x : if y == 0 then [] else go (x - y/y') -            where-                (y,y') = diff' f x-{-# INLINE findZero #-}---- | The 'inverseNewton' function inverts a scalar function using--- Newton's method; its output is a stream of increasingly accurate--- results.  (Modulo the usual caveats.)------ Example:------ > take 10 $ inverseNewton sqrt 1 (sqrt 10)  -- converges to 10----inverse :: (Fractional a, Eq a) => UU a -> a -> a -> [a]-inverse f x0 y = findZero (\x -> f x - lift y) x0-{-# INLINE inverse  #-}---- | The 'fixedPoint' function find a fixedpoint of a scalar--- function using Newton's method; its output is a stream of--- increasingly accurate results.  (Modulo the usual caveats.)--- --- > take 10 $ fixedPoint cos 1 -- converges to 0.7390851332151607-fixedPoint :: (Fractional a, Eq a) => UU a -> a -> [a]-fixedPoint f = findZero (\x -> f x - x)-{-# INLINE fixedPoint #-}---- | The 'extremum' function finds an extremum of a scalar--- function using Newton's method; produces a stream of increasingly--- accurate results.  (Modulo the usual caveats.)------ > take 10 $ extremum cos 1 -- convert to 0 -extremum :: (Fractional a, Eq a) => UU a -> a -> [a]-extremum f = findZero (diff (decomposeMode . f . composeMode))-{-# INLINE extremum #-}---- | The 'gradientDescent' function performs a multivariate--- optimization, based on the naive-gradient-descent in the file--- @stalingrad\/examples\/flow-tests\/pre-saddle-1a.vlad@ from the--- VLAD compiler Stalingrad sources.  Its output is a stream of--- increasingly accurate results.  (Modulo the usual caveats.)------ It uses reverse mode automatic differentiation to compute the gradient.-gradientDescent :: (Traversable f, Fractional a, Ord a) => FU f a -> f a -> [f a]-gradientDescent f x0 = go x0 fx0 xgx0 0.1 (0 :: Int)-    where-        (fx0, xgx0) = gradWith' (,) f x0-        go x fx xgx !eta !i-            | eta == 0     = [] -- step size is 0-            | fx1 > fx     = go x fx xgx (eta/2) 0 -- we stepped too far-            | zeroGrad xgx = [] -- gradient is 0-            | otherwise    = x1 : if i == 10-                                  then go x1 fx1 xgx1 (eta*2) 0-                                  else go x1 fx1 xgx1 eta (i+1)-            where-                zeroGrad = all (\(_,g) -> g == 0)-                x1 = fmap (\(xi,gxi) -> xi - eta * gxi) xgx-                (fx1, xgx1) = gradWith' (,) f x1-{-# INLINE gradientDescent #-}--gradientAscent :: (Traversable f, Fractional a, Ord a) => FU f a -> f a -> [f a]-gradientAscent f = gradientDescent (negate . f)-{-# INLINE gradientAscent #-}

@@ -1,51 +0,0 @@-{-# LANGUAGE Rank2Types #-}--------------------------------------------------------------------------------- |--- Module      :  Numeric.AD.Types--- Copyright   :  (c) Edward Kmett 2010--- License     :  BSD3--- Maintainer  :  ekmett@gmail.com--- Stability   :  experimental--- Portability :  GHC only-----------------------------------------------------------------------------------module Numeric.AD.Types-    ( -      AD(..)-    -- * Differentiable Functions-    , UU, UF, FU, FF-    -- * Tensors-    , Tensors(..)-    , headT-    , tailT-    , tensors-    -- * An Identity Mode. -    , Id(..)-    , probe, unprobe-    , probed, unprobed-    -- * Apply functions that use 'lift'-    , lowerUU, lowerUF, lowerFU, lowerFF-    ) where--import Numeric.AD.Internal.Identity-import Numeric.AD.Internal.Types-import Numeric.AD.Internal.Tensors---- these exploit the 'magic' that is probed to avoid the need for Functor, etc.--lowerUU :: UU a -> a -> a-lowerUU f = unprobe . f . probe-{-# INLINE lowerUU #-}--lowerUF :: UF f a -> a -> f a-lowerUF f = unprobed . f . probe-{-# INLINE lowerUF #-}--lowerFU :: FU f a -> f a -> a-lowerFU f = unprobe . f . probed-{-# INLINE lowerFU #-}--lowerFF :: FF f g a -> f a -> g a-lowerFF f = unprobed . f . probed-{-# INLINE lowerFF #-}

@@ -1,8 +1,8 @@ name:         ad-version:      1.4+version:      1.5 license:      BSD3 license-File: LICENSE-copyright:    (c) Edward Kmett 2010-2011,+copyright:    (c) Edward Kmett 2010-2012,               (c) Barak Pearlmutter and Jeffrey Mark Siskind 2008-2009 author:       Edward Kmett maintainer:   ekmett@gmail.com@@ -30,8 +30,7 @@     .     * @Numeric.AD.Mode.Tower@ computes a dense forward-mode AD tower useful for higher derivatives of single input functions.     .-    * @Numeric.AD.Mode.Mixed@ computes using whichever mode or combination thereof is suitable to each individual combinator. This mode is the default, re-exported by @Numeric.AD@-    .+    * @Numeric.AD@ computes using whichever mode or combination thereof is suitable to each individual combinator.     .     While not every mode can provide all operations, the following basic operations are supported, modified as     appropriate by the suffixes below:@@ -60,54 +59,23 @@     .     * @0@ means that the resulting derivative list is padded with 0s at the end.     .-    Changes since 1.3-    .-    * Dependency bump to be compatible with ghc 7.4.1 and mtl 2.1-    .-    * Work on diff (**2) 0-    .-    Changes since 1.2-    .-    * Compiles with template haskell 2.6, changed interface to comply with the lack of Eq and Show as superclasses of Num in the new GHC.-    .-    Changes since 1.1.3-    .-    * Made primal calculations strict where possible.-    .-    Changes since 1.1.0-    .-    * Introduced a much faster topological sort into the reverse mode AD implementation by Anthony Cowley. This fixes a space leak and a stack overflow problem on very large (>2000 variable) problem sets.-    .-    * Made bound calculations in reverse mode more strict.-    .-    Changes since 1.0.0-    .-    * Changed the way 'Show' was derived to comply with changes in instance resolution in ghc >= 7.0 && <= 7.1-    .-    Changes since 0.45.0-    .-    * Converted 'Stream' to use the external 'comonad' package-    .-    Changes since 0.44.5-    .-    * Added Halley's method+    /Changes since 1.3/:     .-    Changes since 0.40.0+    * Moved the contents of @Numeric.AD.Mode.Mixed@ into @Numeric.AD@     .-    * Fixed bug fix for @'(/)' :: (Mode s, Fractional a) => AD s a@+    * Split off @Numeric.AD.Variadic@ for the variadic combinators     .-    * Improved documentation+    * Removed the @UU@, @FU@, @UF@, and @FF@ type aliases.     .-    * Regularized naming conventions+    * Stopped exporting the types for @Mode@ and @AD@ from almost every module. Import @Numeric.AD.Types@ if necessary.     .-    * Exposed 'Id', probe, and lower methods via @Numeric.AD.Types@+    * Renamed @Tensors@ to @Jet@     .-    * Removed monadic combinators+    * Dependency bump to be compatible with ghc 7.4.1 and mtl 2.1     .-    * Retuned the 'Mixed' mode jacobian calculations to only require a 'Functor'-based result.+    * More aggressive zero tracking.     .-    * Added unsafe variadic 'vgrad', 'vgrad'', and 'vgrads' combinators-+    * @diff (**n) 0@ for constant n and @diff (0**)@ both now yield the correct answer for all modes.  source-repository head   type: git@@ -115,6 +83,7 @@  library   extensions: CPP+  hs-source-dirs: src    other-extensions:     BangPatterns@@ -143,14 +112,23 @@    exposed-modules:     Numeric.AD-    Numeric.AD.Classes     Numeric.AD.Types+     Numeric.AD.Newton     Numeric.AD.Halley +    Numeric.AD.Mode.Directed+    Numeric.AD.Mode.Forward+    Numeric.AD.Mode.Reverse+    Numeric.AD.Mode.Tower+    Numeric.AD.Mode.Sparse++    Numeric.AD.Variadic+    Numeric.AD.Variadic.Reverse+    Numeric.AD.Variadic.Sparse+     Numeric.AD.Internal.Classes     Numeric.AD.Internal.Combinators-     Numeric.AD.Internal.Forward     Numeric.AD.Internal.Tower     Numeric.AD.Internal.Reverse@@ -158,16 +136,9 @@     Numeric.AD.Internal.Dense     Numeric.AD.Internal.Composition -    Numeric.AD.Mode.Directed-    Numeric.AD.Mode.Forward-    Numeric.AD.Mode.Mixed-    Numeric.AD.Mode.Reverse-    Numeric.AD.Mode.Tower-    Numeric.AD.Mode.Sparse-   other-modules:     Numeric.AD.Internal.Types-    Numeric.AD.Internal.Tensors+    Numeric.AD.Internal.Jet     Numeric.AD.Internal.Identity    ghc-options: -Wall -fspec-constr -fdicts-cheap -O2

@@ -0,0 +1,215 @@+{-# LANGUAGE Rank2Types, TypeFamilies, PatternGuards #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.AD+-- Copyright   :  (c) Edward Kmett 2010+-- License     :  BSD3+-- Maintainer  :  ekmett@gmail.com+-- Stability   :  experimental+-- Portability :  GHC only+--+-- Mixed-Mode Automatic Differentiation.+--+-- Each combinator exported from this module chooses an appropriate AD mode.+-- The following basic operations are supported, modified as appropriate by the suffixes below:+--+-- * 'grad' computes the gradient (partial derivatives) of a function at a point+--+-- * 'jacobian' computes the Jacobian matrix of a function at a point+--+-- * 'diff' computes the derivative of a function at a point+--+-- * 'du' computes a directional derivative of a function at a point+--+-- * 'hessian' compute the Hessian matrix (matrix of second partial derivatives) of a function at a point+--+-- The suffixes have the following meanings:+--+-- * @\'@ -- also return the answer+--+-- * @With@ lets the user supply a function to blend the input with the output+--+-- * @F@ is a version of the base function lifted to return a 'Traversable' (or 'Functor') result+--+-- * @s@ means the function returns all higher derivatives in a list or f-branching 'Stream'+--+-- * @T@ means the result is transposed with respect to the traditional formulation.+--+-- * @0@ means that the resulting derivative list is padded with 0s at the end.+-----------------------------------------------------------------------------++module Numeric.AD+    (+    -- * Gradients (Reverse Mode)+      grad+    , grad'+    , gradWith+    , gradWith'++    -- * Higher Order Gradients (Sparse-on-Reverse)+    , grads++    -- * Jacobians (Sparse or Reverse)+    , jacobian+    , jacobian'+    , jacobianWith+    , jacobianWith'++    -- * Higher Order Jacobian (Sparse-on-Reverse)+    , jacobians++    -- * Transposed Jacobians (Forward Mode)+    , jacobianT+    , jacobianWithT++    -- * Hessian (Sparse-On-Reverse)+    , hessian+    , hessian'++    -- * Hessian Tensors (Sparse or Sparse-On-Reverse)+    , hessianF+    -- * Hessian Tensors (Sparse)+    , hessianF'++    -- * Hessian Vector Products (Forward-On-Reverse)+    , hessianProduct+    , hessianProduct'++    -- * Derivatives (Forward Mode)+    , diff+    , diffF++    , diff'+    , diffF'++    -- * Derivatives (Tower)+    , diffs+    , diffsF++    , diffs0+    , diffs0F++    -- * Directional Derivatives (Forward Mode)+    , du+    , du'+    , duF+    , duF'++    -- * Directional Derivatives (Tower)+    , dus+    , dus0+    , dusF+    , dus0F++    -- * Taylor Series (Tower)+    , taylor+    , taylor0++    -- * Maclaurin Series (Tower)+    , maclaurin+    , maclaurin0+    ) where++import Data.Traversable (Traversable)+import Data.Foldable (Foldable, foldr')+import Control.Applicative++import Numeric.AD.Types+import Numeric.AD.Internal.Composition+import Numeric.AD.Internal.Identity++import Numeric.AD.Mode.Forward+    ( diff, diff', diffF, diffF'+    , du, du', duF, duF'+    , jacobianT, jacobianWithT )++import Numeric.AD.Mode.Tower+    ( diffsF, diffs0F, diffs, diffs0+    , taylor, taylor0, maclaurin, maclaurin0+    , dus, dus0, dusF, dus0F )++import qualified Numeric.AD.Mode.Reverse as Reverse+import Numeric.AD.Mode.Reverse+    ( grad, grad', gradWith, gradWith')++-- temporary until we make a full sparse mode+import qualified Numeric.AD.Mode.Sparse as Sparse+import Numeric.AD.Mode.Sparse+    ( grads, jacobians, hessian', hessianF')++-- | Calculate the Jacobian of a non-scalar-to-non-scalar function, automatically choosing between forward and reverse mode AD based on the number of inputs and outputs.+--+-- If you know the relative number of inputs and outputs, consider 'Numeric.AD.Reverse.jacobian' or 'Nuneric.AD.Sparse.jacobian'.+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 bs = snd <$> jacobian' f bs+{-# INLINE jacobian #-}++data Nat = Z | S Nat deriving (Eq, Ord)++size :: Foldable f => f a -> Nat+size = foldr' (\_ b -> S b) Z++big :: Nat -> Bool+big (S (S (S (S (S (S (S (S (S (S _)))))))))) = True+big _ = False++-- | Calculate both the answer and Jacobian of a non-scalar-to-non-scalar function, automatically choosing between forward- and reverse- mode AD based on the relative, based on the number of inputs+--+-- If you know the relative number of inputs and outputs, consider 'Numeric.AD.Reverse.jacobian'' or 'Nuneric.AD.Sparse.jacobian''.+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 bs | Z <- n = fmap (\x -> (unprobe x, bs)) (f (probed bs))+               | big n  = Reverse.jacobian' f bs+               | otherwise = Sparse.jacobian' f bs+    where+        n = size bs+{-# INLINE jacobian' #-}++-- | @'jacobianWith' g f@ calculates the Jacobian of a non-scalar-to-non-scalar function, automatically choosing between forward and reverse mode AD based on the number of inputs and outputs.+--+-- The resulting Jacobian matrix is then recombined element-wise with the input using @g@.+--+-- If you know the relative number of inputs and outputs, consider 'Numeric.AD.Reverse.jacobianWith' or 'Nuneric.AD.Sparse.jacobianWith'.+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 bs = snd <$> jacobianWith' g f bs+{-# INLINE jacobianWith #-}++-- | @'jacobianWith'' g f@ calculates the answer and Jacobian of a non-scalar-to-non-scalar function, automatically choosing between sparse and reverse mode AD based on the number of inputs and outputs.+--+-- The resulting Jacobian matrix is then recombined element-wise with the input using @g@.+--+-- If you know the relative number of inputs and outputs, consider 'Numeric.AD.Reverse.jacobianWith'' or 'Nuneric.AD.Sparse.jacobianWith''.+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 bs+    | Z <- n = fmap (\x -> (unprobe x, undefined <$> bs)) (f (probed bs))+    | big n  = Reverse.jacobianWith' g f bs+    | otherwise = Sparse.jacobianWith' g f bs+    where+        n = size bs+{-# INLINE jacobianWith' #-}++-- | @'hessianProduct' f wv@ computes the product of the hessian @H@ of a non-scalar-to-scalar function @f@ at @w = 'fst' <$> wv@ with a vector @v = snd <$> wv@ using \"Pearlmutter\'s method\" from <http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.29.6143>, which states:+--+-- > H v = (d/dr) grad_w (w + r v) | r = 0+--+-- Or in other words, we take the directional derivative of the gradient. The gradient is calculated in reverse mode, then the directional derivative is calculated in forward mode.+--+hessianProduct :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f (a, a) -> f a+hessianProduct f = duF (grad (decomposeMode . f . fmap composeMode))++-- | @'hessianProduct'' f wv@ computes both the gradient of a non-scalar-to-scalar @f@ at @w = 'fst' <$> wv@ and the product of the hessian @H@ at @w@ with a vector @v = snd <$> wv@ using \"Pearlmutter's method\". The outputs are returned wrapped in the same functor.+--+-- > H v = (d/dr) grad_w (w + r v) | r = 0+--+-- Or in other words, we return the gradient and the directional derivative of the gradient. The gradient is calculated in reverse mode, then the directional derivative is calculated in forward mode.+hessianProduct' :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f (a, a) -> f (a, a)+hessianProduct' f = duF' (grad (decomposeMode . f . fmap composeMode))++-- | Compute the Hessian via the Jacobian of the gradient. gradient is computed in reverse mode and then the Jacobian is computed in sparse (forward) mode.+hessian :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f (f a)+hessian f = Sparse.jacobian (grad (decomposeMode . f . fmap composeMode))++-- | Compute the order 3 Hessian tensor on a non-scalar-to-non-scalar function using Sparse or Sparse-on-Reverse+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 as+    | big (size as) = decomposeFunctor $ Sparse.jacobian (ComposeFunctor . Reverse.jacobian (fmap decomposeMode . f . fmap composeMode)) as+    | otherwise = Sparse.hessianF f as

@@ -0,0 +1,80 @@+{-# LANGUAGE Rank2Types, BangPatterns, ScopedTypeVariables #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.AD.Halley+-- Copyright   :  (c) Edward Kmett 2010+-- License     :  BSD3+-- Maintainer  :  ekmett@gmail.com+-- Stability   :  experimental+-- Portability :  GHC only+--+-- Root finding using Halley's rational method (the second in+-- the class of Householder methods). Assumes the function is three+-- times continuously differentiable and converges cubically when+-- progress can be made.+--+-----------------------------------------------------------------------------++module Numeric.AD.Halley+    (+    -- * Halley's Method (Tower AD)+      findZero+    , inverse+    , fixedPoint+    , extremum+    ) where++import Prelude hiding (all)+import Numeric.AD.Types+import Numeric.AD.Mode.Tower (diffs0)+import Numeric.AD.Mode.Forward (diff) -- , diff')+import Numeric.AD.Internal.Composition++-- | The 'findZero' function finds a zero of a scalar function using+-- Halley's method; its output is a stream of increasingly accurate+-- results.  (Modulo the usual caveats.)+--+-- Examples:+--+--  > take 10 $ findZero (\\x->x^2-4) 1  -- converge to 2.0+--+--  > module Data.Complex+--  > take 10 $ findZero ((+1).(^2)) (1 :+ 1)  -- converge to (0 :+ 1)@+--+findZero :: (Fractional a, Eq a) => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]+findZero f = go+    where+        go x = x : if y == 0 then [] else go (x - 2*y*y'/(2*y'*y'-y*y''))+            where+                (y:y':y'':_) = diffs0 f x+{-# INLINE findZero #-}++-- | The 'inverse' function inverts a scalar function using+-- Halley's method; its output is a stream of increasingly accurate+-- results.  (Modulo the usual caveats.)+--+-- 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+{-# INLINE inverse  #-}++-- | The 'fixedPoint' function find a fixedpoint of a scalar+-- function using Halley's method; its output is a stream of+-- increasingly accurate results.  (Modulo the usual caveats.)+--+-- > take 10 $ fixedPoint cos 1 -- converges to 0.7390851332151607+fixedPoint :: (Fractional a, Eq a) => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]+fixedPoint f = findZero (\x -> f x - x)+{-# INLINE fixedPoint #-}++-- | The 'extremum' function finds an extremum of a scalar+-- function using Halley's method; produces a stream of increasingly+-- accurate results.  (Modulo the usual caveats.)+--+-- > take 10 $ extremum cos 1 -- convert to 0+extremum :: (Fractional a, Eq a) => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]+extremum f = findZero (diff (decomposeMode . f . composeMode))+{-# INLINE extremum #-}+

@@ -0,0 +1,328 @@+{-# LANGUAGE Rank2Types, TypeFamilies, FlexibleInstances, MultiParamTypeClasses, PatternGuards, CPP #-}+{-# LANGUAGE FlexibleContexts, FunctionalDependencies, UndecidableInstances, GeneralizedNewtypeDeriving, TemplateHaskell #-}+-- {-# OPTIONS_HADDOCK hide #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.AD.Internal.Classes+-- Copyright   :  (c) Edward Kmett 2010+-- License     :  BSD3+-- Maintainer  :  ekmett@gmail.com+-- Stability   :  experimental+-- Portability :  GHC only+--+-----------------------------------------------------------------------------++module Numeric.AD.Internal.Classes+    (+    -- * AD modes+      Mode(..)+    , one+    -- * Automatically Deriving AD+    , Jacobian(..)+    , Primal(..)+    , deriveLifted+    , deriveNumeric+    , Lifted(..)+    , Iso(..)+    ) where++import Control.Applicative hiding ((<**>))+import Data.Char+import Language.Haskell.TH+import Numeric.AD.Internal.Combinators (on)++infixr 8 **!, <**>+infixl 7 *!, /!, ^*, *^, ^/+infixl 6 +!, -!, <+>+infix 4 ==!++class Iso a b where+    iso :: f a -> f b+    osi :: f b -> f a++instance Iso a a where+    iso = id+    osi = id++class Lifted t where+    showsPrec1          :: (Num a, Show a) => Int -> t a -> ShowS+    (==!)               :: (Num a, Eq a) => t a -> t a -> Bool+    compare1            :: (Num a, Ord a) => t a -> t a -> Ordering+    fromInteger1        :: Num a => Integer -> t a+    (+!),(-!),(*!)      :: Num a => t a -> t a -> t a+    negate1, abs1, signum1 :: Num a => t a -> t a+    (/!)                :: Fractional a => t a -> t a -> t a+    recip1              :: Fractional a => t a -> t a+    fromRational1       :: Fractional a => Rational -> t a+    toRational1         :: Real a => t a -> Rational -- unsafe+    pi1                 :: Floating a => t a+    exp1, log1, sqrt1   :: Floating a => t a -> t a+    (**!), logBase1     :: Floating a => t a -> t a -> t a+    sin1, cos1, tan1, asin1, acos1, atan1 :: Floating a => t a -> t a+    sinh1, cosh1, tanh1, asinh1, acosh1, atanh1 :: Floating a => t a -> t a+    properFraction1 :: (RealFrac a, Integral b) => t a -> (b, t a)+    truncate1, round1, ceiling1, floor1 :: (RealFrac a, Integral b) => t a -> b+    floatRadix1     :: RealFloat a => t a -> Integer+    floatDigits1    :: RealFloat a => t a -> Int+    floatRange1     :: RealFloat a => t a -> (Int, Int)+    decodeFloat1    :: RealFloat a => t a -> (Integer, Int)+    encodeFloat1    :: RealFloat a => Integer -> Int -> t a+    exponent1       :: RealFloat a => t a -> Int+    significand1    :: RealFloat a => t a -> t a+    scaleFloat1     :: RealFloat a => Int -> t a -> t a+    isNaN1, isInfinite1, isDenormalized1, isNegativeZero1, isIEEE1 :: RealFloat a => t a -> Bool+    atan21          :: RealFloat a => t a -> t a -> t a+    succ1, pred1    :: (Num a, Enum a) => t a -> t a+    toEnum1         :: (Num a, Enum a) => Int -> t a+    fromEnum1       :: (Num a, Enum a) => t a -> Int+    enumFrom1       :: (Num a, Enum a) => t a -> [t a]+    enumFromThen1   :: (Num a, Enum a) => t a -> t a -> [t a]+    enumFromTo1     :: (Num a, Enum a) => t a -> t a -> [t a]+    enumFromThenTo1 :: (Num a, Enum a) => t a -> t a -> t a -> [t a]+    minBound1       :: (Num a, Bounded a) => t a+    maxBound1       :: (Num a, Bounded a) => t a++class Lifted t => Mode t where+    -- | allowed to return False for items with a zero derivative, but we'll give more NaNs than strictly necessary+    isKnownConstant :: t a -> Bool+    isKnownConstant _ = False++    -- | allowed to return False for zero, but we give more NaN's than strictly necessary then+    isKnownZero :: Num a => t a -> Bool+    isKnownZero _ = False++    -- | Embed a constant+    lift  :: Num a => a -> t a++    -- | Vector sum+    (<+>) :: Num a => t a -> t a -> t a++    -- | Scalar-vector multiplication+    (*^) :: Num a => a -> t a -> t a++    -- | Vector-scalar multiplication+    (^*) :: Num a => t a -> a -> t a++    -- | Scalar division+    (^/) :: Fractional a => t a -> a -> t a++    -- | Exponentiation, this should be overloaded if you can figure out anything about what is constant!+    (<**>) :: Floating a => t a -> t a -> t a+--  x <**> y = lift2_ (**) (\z xi yi -> (yi *! z /! xi, z *! log1 xi)) x y++    -- | > 'zero' = 'lift' 0+    zero :: Num a => t a++    a *^ b = lift a *! b+    a ^* b = a *! lift b++    a ^/ b = a ^* recip b++    zero = lift 0++one :: (Mode t, Num a) => t a+one = lift 1+{-# INLINE one #-}++negOne :: (Mode t, Num a) => t a+negOne = lift (-1)+{-# INLINE negOne #-}++-- | 'Primal' is used by 'deriveMode' but is not exposed+-- via the 'Mode' class to prevent its abuse by end users+-- 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+-- information. The end user is protected from accidentally using this function+-- by the universal quantification on the various combinators we expose.++class Primal t where+    primal :: Num a => t a -> a++-- | 'Jacobian' is used by 'deriveMode' but is not exposed+-- via 'Mode' to prevent its abuse by end users+-- via the 'AD' data type.+class (Mode t, Mode (D t)) => Jacobian t where+    type D t :: * -> *++    unary  :: Num a => (a -> a) -> D t a -> t a -> t a+    lift1  :: Num a => (a -> a) -> (D t a -> D t a) -> t a -> t a+    lift1_ :: Num a => (a -> a) -> (D t a -> D t a -> D t a) -> t a -> t a++    binary :: Num a => (a -> a -> a) -> D t a -> D t a -> t a -> t a -> t a+    lift2  :: Num a => (a -> a -> a) -> (D t a -> D t a -> (D t a, D t a)) -> t a -> t a -> t a+    lift2_ :: Num a => (a -> a -> a) -> (D t a -> D t a -> D t a -> (D t a, D t a)) -> t a -> t a -> t a++withPrimal :: (Jacobian t, Num a) => t a -> a -> t a+withPrimal t a = unary (const a) one t+{-# INLINE withPrimal #-}++fromBy :: (Jacobian t, Num a) => t a -> t a -> Int -> a -> t a+fromBy a delta n x = binary (\_ _ -> x) one (fromIntegral1 n) a delta++fromIntegral1 :: (Integral n, Lifted t, Num a) => n -> t a+fromIntegral1 = fromInteger1 . fromIntegral+{-# INLINE fromIntegral1 #-}++square1 :: (Lifted t, Num a) => t a -> t a+square1 x = x *! x+{-# INLINE square1 #-}++discrete1 :: (Primal t, Num a) => (a -> c) -> t a -> c+discrete1 f x = f (primal x)+{-# INLINE discrete1 #-}++discrete2 :: (Primal t, Num a) => (a -> a -> c) -> t a -> t a -> c+discrete2 f x y = f (primal x) (primal y)+{-# INLINE discrete2 #-}++discrete3 :: (Primal t, Num a) => (a -> a -> a -> d) -> t a -> t a -> t a -> d+discrete3 f x y z = f (primal x) (primal y) (primal z)+{-# INLINE discrete3 #-}++-- | @'deriveLifted' t@ provides+--+-- > instance Lifted $t+--+-- given supplied instances for+--+-- > instance Lifted $t => Primal $t where ...+-- > instance Lifted $t => Jacobian $t where ...+--+-- The seemingly redundant @'Lifted' $t@ constraints are caused by Template Haskell staging restrictions.+deriveLifted :: ([Q Pred] -> [Q Pred]) -> Q Type -> Q [Dec]+deriveLifted f _t = do+        [InstanceD cxt0 type0 dec0] <- lifted+        return <$> instanceD (cxt (f (return <$> cxt0))) (return type0) (return <$> dec0)+    where+      lifted = [d|+       instance Lifted $_t where+        (==!)         = (==) `on` primal+        compare1      = compare `on` primal+        maxBound1     = lift maxBound+        minBound1     = lift minBound+        showsPrec1 d  = showsPrec d . primal+        fromInteger1 0 = zero+        fromInteger1 n = lift (fromInteger n)+        (+!)          = (<+>) -- binary (+) one one+        (-!)          = binary (-) one negOne -- TODO: <-> ? as it is, this might be pretty bad for Tower+        (*!)          = lift2 (*) (\x y -> (y, x))+        negate1       = lift1 negate (const negOne)+        abs1          = lift1 abs signum1+        signum1       = lift1 signum (const zero)+        fromRational1 0 = zero+        fromRational1 r = lift (fromRational r)+        x /! y        = x *! recip1 y+        recip1        = lift1_ recip (const . negate1 . square1)+        pi1       = lift pi+        exp1      = lift1_ exp const+        log1      = lift1 log recip1+        logBase1 x y = log1 y /! log1 x+        sqrt1     = lift1_ sqrt (\z _ -> recip1 (lift 2 *! z))+        (**!)     = (<**>)+        --x **! y+        --   | isKnownZero y     = 1+        --   | isKnownConstant y, y' <- primal y = lift1 (** y') ((y'*) . (**(y'-1))) x+        --   | otherwise         = lift2_ (**) (\z xi yi -> (yi *! z /! xi, z *! log1 xi)) x y+        sin1      = lift1 sin cos1+        cos1      = lift1 cos $ negate1 . sin1+        tan1 x    = sin1 x /! cos1 x+        asin1     = lift1 asin $ \x -> recip1 (sqrt1 (one -! square1 x))+        acos1     = lift1 acos $ \x -> negate1 (recip1 (sqrt1 (one -! square1 x)))+        atan1     = lift1 atan $ \x -> recip1 (one +! square1 x)+        sinh1     = lift1 sinh cosh1+        cosh1     = lift1 cosh sinh1+        tanh1 x   = sinh1 x /! cosh1 x+        asinh1    = lift1 asinh $ \x -> recip1 (sqrt1 (one +! square1 x))+        acosh1    = lift1 acosh $ \x -> recip1 (sqrt1 (square1 x -! one))+        atanh1    = lift1 atanh $ \x -> recip1 (one -! square1 x)++        succ1                 = lift1 succ (const one)+        pred1                 = lift1 pred (const one)+        toEnum1               = lift . toEnum+        fromEnum1             = discrete1 fromEnum+        enumFrom1 a           = withPrimal a <$> discrete1 enumFrom a+        enumFromTo1 a b       = withPrimal a <$> discrete2 enumFromTo a b+        enumFromThen1 a b     = zipWith (fromBy a delta) [0..] $ discrete2 enumFromThen a b where delta = b -! a+        enumFromThenTo1 a b c = zipWith (fromBy a delta) [0..] $ discrete3 enumFromThenTo a b c where delta = b -! a++        toRational1      = discrete1 toRational+        floatRadix1      = discrete1 floatRadix+        floatDigits1     = discrete1 floatDigits+        floatRange1      = discrete1 floatRange+        decodeFloat1     = discrete1 decodeFloat+        encodeFloat1 m e = lift (encodeFloat m e)+        isNaN1           = discrete1 isNaN+        isInfinite1      = discrete1 isInfinite+        isDenormalized1  = discrete1 isDenormalized+        isNegativeZero1  = discrete1 isNegativeZero+        isIEEE1          = discrete1 isIEEE+        exponent1 = exponent . primal+        scaleFloat1 n = unary (scaleFloat n) (scaleFloat1 n one)+        significand1 x =  unary significand (scaleFloat1 (- floatDigits1 x) one) x+        atan21 = lift2 atan2 $ \vx vy -> let r = recip1 (square1 vx +! square1 vy) in (vy *! r, negate1 vx *! r)+        properFraction1 a = (w, a `withPrimal` pb) where+             pa = primal a+             (w, pb) = properFraction pa+        truncate1 = discrete1 truncate+        round1    = discrete1 round+        ceiling1  = discrete1 ceiling+        floor1    = discrete1 floor |]++varA :: Q Type+varA = varT (mkName "a")++-- | Find all the members defined in the 'Lifted' data type+liftedMembers :: Q [String]+liftedMembers = do+#ifdef OldClassI+    ClassI (ClassD _ _ _ _ ds) <- reify ''Lifted+#else+    ClassI (ClassD _ _ _ _ ds) _ <- reify ''Lifted+#endif+    return [ nameBase n | SigD n _ <- ds]++-- | @'deriveNumeric' f g@ provides the following instances:+--+-- > instance ('Lifted' $f, 'Num' a, 'Enum' a) => 'Enum' ($g a)+-- > instance ('Lifted' $f, 'Num' a, 'Eq' a) => 'Eq' ($g a)+-- > instance ('Lifted' $f, 'Num' a, 'Ord' a) => 'Ord' ($g a)+-- > instance ('Lifted' $f, 'Num' a, 'Bounded' a) => 'Bounded' ($g a)+--+-- > instance ('Lifted' $f, 'Show' a) => 'Show' ($g a)+-- > instance ('Lifted' $f, 'Num' a) => 'Num' ($g a)+-- > instance ('Lifted' $f, 'Fractional' a) => 'Fractional' ($g a)+-- > instance ('Lifted' $f, 'Floating' a) => 'Floating' ($g a)+-- > instance ('Lifted' $f, 'RealFloat' a) => 'RealFloat' ($g a)+-- > instance ('Lifted' $f, 'RealFrac' a) => 'RealFrac' ($g a)+-- > instance ('Lifted' $f, 'Real' a) => 'Real' ($g a)+deriveNumeric :: ([Q Pred] -> [Q Pred]) -> Q Type -> Q [Dec]+deriveNumeric f t = do+    members <- liftedMembers+    let keep n = nameBase n `elem` members+    xs <- lowerInstance keep ((classP ''Num [varA]:) . f) t `mapM` [''Enum, ''Eq, ''Ord, ''Bounded, ''Show]+    ys <- lowerInstance keep f                            t `mapM` [''Num, ''Fractional, ''Floating, ''RealFloat,''RealFrac, ''Real]+    return (xs ++ ys)++lowerInstance :: (Name -> Bool) -> ([Q Pred] -> [Q Pred]) -> Q Type -> Name -> Q Dec+lowerInstance p f t n = do+#ifdef OldClassI+    ClassI (ClassD _ _ _ _ ds) <- reify n+#else+    ClassI (ClassD _ _ _ _ ds) _ <- reify n+#endif+    instanceD (cxt (f [classP n [varA]]))+              (conT n `appT` (t `appT` varA))+              (concatMap lower1 ds)+    where+        lower1 :: Dec -> [Q Dec]+        lower1 (SigD n' _) | p n'' = [valD (varP n') (normalB (varE n'')) []] where n'' = primed n'+        lower1 _          = []++        primed n' = mkName $ base ++ [prime]+            where+                base = nameBase n'+                h = head base+                prime | isSymbol h || h `elem` "/*-<>" = '!'+                      | otherwise = '1'

@@ -0,0 +1,28 @@+{-# OPTIONS_GHC -fno-warn-incomplete-patterns #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.AD.Internal.Combinators+-- Copyright   :  (c) Edward Kmett 2010+-- License     :  BSD3+-- Maintainer  :  ekmett@gmail.com+-- Stability   :  experimental+-- Portability :  GHC only+--+-----------------------------------------------------------------------------+module Numeric.AD.Internal.Combinators+    ( zipWithT+    , zipWithDefaultT+    , on+    ) where++import Data.Traversable (Traversable, mapAccumL)+import Data.Foldable (Foldable, toList)++on :: (a -> a -> b) -> (c -> a) -> c -> c -> b+on f g a b = f (g a) (g b)++zipWithT :: (Foldable f, Traversable g) => (a -> b -> c) -> f a -> g b -> g c+zipWithT f as = snd . mapAccumL (\(a:as') b -> (as', f a b)) (toList as)++zipWithDefaultT :: (Foldable f, Traversable g) => a -> (a -> b -> c) -> f a -> g b -> g c+zipWithDefaultT z f as = zipWithT f (toList as ++ repeat z)

@@ -0,0 +1,183 @@+{-# LANGUAGE Rank2Types, TypeFamilies, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, FlexibleContexts, TemplateHaskell, UndecidableInstances, TypeOperators #-}+-- {-# OPTIONS_HADDOCK hide, prune #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.AD.Internal.Composition+-- Copyright   :  (c) Edward Kmett 2010+-- License     :  BSD3+-- Maintainer  :  ekmett@gmail.com+-- Stability   :  experimental+-- Portability :  GHC only+--+-----------------------------------------------------------------------------++module Numeric.AD.Internal.Composition+    ( ComposeFunctor(..)+    , ComposeMode(..)+    , composeMode+    , decomposeMode+    ) where++import Control.Applicative hiding ((<**>))+import Data.Data (Data(..), mkDataType, DataType, mkConstr, Constr, constrIndex, Fixity(..))+import Data.Typeable (Typeable1(..), Typeable(..), TyCon, mkTyCon3, mkTyConApp, typeOfDefault, gcast1)+import Data.Foldable (Foldable(foldMap))+import Data.Traversable (Traversable(traverse))+import Numeric.AD.Internal.Classes+import Numeric.AD.Internal.Types++-- | Functor composition, used to nest the use of jacobian and grad+newtype ComposeFunctor f g a = ComposeFunctor { decomposeFunctor :: f (g a) }++instance (Functor f, Functor g) => Functor (ComposeFunctor f g) where+    fmap f (ComposeFunctor a) = ComposeFunctor (fmap (fmap f) a)++instance (Foldable f, Foldable g) => Foldable (ComposeFunctor f g) where+    foldMap f (ComposeFunctor a) = foldMap (foldMap f) a++instance (Traversable f, Traversable g) => Traversable (ComposeFunctor f g) where+    traverse f (ComposeFunctor a) = ComposeFunctor <$> traverse (traverse f) a++instance (Typeable1 f, Typeable1 g) => Typeable1 (ComposeFunctor f g) where+    typeOf1 tfga = mkTyConApp composeFunctorTyCon [typeOf1 (fa tfga), typeOf1 (ga tfga)]+        where fa :: t f (g :: * -> *) a -> f a+              fa = undefined+              ga :: t (f :: * -> *) g a -> g a+              ga = undefined++composeFunctorTyCon :: TyCon+composeFunctorTyCon = mkTyCon3 "ad" "Numeric.AD.Internal.Composition" "ComposeFunctor"+{-# NOINLINE composeFunctorTyCon #-}++composeFunctorConstr :: Constr+composeFunctorConstr = mkConstr composeFunctorDataType "ComposeFunctor" [] Prefix+{-# NOINLINE composeFunctorConstr #-}++composeFunctorDataType :: DataType+composeFunctorDataType = mkDataType "Numeric.AD.Internal.Composition.ComposeFunctor" [composeFunctorConstr]+{-# NOINLINE composeFunctorDataType #-}++instance (Typeable1 f, Typeable1 g, Data (f (g a)), Data a) => Data (ComposeFunctor f g a) where+    gfoldl f z (ComposeFunctor a) = z ComposeFunctor `f` a+    toConstr _ = composeFunctorConstr+    gunfold k z c = case constrIndex c of+        1 -> k (z ComposeFunctor)+        _ -> error "gunfold"+    dataTypeOf _ = composeFunctorDataType+    dataCast1 f = gcast1 f++-- | The composition of two AD modes is an AD mode in its own right+newtype ComposeMode f g a = ComposeMode { runComposeMode :: f (AD g a) }++composeMode :: AD f (AD g a) -> AD (ComposeMode f g) a+composeMode (AD a) = AD (ComposeMode a)++decomposeMode :: AD (ComposeMode f g) a -> AD f (AD g a)+decomposeMode (AD (ComposeMode a)) = AD a++instance (Primal f, Mode g, Primal g) => Primal (ComposeMode f g) where+    primal = primal . primal . runComposeMode++instance (Mode f, Mode g) => Mode (ComposeMode f g) where+    lift = ComposeMode . lift . lift+    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)+    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+    ComposeMode a +! ComposeMode b = ComposeMode (a +! b)+    ComposeMode a -! ComposeMode b = ComposeMode (a -! b)+    ComposeMode a *! ComposeMode b = ComposeMode (a *! b)+    negate1 (ComposeMode a) = ComposeMode (negate1 a)+    abs1 (ComposeMode a) = ComposeMode (abs1 a)+    signum1 (ComposeMode a) = ComposeMode (signum1 a)+    ComposeMode a /! ComposeMode b = ComposeMode (a /! b)+    recip1 (ComposeMode a) = ComposeMode (recip1 a)+    fromRational1 = ComposeMode . lift . fromRational1+    toRational1 (ComposeMode a) = toRational1 a+    pi1 = ComposeMode pi1+    exp1 (ComposeMode a) = ComposeMode (exp1 a)+    log1 (ComposeMode a) = ComposeMode (log1 a)+    sqrt1 (ComposeMode a) = ComposeMode (sqrt1 a)+    ComposeMode a **! ComposeMode b = ComposeMode (a **! b)+    logBase1 (ComposeMode a) (ComposeMode b) = ComposeMode (logBase1 a b)+    sin1 (ComposeMode a) = ComposeMode (sin1 a)+    cos1 (ComposeMode a) = ComposeMode (cos1 a)+    tan1 (ComposeMode a) = ComposeMode (tan1 a)+    asin1 (ComposeMode a) = ComposeMode (asin1 a)+    acos1 (ComposeMode a) = ComposeMode (acos1 a)+    atan1 (ComposeMode a) = ComposeMode (atan1 a)+    sinh1 (ComposeMode a) = ComposeMode (sinh1 a)+    cosh1 (ComposeMode a) = ComposeMode (cosh1 a)+    tanh1 (ComposeMode a) = ComposeMode (tanh1 a)+    asinh1 (ComposeMode a) = ComposeMode (asinh1 a)+    acosh1 (ComposeMode a) = ComposeMode (acosh1 a)+    atanh1 (ComposeMode a) = ComposeMode (atanh1 a)+    properFraction1 (ComposeMode a) = (b, ComposeMode c) where+        (b, c) = properFraction1 a+    truncate1 (ComposeMode a) = truncate1 a+    round1 (ComposeMode a) = round1 a+    ceiling1 (ComposeMode a) = ceiling1 a+    floor1 (ComposeMode a) = floor1 a+    floatRadix1 (ComposeMode a) = floatRadix1 a+    floatDigits1 (ComposeMode a) = floatDigits1 a+    floatRange1 (ComposeMode a) = floatRange1 a+    decodeFloat1 (ComposeMode a) = decodeFloat1 a+    encodeFloat1 m e = ComposeMode (encodeFloat1 m e)+    exponent1 (ComposeMode a) = exponent1 a+    significand1 (ComposeMode a) = ComposeMode (significand1 a)+    scaleFloat1 n (ComposeMode a) = ComposeMode (scaleFloat1 n a)+    isNaN1 (ComposeMode a) = isNaN1 a+    isInfinite1 (ComposeMode a) = isInfinite1 a+    isDenormalized1 (ComposeMode a) = isDenormalized1 a+    isNegativeZero1 (ComposeMode a) = isNegativeZero1 a+    isIEEE1 (ComposeMode a) = isIEEE1 a+    atan21 (ComposeMode a) (ComposeMode b) = ComposeMode (atan21 a b)+    succ1 (ComposeMode a) = ComposeMode (succ1 a)+    pred1 (ComposeMode a) = ComposeMode (pred1 a)+    toEnum1 n = ComposeMode (toEnum1 n)+    fromEnum1 (ComposeMode a) = fromEnum1 a+    enumFrom1 (ComposeMode a) = map ComposeMode $ enumFrom1 a+    enumFromThen1 (ComposeMode a) (ComposeMode b) = map ComposeMode $ enumFromThen1 a b+    enumFromTo1 (ComposeMode a) (ComposeMode b) = map ComposeMode $ enumFromTo1 a b+    enumFromThenTo1 (ComposeMode a) (ComposeMode b) (ComposeMode c) = map ComposeMode $ enumFromThenTo1 a b c+    minBound1 = ComposeMode minBound1+    maxBound1 = ComposeMode maxBound1++instance (Typeable1 f, Typeable1 g) => Typeable1 (ComposeMode f g) where+    typeOf1 tfga = mkTyConApp composeModeTyCon [typeOf1 (fa tfga), typeOf1 (ga tfga)]+        where fa :: t f (g :: * -> *) a -> f a+              fa = undefined+              ga :: t (f :: * -> *) g a -> g a+              ga = undefined++instance (Typeable1 f, Typeable1 g, Typeable a) => Typeable (ComposeMode f g a) where+    typeOf = typeOfDefault+    +composeModeTyCon :: TyCon+composeModeTyCon = mkTyCon3 "ad" "Numeric.AD.Internal.Composition" "ComposeMode"+{-# NOINLINE composeModeTyCon #-}++composeModeConstr :: Constr+composeModeConstr = mkConstr composeModeDataType "ComposeMode" [] Prefix+{-# NOINLINE composeModeConstr #-}++composeModeDataType :: DataType+composeModeDataType = mkDataType "Numeric.AD.Internal.Composition.ComposeMode" [composeModeConstr]+{-# NOINLINE composeModeDataType #-}++instance (Typeable1 f, Typeable1 g, Data (f (AD g a)), Data a) => Data (ComposeMode f g a) where+    gfoldl f z (ComposeMode a) = z ComposeMode `f` a+    toConstr _ = composeModeConstr+    gunfold k z c = case constrIndex c of+        1 -> k (z ComposeMode)+        _ -> error "gunfold"+    dataTypeOf _ = composeModeDataType+    dataCast1 f = gcast1 f+

@@ -0,0 +1,186 @@+{-# LANGUAGE Rank2Types, TypeFamilies, FlexibleContexts, UndecidableInstances, TemplateHaskell, DeriveDataTypeable, BangPatterns #-}+-- {-# OPTIONS_HADDOCK hide, prune #-}+-----------------------------------------------------------------------------+-- |+-- Module      : Numeric.AD.Internal.Dense+-- Copyright   : (c) Edward Kmett 2010+-- License     : BSD3+-- Maintainer  : ekmett@gmail.com+-- Stability   : experimental+-- Portability : GHC only+--+-- Dense Forward AD. Useful when the result involves the majority of the input+-- elements. Do not use for 'Numeric.AD.Mode.Mixed.hessian' and beyond, since+-- they only contain a small number of unique @n@th derivatives --+-- @(n + k - 1) `choose` k@ for functions of @k@ inputs rather than the+-- @k^n@ that would be generated by using 'Dense', not to mention the redundant+-- intermediate derivatives that would be+-- calculated over and over during that process!+--+-- Assumes all instances of 'f' have the same number of elements.+--+-- NB: We don't need the full power of 'Traversable' here, we could get+-- by with a notion of zippable that can plug in 0's for the missing+-- entries. This might allow for gradients where @f@ has exponentials like @((->) a)@+-----------------------------------------------------------------------------++module Numeric.AD.Internal.Dense+    ( Dense(..)+    , ds+    , ds'+    , vars+    , apply+    ) where++import Language.Haskell.TH+import Data.Typeable ()+import Data.Traversable (Traversable, mapAccumL)+import Data.Data ()+import Numeric.AD.Internal.Types+import Numeric.AD.Internal.Combinators+import Numeric.AD.Internal.Classes+import Numeric.AD.Internal.Identity++data Dense f a+    = Lift !a+    | Dense !a (f a)+    | Zero++instance Show a => Show (Dense f a) where+    showsPrec d (Lift a)    = showsPrec d a+    showsPrec d (Dense a _) = showsPrec d a+    showsPrec _ Zero        = showString "0"++ds :: f a -> AD (Dense f) a -> f a+ds _ (AD (Dense _ da)) = da+ds z _ = z+{-# INLINE ds #-}++ds' :: Num a => f a -> AD (Dense f) a -> (a, f a)+ds' _ (AD (Dense a da)) = (a, da)+ds' z (AD (Lift a)) = (a, z)+ds' z (AD Zero) = (0, z)+{-# INLINE ds' #-}++-- Bind variables and count inputs+vars :: (Traversable f, Num a) => f a -> f (AD (Dense f) a)+vars as = snd $ mapAccumL outer (0 :: Int) as+    where+        outer !i a = (i + 1, AD $ Dense a $ snd $ mapAccumL (inner i) 0 as)+        inner !i !j _ = (j + 1, if i == j then 1 else 0)+{-# INLINE vars #-}++apply :: (Traversable f, Num a) => (f (AD (Dense f) a) -> b) -> f a -> b+apply f as = f (vars as)+{-# INLINE apply #-}++instance Primal (Dense f) where+    primal Zero = 0+    primal (Lift a) = a+    primal (Dense a _) = a++instance (Traversable f, Lifted (Dense f)) => Mode (Dense f) where+    lift = Lift+    zero = Zero++    Zero <+> a = a+    a <+> Zero = a+    Lift a     <+> Lift b     = Lift (a + b)+    Lift a     <+> Dense b db = Dense (a + b) db+    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+    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++    _ *^ Zero       = Zero+    a *^ Lift b     = Lift (a * b)+    a *^ Dense b db = Dense (a * b) $ fmap (a*) db+    Zero       ^* _ = Zero+    Lift a     ^* b = Lift (a * b)+    Dense a da ^* b = Dense (a * b) $ fmap (*b) da+    Zero       ^/ _ = Zero+    Lift a     ^/ b = Lift (a / b)+    Dense a da ^/ b = Dense (a / b) $ fmap (/b) da++instance (Traversable f, Lifted (Dense f)) => Jacobian (Dense f) where+    type D (Dense f) = Id+    unary f _         Zero        = Lift (f 0)+    unary f _         (Lift b)    = Lift (f b)+    unary f (Id dadb) (Dense b db) = Dense (f b) (fmap (dadb *) db)++    lift1 f _  Zero        = Lift (f 0)+    lift1 f _  (Lift b)    = Lift (f b)+    lift1 f df (Dense b db) = Dense (f b) (fmap (dadb *) db)+        where+            Id dadb = df (Id b)++    lift1_ f _  Zero         = Lift (f 0)+    lift1_ f _  (Lift b)     = Lift (f b)+    lift1_ f df (Dense b db) = Dense a (fmap (dadb *) db)+        where+            a = f b+            Id dadb = df (Id a) (Id b)++    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         (Dense c dc) = Dense (f 0 c) $ fmap (* dadc) dc+    binary f _         (Id dadc) (Lift b)     (Dense c dc) = Dense (f b c) $ fmap (* dadc) dc+    binary f (Id dadb) _         (Dense b db) Zero         = Dense (f b 0) $ fmap (dadb *) db+    binary f (Id dadb) _         (Dense b db) (Lift c)     = Dense (f b c) $ fmap (dadb *) db+    binary f (Id dadb) (Id dadc) (Dense b db) (Dense c dc) = Dense (f b c) $ zipWithT productRule db dc+        where productRule dbi dci = dadb * dbi + dci * dadc++    lift2 f _  Zero         Zero         = Lift (f 0 0)+    lift2 f _  Zero         (Lift c)     = Lift (f 0 c)+    lift2 f _  (Lift b)     Zero         = Lift (f b 0)+    lift2 f _  (Lift b)     (Lift c)     = Lift (f b c)+    lift2 f df Zero         (Dense c dc) = Dense (f 0 c) $ fmap (*dadc) dc where dadc = runId (snd (df (Id 0) (Id c)))+    lift2 f df (Lift b)     (Dense c dc) = Dense (f b c) $ fmap (*dadc) dc where dadc = runId (snd (df (Id b) (Id c)))+    lift2 f df (Dense b db) Zero         = Dense (f b 0) $ fmap (dadb*) db where dadb = runId (fst (df (Id b) (Id 0)))+    lift2 f df (Dense b db) (Lift c)     = Dense (f b c) $ fmap (dadb*) db where dadb = runId (fst (df (Id b) (Id c)))+    lift2 f df (Dense b db) (Dense c dc) = Dense (f b c) da+        where+            (Id dadb, Id dadc) = df (Id b) (Id c)+            da = zipWithT productRule db dc+            productRule dbi dci = dadb * dbi + dci * dadc++    lift2_ f _  Zero     Zero     = Lift (f 0 0)+    lift2_ f _  Zero     (Lift c) = Lift (f 0 c)+    lift2_ f _  (Lift b) Zero     = Lift (f b 0)+    lift2_ f _  (Lift b) (Lift c) = Lift (f b c)+    lift2_ f df Zero     (Dense c dc)+        = Dense a $ fmap (*dadc) dc+        where+            a = f 0 c+            (_, Id dadc) = df (Id a) (Id 0) (Id c)+    lift2_ f df (Lift b) (Dense c dc)+        = Dense a $ fmap (*dadc) dc+        where+            a = f b c+            (_, Id dadc) = df (Id a) (Id b) (Id c)+    lift2_ f df (Dense b db) Zero+        = Dense a $ fmap (dadb*) db+        where+            a = f b 0+            (Id dadb, _) = df (Id a) (Id b) (Id 0)+    lift2_ f df (Dense b db) (Lift c)+        = Dense a $ fmap (dadb*) db+        where+            a = f b c+            (Id dadb, _) = df (Id a) (Id b) (Id c)+    lift2_ f df (Dense b db) (Dense c dc)+        = Dense a $ zipWithT productRule db dc+        where+            a = f b c+            (Id dadb, Id dadc) = df (Id a) (Id b) (Id c)+            productRule dbi dci = dadb * dbi + dci * dadc++let f = varT (mkName "f") in+    deriveLifted+        (classP ''Traversable [f]:)+        (conT ''Dense `appT` f)

@@ -0,0 +1,200 @@+{-# LANGUAGE Rank2Types, TypeFamilies, DeriveDataTypeable, TemplateHaskell, UndecidableInstances, BangPatterns #-}+-- {-# OPTIONS_HADDOCK hide, prune #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.AD.Internal.Forward+-- Copyright   :  (c) Edward Kmett 2010+-- License     :  BSD3+-- Maintainer  :  ekmett@gmail.com+-- Stability   :  experimental+-- Portability :  GHC only+--+-- Unsafe and often partial combinators intended for internal usage.+--+-- Handle with care.+-----------------------------------------------------------------------------++module Numeric.AD.Internal.Forward+    ( Forward(..)+    , tangent+    , bundle+    , unbundle+    , apply+    , bind+    , bind'+    , bindWith+    , bindWith'+    , transposeWith+    ) where++import Language.Haskell.TH+import Data.Typeable+import Data.Traversable (Traversable, mapAccumL)+import Data.Foldable (Foldable, toList)+import Data.Data+import Control.Applicative+import Numeric.AD.Internal.Types+import Numeric.AD.Internal.Classes+import Numeric.AD.Internal.Identity++data Forward a+  = Forward !a a+  | Lift !a+  | Zero+  deriving (Show, Data, Typeable)++tangent :: Num a => AD Forward a -> a+tangent (AD (Forward _ da)) = da+tangent _ = 0+{-# INLINE tangent #-}++unbundle :: Num a => AD Forward a -> (a, a)+unbundle (AD (Forward a da)) = (a, da)+unbundle (AD Zero) = (0,0)+unbundle (AD (Lift a)) = (a, 0)+{-# INLINE unbundle #-}++bundle :: a -> a -> AD Forward a+bundle a da = AD (Forward a da)+{-# INLINE bundle #-}++apply :: Num a => (AD Forward a -> b) -> a -> b+apply f a = f (bundle a 1)+{-# INLINE apply #-}++instance Primal Forward where+    primal (Forward a _) = a+    primal (Lift a) = a+    primal Zero = 0++instance Lifted Forward => Mode Forward where+    lift = Lift+    zero = Zero++    isKnownZero Zero = True+    isKnownZero _    = False++    isKnownConstant Forward{} = False+    isKnownConstant _ = True++    Zero <+> a = a+    a <+> Zero = a+    Forward a da <+> Forward b db = Forward (a + b) (da + db)+    Forward a da <+> Lift b = Forward (a + b) da+    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+    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++    a *^ Forward b db = Forward (a * b) (a * db)+    a *^ Lift b = Lift (a * b)+    _ *^ Zero = Zero++    Forward a da ^* b = Forward (a * b) (da * b)+    Lift a ^* b = Lift (a * b)+    Zero ^* _ = Zero++    Forward a da ^/ b = Forward (a / b) (da / b)+    Lift a ^/ b = Lift (a / b)+    Zero ^/ _ = Zero++instance Lifted Forward => Jacobian Forward where+    type D Forward = Id+++    unary f (Id dadb) (Forward b db) = Forward (f b) (dadb * db)+    unary f _         (Lift b)       = Lift (f b)+    unary f _         Zero           = Lift (f 0)++    lift1 f _ Zero            = Lift (f 0)+    lift1 f _  (Lift b)       = Lift (f b)+    lift1 f df (Forward b db) = Forward (f b) (dadb * db)+        where+            Id dadb = df (Id b)++    lift1_ f _  Zero           = Lift (f 0)+    lift1_ f _  (Lift b)       = Lift (f b)+    lift1_ f df (Forward b db) = Forward a da+        where+            a = f b+            Id da = df (Id a) (Id b) ^* db++    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           (Forward c dc) = Forward (f 0 c) $ dc * dadc+    binary f _         (Id dadc) (Lift b)       (Forward c dc) = Forward (f b c) $ dc * dadc+    binary f (Id dadb) _         (Forward b db) Zero           = Forward (f b 0) $ dadb * db+    binary f (Id dadb) _         (Forward b db) (Lift c)       = Forward (f b c) $ dadb * db+    binary f (Id dadb) (Id dadc) (Forward b db) (Forward c dc) = Forward (f b c) $ dadb * db + dc * dadc++    lift2 f _  Zero           Zero           = Lift (f 0 0)+    lift2 f _  Zero           (Lift c)       = Lift (f 0 c)+    lift2 f _  (Lift b)       Zero           = Lift (f b 0)+    lift2 f _  (Lift b)       (Lift c)       = Lift (f b c)+    lift2 f df Zero           (Forward c dc) = Forward (f 0 c) $ dc * runId (snd (df (Id 0) (Id c)))+    lift2 f df (Lift b)       (Forward c dc) = Forward (f b c) $ dc * runId (snd (df (Id b) (Id c)))+    lift2 f df (Forward b db) Zero           = Forward (f b 0) $ runId (fst (df (Id b) (Id 0))) * db+    lift2 f df (Forward b db) (Lift c)       = Forward (f b c) $ runId (fst (df (Id b) (Id c))) * db+    lift2 f df (Forward b db) (Forward c dc) = Forward a da+        where+            a = f b c+            (Id dadb, Id dadc) = df (Id b) (Id c)+            da = dadb * db + dc * dadc++    lift2_ f _  Zero           Zero           = Lift (f 0 0)+    lift2_ f _  Zero           (Lift c)       = Lift (f 0 c)+    lift2_ f _  (Lift b)       Zero           = Lift (f b 0)+    lift2_ f _  (Lift b)       (Lift c)       = Lift (f b c)+    lift2_ f df Zero           (Forward c dc) = Forward a $ dc * runId (snd (df (Id a) (Id 0) (Id c))) where a = f 0 c+    lift2_ f df (Lift b)       (Forward c dc) = Forward a $ dc * runId (snd (df (Id a) (Id b) (Id c))) where a = f b c+    lift2_ f df (Forward b db) Zero           = Forward a $ runId (fst (df (Id a) (Id b) (Id 0))) * db where a = f b 0+    lift2_ f df (Forward b db) (Lift c)       = Forward a $ runId (fst (df (Id a) (Id b) (Id c))) * db where a = f b c+    lift2_ f df (Forward b db) (Forward c dc) = Forward a da+        where+            a = f b c+            (Id dadb, Id dadc) = df (Id a) (Id b) (Id c)+            da = dadb * db + dc * dadc++deriveLifted id $ conT ''Forward++bind :: (Traversable f, Num a) => (f (AD Forward a) -> b) -> f a -> f b+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 AD Zero)++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 AD Zero)+        b0 = f (lift <$> 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)++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)+        dropIx ((_,b),bs) = (b,bs)++-- we can't transpose arbitrary traversables, since we can't construct one out of whole cloth, and the outer+-- traversable could be empty. So instead we use one as a 'skeleton'+transposeWith :: (Functor f, Foldable f, Traversable g) => (b -> f a -> c) -> f (g a) -> g b -> g c+transposeWith f as = snd . mapAccumL go xss0+    where+        go xss b = (tail <$> xss, f b (head <$> xss))+        xss0 = toList <$> as+

@@ -0,0 +1,139 @@+{-# LANGUAGE GeneralizedNewtypeDeriving, FlexibleContexts, FlexibleInstances, MultiParamTypeClasses, DeriveDataTypeable #-}+{-# OPTIONS_HADDOCK hide #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.AD.Internal.Identity+-- Copyright   :  (c) Edward Kmett 2010+-- License     :  BSD3+-- Maintainer  :  ekmett@gmail.com+-- Stability   :  experimental+-- Portability :  GHC only+--+-----------------------------------------------------------------------------+module Numeric.AD.Internal.Identity+    ( Id(..)+    , probe+    , unprobe+    , probed+    , unprobed+    ) where++import Control.Applicative+import Numeric.AD.Internal.Classes+import Numeric.AD.Internal.Types+import Data.Monoid+import Data.Data (Data)+import Data.Typeable (Typeable)+import Data.Traversable (Traversable, traverse)+import Data.Foldable (Foldable, foldMap)++newtype Id a = Id { runId :: a } deriving+    (Iso a, Eq, Ord, Show, Enum, Bounded, Num, Real, Fractional, Floating, RealFrac, RealFloat, Monoid, Data, Typeable)++probe :: a -> AD Id a+probe a = AD (Id a)++unprobe :: AD Id a -> a+unprobe (AD (Id a)) = a++pid :: f a -> f (Id a)+pid = iso++unpid :: f (Id a) -> f a+unpid = osi++probed :: f a -> f (AD Id a)+probed = iso . pid++unprobed :: f (AD Id a) -> f a+unprobed = unpid . osi++instance Functor Id where+    fmap f (Id a) = Id (f a)++instance Foldable Id where+    foldMap f (Id a) = f a++instance Traversable Id where+    traverse f (Id a) = Id <$> f a++instance Applicative Id where+    pure = Id+    Id f <*> Id a = Id (f a)++instance Monad Id where+    return = Id+    Id a >>= f = f a++instance Lifted Id where+    (==!) = (==)+    compare1 = compare+    showsPrec1 = showsPrec+    fromInteger1 = fromInteger+    (+!) = (+)+    (-!) = (-)+    (*!) = (*)+    negate1 = negate+    abs1 = abs+    signum1 = signum+    (/!) = (/)+    recip1 = recip+    fromRational1 = fromRational+    toRational1 = toRational+    pi1 = pi+    exp1 = exp+    log1 = log+    sqrt1 = sqrt+    (**!) = (**)+    logBase1 = logBase+    sin1 = sin+    cos1 = cos+    tan1 = tan+    asin1 = asin+    acos1 = acos+    atan1 = atan+    sinh1 = sinh+    cosh1 = cosh+    tanh1 = tanh+    asinh1 = asinh+    acosh1 = acosh+    atanh1 = atanh+    properFraction1 = properFraction+    truncate1 = truncate+    round1 = round+    ceiling1 = ceiling+    floor1 = floor+    floatRadix1 = floatRadix+    floatDigits1 = floatDigits+    floatRange1 = floatRange+    decodeFloat1 = decodeFloat+    encodeFloat1 = encodeFloat+    exponent1 = exponent+    significand1 = significand+    scaleFloat1 = scaleFloat+    isNaN1 = isNaN+    isInfinite1 = isInfinite+    isDenormalized1 = isDenormalized+    isNegativeZero1 = isNegativeZero+    isIEEE1 = isIEEE+    atan21 = atan2+    succ1 = succ+    pred1 = pred+    toEnum1 = toEnum+    fromEnum1 = fromEnum+    enumFrom1 = enumFrom+    enumFromThen1 = enumFromThen+    enumFromTo1 = enumFromTo+    enumFromThenTo1 = enumFromThenTo+    minBound1 = minBound+    maxBound1 = maxBound++instance Mode Id where+    lift = Id+    Id a ^* b = Id (a * b)+    a *^ Id b = Id (a * b)+    Id a <+> Id b = Id (a + b)+    Id a <**> Id b = Id (a ** b)++instance Primal Id where+    primal (Id a) = a

@@ -0,0 +1,86 @@+{-# LANGUAGE TypeOperators, TemplateHaskell, ScopedTypeVariables, FlexibleContexts #-}+{-# OPTIONS_HADDOCK hide #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.AD.Internal.Jet+-- Copyright   :  (c) Edward Kmett 2010+-- License     :  BSD3+-- Maintainer  :  ekmett@gmail.com+-- Stability   :  experimental+-- Portability :  GHC only+--+-----------------------------------------------------------------------------++module Numeric.AD.Internal.Jet+    ( Jet(..)+    , headJet+    , tailJet+    , jet+    ) where++import Control.Applicative+import Data.Foldable+import Data.Traversable+import Data.Monoid+#if __GLASGOW_HASKELL__ < 704+import Data.Typeable (Typeable1(..), TyCon, mkTyCon, mkTyConApp)+#else+import Data.Typeable (Typeable1(..), TyCon, mkTyCon3, mkTyConApp)+#endif+import Control.Comonad.Cofree++infixl 3 :-++-- | A jet is a tower of all (higher order) partial derivatives of a function+data Jet f a = a :- Jet f (f a)++newtype Showable = Showable (Int -> String -> String)++instance Show Showable where+  showsPrec d (Showable f) = f d++showable :: Show a => a -> Showable+showable a = Showable (\d -> showsPrec d a)++-- Polymorphic recursion precludes 'Data' in its current form, as no Data1 class exists+-- Polymorphic recursion also breaks 'show' for 'Jet'!+-- factor Show1 out of Lifted?+instance (Functor f, Show (f Showable), Show a) => Show (Jet f a) where+  showsPrec d (a :- as) = showParen (d > 3) $+    showsPrec 4 a . showString " :- " . showsPrec 3 (fmap showable <$> as)++instance Functor f => Functor (Jet f) where+    fmap f (a :- as) = f a :- fmap (fmap f) as++instance Foldable f => Foldable (Jet f) where+    foldMap f (a :- as) = f a `mappend` foldMap (foldMap f) as++instance Traversable f => Traversable (Jet f) where+    traverse f (a :- as) = (:-) <$> f a <*> traverse (traverse f) as++tailJet :: Jet f a -> Jet f (f a)+tailJet (_ :- as) = as+{-# INLINE tailJet #-}++headJet :: Jet f a -> a+headJet (a :- _) = a+{-# INLINE headJet #-}++jet :: Functor f => Cofree f a -> Jet f a+jet (a :< as) = a :- dist (jet <$> as)+    where+        dist :: Functor f => f (Jet f a) -> Jet f (f a)+        dist x = (headJet <$> x) :- dist (tailJet <$> x)++instance Typeable1 f => Typeable1 (Jet f) where+    typeOf1 tfa = mkTyConApp jetTyCon [typeOf1 (undefined `asArgsType` tfa)]+        where asArgsType :: f a -> t f a -> f a+              asArgsType = const++jetTyCon :: TyCon+#if __GLASGOW_HASKELL__ < 704+jetTyCon = mkTyCon "Numeric.AD.Internal.Jet.Jet"+#else+jetTyCon = mkTyCon3 "ad" "Numeric.AD.Internal.Jet" "Jet"+#endif+{-# NOINLINE jetTyCon #-}

@@ -0,0 +1,281 @@+{-# LANGUAGE Rank2Types, TypeFamilies, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, FlexibleContexts, TemplateHaskell, UndecidableInstances, DeriveDataTypeable #-}+-- {-# OPTIONS_HADDOCK hide, prune #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.AD.Internal.Reverse+-- Copyright   :  (c) Edward Kmett 2010+-- License     :  BSD3+-- Maintainer  :  ekmett@gmail.com+-- Stability   :  experimental+-- Portability :  GHC only+--+-- Reverse-Mode Automatic Differentiation implementation details+--+-- 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.Internal.Reverse+    ( Reverse(..)+    , Tape(..)+    , partials+    , partialArray+    , partialMap+    , derivative+    , derivative'+    , Var(..)+    , bind+    , unbind+    , unbindMap+    , unbindWith+    , unbindMapWithDefault+    , vgrad, vgrad'+    , Grad(..)+    ) where++import Prelude hiding (mapM)+import Control.Applicative (Applicative(..),(<$>))+import Control.Monad.ST+import Control.Monad (forM_)+import Data.List (foldl', delete)+import Data.Array.ST+import Data.Array+import Data.IntMap (IntMap, fromListWith, findWithDefault, fromAscList, +                    updateLookupWithKey)+import qualified Data.IntSet as IS+import Data.Graph (graphFromEdges', Vertex, vertices, edges, transposeG, Graph)+import Data.Reify (reifyGraph, MuRef(..))+import qualified Data.Reify.Graph as Reified+import Data.Traversable (Traversable, mapM)+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++-- | 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)++-- | @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)++-- deriving instance (Data (Tape a (Reverse a)) => Data (Reverse a)++instance MuRef (Reverse a) where+    type DeRef (Reverse a) = Tape a++    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++instance Lifted Reverse => Mode Reverse where+    lift 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++    Reverse Zero <**> y                = lift (0 ** primal y)+    _            <**> Reverse Zero     = lift 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++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++instance Lifted Reverse => Jacobian Reverse where+    type D Reverse = Id++    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)++    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 _         _         (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)++    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 ''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 = go (fromAscList . assocs $ transposeG g) starters+  where starters = IS.toList $ foldl' (flip IS.delete)+                                      (IS.fromList $ vertices g)+                                      (map snd $ edges g)+        go _ [] = []+        go g' (n:ns) = let (g'',ns') = foldl' (uncurry (prune n)) (g',[]) (g!n)+                       in n : go g'' (ns'++ns)+        prune n g' acc m = let f _ = Just . delete n+                               (Just ns, g'') = updateLookupWithKey f m g'+                           in g'' `seq` (g'', if null (tail ns) then m:acc else acc)+++-- | This returns a list of contributions to the partials.+-- The variable ids returned in the list are likely /not/ unique!+partials :: Num a => AD Reverse a -> [(Int, a)]+partials (AD tape) = [ (ident, sensitivities ! ix) | (ix, Var _ ident) <- xs ]+    where+        Reified.Graph xs start = unsafePerformIO $ reifyGraph tape+        (g, vmap) = graphFromEdges' (edgeSet <$> filter nonConst xs)+        sensitivities = runSTArray $ do+            ss <- newArray (sbounds xs) 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!+        edgeSet (i, t) = (t, i, successors t)+        nonConst (_, Lift{}) = False+        nonConst _ = True+        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 Reverse a -> Array Int a+partialArray vbounds tape = accumArray (+) 0 vbounds (partials tape)+{-# INLINE partialArray #-}++-- | Return an 'IntMap' of sparse partials+partialMap :: Num a => AD Reverse 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')++-- | Used to mark variables for inspection during the reverse pass+class Primal v => Var v where+    var   :: a -> Int -> v a+    varId :: v a -> Int++instance Var Reverse where+    var a v = Reverse (Var a v)+    varId (Reverse (Var _ v)) = v+    varId _ = error "varId: not a Var"++instance Var (AD Reverse) where+    var a v = AD (var a v)+    varId (AD v) = varId v++bind :: (Traversable f, Var v) => f a -> (f (v a), (Int,Int))+bind xs = (r,(0,hi))+    where+        (r,hi) = runS (mapM freshVar xs) 0+        freshVar a = S (\s -> let s' = s + 1 in s' `seq` (var a s, s'))++unbind :: (Functor f, Var v)  => f (v a) -> Array Int a -> f a+unbind xs ys = fmap (\v -> ys ! varId v) xs++unbindWith :: (Functor f, Var v, Num a) => (a -> b -> c) -> f (v a) -> Array Int b -> f c+unbindWith f xs ys = fmap (\v -> f (primal v) (ys ! varId v)) xs++unbindMap :: (Functor f, Var v, Num a) => f (v a) -> IntMap a -> f a+unbindMap xs ys = fmap (\v -> findWithDefault 0 (varId v) ys) xs++unbindMapWithDefault :: (Functor f, Var v, Num a) => b -> (a -> b -> c) -> f (v a) -> IntMap b -> f c+unbindMapWithDefault z f xs ys = fmap (\v -> f (primal v) $ findWithDefault z (varId v) ys) xs++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+

@@ -0,0 +1,256 @@+{-# LANGUAGE BangPatterns, TemplateHaskell, TypeFamilies, TypeOperators, FlexibleContexts, UndecidableInstances, DeriveDataTypeable, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances #-}+{-# OPTIONS_GHC -fno-warn-name-shadowing #-}+module Numeric.AD.Internal.Sparse+    ( Index(..)+    , emptyIndex+    , addToIndex+    , indices+    , Sparse(..)+    , apply+    , vars+    , d, d', ds+    , skeleton+    , spartial+    , partial+    , vgrad+    , vgrad'+    , vgrads+    , Grad(..)+    , Grads(..)+    ) where++import Prelude hiding (lookup)+import Control.Applicative hiding ((<**>))+import Numeric.AD.Internal.Classes+import Control.Comonad.Cofree+import Numeric.AD.Internal.Types+import Data.Data+import Data.Typeable ()+import qualified Data.IntMap as IntMap+import Data.IntMap (IntMap, mapWithKey, unionWith, findWithDefault, toAscList, singleton, insertWith, lookup)+import Data.Traversable+import Language.Haskell.TH++newtype Index = Index (IntMap Int)++emptyIndex :: Index+emptyIndex = Index IntMap.empty+{-# INLINE emptyIndex #-}++addToIndex :: Int -> Index -> Index+addToIndex k (Index m) = Index (insertWith (+) k 1 m)+{-# INLINE addToIndex #-}++indices :: Index -> [Int]+indices (Index as) = uncurry (flip replicate) `concatMap` toAscList as+{-# INLINE indices #-}++-- | We only store partials in sorted order, so the map contained in a partial+-- will only contain partials with equal or greater keys to that of the map in+-- which it was found. This should be key for efficiently computing sparse hessians.+-- there are only (n + k - 1) choose k distinct nth partial derivatives of a+-- function with k inputs.+data Sparse a+  = Sparse !a (IntMap (Sparse a))+  | Zero+  deriving (Show, Data, Typeable)++-- | drop keys below a given value+dropMap :: Int -> IntMap a -> IntMap a+dropMap n = snd . IntMap.split (n - 1)+{-# INLINE dropMap #-}++times :: Num a => Sparse a -> Int -> Sparse a -> Sparse a+times Zero _ _ = Zero+times _ _ Zero = Zero+times (Sparse a as) n (Sparse b bs) = Sparse (a * b) $+    unionWith (<+>)+        (fmap (^* b) (dropMap n as))+        (fmap (a *^) (dropMap n bs))+{-# INLINE times #-}++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)+{-# INLINE vars #-}++apply :: (Traversable f, Num a) => (f (AD Sparse a) -> b) -> f a -> b+apply f = f . vars+{-# INLINE apply #-}++skeleton :: Traversable f => f a -> f Int+skeleton = snd . mapAccumL (\ !n _ -> (n + 1, n)) 0+{-# INLINE skeleton #-}++d :: (Traversable f, Num a) => f b -> AD Sparse a -> f a+d fs (AD Zero) = 0 <$ fs+d fs (AD (Sparse _ da)) = snd $ mapAccumL (\ !n _ -> (n + 1, maybe 0 primal $ lookup n da)) 0 fs+{-# INLINE d #-}++d' :: (Traversable f, Num a) => f a -> AD Sparse a -> (a, f a)+d' fs (AD Zero) = (0, 0 <$ fs)+d' fs (AD (Sparse a da)) = (a, snd $ mapAccumL (\ !n _ -> (n + 1, maybe 0 primal $ lookup n da)) 0 fs)+{-# INLINE d' #-}++ds :: (Traversable f, Num a) => f b -> AD Sparse a -> Cofree f a+ds fs (AD Zero) = r where r = 0 :< (r <$ fs)+ds fs (AD as@(Sparse a _)) = a :< (go emptyIndex <$> fns)+    where+        fns = skeleton fs+        -- go :: Index -> Int -> Cofree f a+        go ix i = partial (indices ix') as :< (go ix' <$> fns)+            where ix' = addToIndex i ix+{-# INLINE ds #-}++{-+vvars :: Num a => Vector a -> Vector (AD Sparse a)+vvars = Vector.imap (\n a -> AD $ Sparse a $ singleton n $ lift 1)+{-# INLINE vvars #-}++vapply :: Num a => (Vector (AD Sparse a) -> b) -> Vector a -> b+vapply f = f . vvars+{-# INLINE vapply #-}+++vd :: Num a => Int -> AD Sparse a -> Vector a+vd n (AD (Sparse _ da)) = Vector.generate n $ \i -> maybe 0 primal $ lookup i da+{-# INLINE vd #-}++vd' :: Num a => Int -> AD Sparse a -> (a, Vector a)+vd' n (AD (Sparse a da)) = (a , Vector.generate n $ \i -> maybe 0 primal $ lookup i da)+{-# INLINE vd' #-}++vds :: Num a => Int -> AD Sparse a -> Cofree Vector a+vds n (AD as@(Sparse a _)) = a :< Vector.generate n (go emptyIndex)+    where+        go ix i = partial (indices ix') as :< Vector.generate n (go ix')+            where ix' = addToIndex i ix+{-# INLINE vds #-}+-}++partial :: Num a => [Int] -> Sparse a -> a+partial []     (Sparse a _)  = a+partial (n:ns) (Sparse _ da) = partial ns $ findWithDefault (lift 0) n da+partial _      Zero          = 0+{-# INLINE partial #-}++spartial :: Num a => [Int] -> Sparse a -> Maybe a+spartial [] (Sparse a _) = Just a+spartial (n:ns) (Sparse _ da) = do+    a' <- lookup n da+    spartial ns a'+spartial _  Zero         = Nothing+{-# INLINE spartial #-}++instance Primal Sparse where+    primal (Sparse a _) = a+    primal Zero = 0++instance Lifted Sparse => Mode Sparse where+    lift a = Sparse a IntMap.empty+    zero = Zero+    Zero <**> y    = lift (0 ** primal y)+    _    <**> Zero = lift 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+    Zero <+> a = a+    a <+> Zero = a+    Sparse a as <+> Sparse b bs = Sparse (a + b) $ unionWith (<+>) as bs+    Zero        ^* _ = Zero+    Sparse a as ^* b = Sparse (a * b) $ fmap (^* b) as+    _ *^ Zero        = Zero+    a *^ Sparse b bs = Sparse (a * b) $ fmap (a *^) bs+    Zero        ^/ _ = Zero+    Sparse a as ^/ b = Sparse (a / b) $ fmap (^/ b) as++instance Lifted Sparse => Jacobian Sparse where+    type D Sparse = Sparse+    unary f _ Zero = lift (f 0)+    unary f dadb (Sparse pb bs) = Sparse (f pb) $ mapWithKey (times dadb) bs++    lift1 f _ Zero = lift (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 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 _    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) $+        unionWith (<+>)+            (mapWithKey (times dadb) db)+            (mapWithKey (times dadc) dc)++    lift2 f _  Zero             Zero = lift (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+        (dadb, dadc) = df b c+        da = unionWith (<+>)+            (mapWithKey (times dadb) db)+            (mapWithKey (times dadc) dc)++    lift2_ f _  Zero             Zero = lift (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+        (dadb, dadc) = df a b c+        a = Sparse (f pb pc) da+        da = unionWith (<+>)+            (mapWithKey (times dadb) db)+            (mapWithKey (times dadc) dc)++deriveLifted id $ conT ''Sparse+++class Num a => Grad i o o' a | i -> a o o', o -> a i o', o' -> a i o where+    pack :: i -> [AD Sparse a] -> AD Sparse a+    unpack :: ([a] -> [a]) -> o+    unpack' :: ([a] -> (a, [a])) -> o'++instance Num a => Grad (AD Sparse a) [a] (a, [a]) a where+    pack i _ = i+    unpack f = f []+    unpack' f = f []++instance Grad i o o' a => Grad (AD Sparse 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 = d as $ apply f as+{-# INLINE vgrad #-}++vgrad' :: Grad i o o' a => i -> o'+vgrad' i = unpack' (unsafeGrad' (pack i))+    where+        unsafeGrad' f as = d' as $ apply f as+{-# INLINE vgrad' #-}++class Num a => Grads i o a | i -> a o, o -> a i where+    packs :: i -> [AD Sparse a] -> AD Sparse a+    unpacks :: ([a] -> Cofree [] a) -> o++instance Num a => Grads (AD Sparse a) (Cofree [] a) a where+    packs i _ = i+    unpacks f = f []++instance Grads i o a => Grads (AD Sparse a -> i) (a -> o) a where+    packs f (a:as) = packs (f a) as+    packs _ [] = error "Grad.pack: logic error"+    unpacks f a = unpacks (f . (a:))++vgrads :: Grads i o a => i -> o+vgrads i = unpacks (unsafeGrads (packs i))+    where+        unsafeGrads f as = ds as $ apply f as+{-# INLINE vgrads #-}+

@@ -0,0 +1,140 @@+{-# LANGUAGE Rank2Types, TypeFamilies, FlexibleContexts, UndecidableInstances, TemplateHaskell, DeriveDataTypeable #-}+{-# OPTIONS_GHC -fno-warn-name-shadowing #-}+-- {-# OPTIONS_HADDOCK hide, prune #-}+-----------------------------------------------------------------------------+-- |+-- Module      : Numeric.AD.Tower.Internal+-- Copyright   : (c) Edward Kmett 2010+-- License     : BSD3+-- Maintainer  : ekmett@gmail.com+-- Stability   : experimental+-- Portability : GHC only+--+-----------------------------------------------------------------------------++module Numeric.AD.Internal.Tower+    ( Tower(..)+    , zeroPad+    , zeroPadF+    , transposePadF+    , d+    , d'+    , withD+    , tangents+    , bundle+    , apply+    , getADTower+    , tower+    ) where++import Prelude hiding (all)+import Control.Applicative hiding ((<**>))+import Data.Foldable+import Data.Data (Data)+import Data.Typeable (Typeable)+import Language.Haskell.TH+import Numeric.AD.Internal.Types+import Numeric.AD.Internal.Classes++-- | @Tower@ is an AD 'Mode' that calculates a tangent tower by forward AD, and provides fast 'diffsUU', 'diffsUF'+newtype Tower a = Tower { getTower :: [a] } deriving (Data, Typeable)++instance Show a => Show (Tower a) where+    showsPrec n (Tower as) = showParen (n > 10) $ showString "Tower " . showList as++-- Local combinators++zeroPad :: Num a => [a] -> [a]+zeroPad xs = xs ++ repeat 0+{-# INLINE zeroPad #-}++zeroPadF :: (Functor f, Num a) => [f a] -> [f a]+zeroPadF fxs@(fx:_) = fxs ++ repeat (const 0 <$> fx)+zeroPadF _ = error "zeroPadF :: empty list"+{-# INLINE zeroPadF #-}++transposePadF :: (Foldable f, Functor f) => a -> f [a] -> [f a]+transposePadF pad fx+    | all null fx = []+    | otherwise = fmap headPad fx : transposePadF pad (drop1 <$> fx)+    where+        headPad [] = pad+        headPad (x:_) = x+        drop1 (_:xs) = xs+        drop1 xs = xs++d :: Num a => [a] -> a+d (_:da:_) = da+d _ = 0+{-# INLINE d #-}++d' :: Num a => [a] -> (a, a)+d' (a:da:_) = (a, da)+d' (a:_)    = (a, 0)+d' _        = (0, 0)+{-# INLINE d' #-}++tangents :: Tower a -> Tower a+tangents (Tower []) = Tower []+tangents (Tower (_:xs)) = Tower xs+{-# INLINE tangents #-}++bundle :: a -> Tower a -> Tower a+bundle a (Tower as) = Tower (a:as)+{-# INLINE bundle #-}++withD :: (a, a) -> AD Tower a+withD (a, da) = AD (Tower [a,da])+{-# INLINE withD #-}++apply :: Num a => (AD Tower a -> b) -> a -> b+apply f a = f (AD (Tower [a,1]))+{-# INLINE apply #-}++getADTower :: AD Tower a -> [a]+getADTower (AD t) = getTower t+{-# INLINE getADTower #-}++tower :: [a] -> AD Tower a+tower as = AD (Tower as)++instance Primal Tower where+    primal (Tower (x:_)) = x+    primal _ = 0++instance Lifted Tower => Mode Tower where+    lift a = Tower [a]+    zero = Tower []+    Tower [] <**> y         = lift (0 ** primal y)+    _        <**> Tower []  = lift 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++    Tower [] <+> bs = bs+    as <+> Tower [] = as+    Tower (a:as) <+> Tower (b:bs) = Tower (c:cs)+        where+            c = a + b+            Tower cs = Tower as <+> Tower bs++    a *^ Tower bs = Tower (map (a*) bs)+    Tower as ^* b = Tower (map (*b) as)+    Tower as ^/ b = Tower (map (/b) as)++instance Lifted Tower => Jacobian Tower where+    type D Tower = Tower+    unary f dadb b = bundle (f (primal b)) (tangents b *! dadb)+    lift1 f df b   = bundle (f (primal b)) (tangents b *! df b)+    lift1_ f df b = a where+        a = bundle (f (primal b)) (tangents b *! df a b)++    binary f dadb dadc b c = bundle (f (primal b) (primal c)) (tangents b *! dadb +! tangents c *! dadc)+    lift2 f df b c = bundle (f (primal b) (primal c)) (tangents b *! dadb +! tangents c *! dadc) where+        (dadb, dadc) = df b c+    lift2_ f df b c = a where+        a0 = f (primal b) (primal c)+        da = tangents b *! dadb +! tangents c *! dadc+        a = bundle a0 da+        (dadb, dadc) = df a b c++deriveLifted id (conT ''Tower)

@@ -0,0 +1,59 @@+{-# LANGUAGE Rank2Types, GeneralizedNewtypeDeriving, TemplateHaskell, FlexibleContexts, FlexibleInstances, MultiParamTypeClasses, UndecidableInstances #-}+{-# OPTIONS_HADDOCK hide #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.AD.Internal.Types+-- Copyright   :  (c) Edward Kmett 2010+-- License     :  BSD3+-- Maintainer  :  ekmett@gmail.com+-- Stability   :  experimental+-- Portability :  GHC only+--+-----------------------------------------------------------------------------+module Numeric.AD.Internal.Types+    ( AD(..)+    ) where++import Data.Data (Data(..), mkDataType, DataType, mkConstr, Constr, constrIndex, Fixity(..))+import Data.Typeable (Typeable1(..), Typeable(..), TyCon, mkTyCon3, mkTyConApp, gcast1)+import Language.Haskell.TH+import Numeric.AD.Internal.Classes++-- | 'AD' serves as a common wrapper for different 'Mode' instances, exposing a traditional+-- numerical tower. Universal quantification is used to limit the actions in user code to+-- machinery that will return the same answers under all AD modes, allowing us to use modes+-- interchangeably as both the type level \"brand\" and dictionary, providing a common API.+newtype AD f a = AD { runAD :: f a } deriving (Iso (f a), Lifted, Mode, Primal)++-- > instance (Lifted f, Num a) => Num (AD f a)+-- etc.+let f = varT (mkName "f") in+    deriveNumeric+        (classP ''Lifted [f]:)+        (conT ''AD `appT` f)++instance Typeable1 f => Typeable1 (AD f) where+    typeOf1 tfa = mkTyConApp adTyCon [typeOf1 (undefined `asArgsType` tfa)]+        where asArgsType :: f a -> t f a -> f a+              asArgsType = const++adTyCon :: TyCon+adTyCon = mkTyCon3 "ad" "Numeric.AD.Internal.Types" "AD"+{-# NOINLINE adTyCon #-}++adConstr :: Constr+adConstr = mkConstr adDataType "AD" [] Prefix+{-# NOINLINE adConstr #-}++adDataType :: DataType+adDataType = mkDataType "Numeric.AD.Internal.Types.AD" [adConstr]+{-# NOINLINE adDataType #-}++instance (Typeable1 f, Typeable a, Data (f a), Data a) => Data (AD f a) where+    gfoldl f z (AD a) = z AD `f` a+    toConstr _ = adConstr+    gunfold k z c = case constrIndex c of+        1 -> k (z AD)+        _ -> error "gunfold"+    dataTypeOf _ = adDataType+    dataCast1 f = gcast1 f

@@ -0,0 +1,89 @@+{-# LANGUAGE Rank2Types #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.AD.Mode.Directed+-- Copyright   :  (c) Edward Kmett 2010+-- License     :  BSD3+-- Maintainer  :  ekmett@gmail.com+-- Stability   :  experimental+-- Portability :  GHC only+--+-- Allows the choice of AD 'Mode' to be specified at the term level for+-- benchmarking or more complicated usage patterns.+-----------------------------------------------------------------------------++module Numeric.AD.Mode.Directed+    (+    -- * Gradients+      grad+    , grad'+    -- * Jacobians+    , jacobian+    , jacobian'+    -- * Derivatives+    , diff+    , diff'+    -- * Exposed Types+    , Direction(..)+    ) where++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.Forward as F+import qualified Numeric.AD.Mode.Tower as T+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+    | Reverse+    | 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 Reverse = R.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' Reverse = R.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 Reverse = R.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' Reverse = R.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 Reverse = R.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' Reverse = R.grad'+grad' Tower   = F.grad' -- error "grad' Tower: unimplemented"+grad' Mixed   = M.grad'+{-# INLINE grad' #-}+

@@ -0,0 +1,161 @@+{-# LANGUAGE Rank2Types #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.AD.Mode.Forward+-- Copyright   :  (c) Edward Kmett 2010+-- License     :  BSD3+-- Maintainer  :  ekmett@gmail.com+-- Stability   :  experimental+-- Portability :  GHC only+--+-- Forward mode automatic differentiation+--+-----------------------------------------------------------------------------++module Numeric.AD.Mode.Forward+    (+    -- * Gradient+      grad+    , grad'+    , gradWith+    , gradWith'+    -- * Jacobian+    , jacobian+    , jacobian'+    , jacobianWith+    , jacobianWith'+    -- * Transposed Jacobian+    , jacobianT+    , jacobianWithT+    -- * Hessian Product+    , hessianProduct+    , hessianProduct'+    -- * Derivatives+    , diff+    , diff'+    , diffF+    , diffF'+    -- * Directional Derivatives+    , du+    , du'+    , duF+    , duF'+    ) where++import Data.Traversable (Traversable)+import Control.Applicative+import Numeric.AD.Types+import Numeric.AD.Internal.Classes+import Numeric.AD.Internal.Composition+import Numeric.AD.Internal.Forward++du :: (Functor f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f (a, a) -> a+du f = tangent . f . fmap (uncurry bundle)+{-# INLINE du #-}++du' :: (Functor f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f (a, a) -> (a, a)+du' f = unbundle . f . fmap (uncurry bundle)+{-# INLINE du' #-}++duF :: (Functor f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f (a, a) -> g a+duF f = fmap tangent . f . fmap (uncurry bundle)+{-# INLINE duF #-}++duF' :: (Functor f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f (a, a) -> g (a, a)+duF' f = fmap unbundle . f . fmap (uncurry bundle)+{-# INLINE duF' #-}++-- | The 'diff' function calculates the first derivative of a scalar-to-scalar function by forward-mode 'AD'+--+-- > diff sin == cos+diff :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> a+diff f a = tangent $ apply f a+{-# INLINE diff #-}++-- | The 'd'' function calculates the result and first derivative of scalar-to-scalar function by F'orward' 'AD'+-- +-- > d' sin == sin &&& cos+-- > d' f = f &&& d f+diff' :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a)+diff' f a = unbundle $ apply f a+{-# INLINE diff' #-}++-- | The 'diffF' function calculates the first derivative of scalar-to-nonscalar function by F'orward' 'AD'+diffF :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f a+diffF f a = tangent <$> apply f a+{-# INLINE diffF #-}++-- | The 'diffF'' function calculates the result and first derivative of a scalar-to-non-scalar function by F'orward' 'AD'+diffF' :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f (a, a)+diffF' f a = unbundle <$> apply f a+{-# INLINE diffF' #-}++-- | A fast, simple transposed Jacobian computed with forward-mode AD.+jacobianT :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> f (g a)+jacobianT f = bind (fmap tangent . f)+{-# INLINE jacobianT #-}++-- | A fast, simple transposed Jacobian computed with forward-mode AD.+jacobianWithT :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> f (g b)+jacobianWithT g f = bindWith g' f+    where g' a ga = g a . tangent <$> ga+{-# INLINE jacobianWithT #-}++jacobian :: (Traversable f, Traversable g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f a)+jacobian f as = transposeWith (const id) t p+    where+        (p, t) = bind' (fmap tangent . f) as+{-# INLINE jacobian #-}++jacobianWith :: (Traversable f, Traversable 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 = transposeWith (const id) t p+    where+        (p, t) = bindWith' g' f as+        g' a ga = g a . tangent <$> ga+{-# INLINE jacobianWith #-}++jacobian' :: (Traversable f, Traversable g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f a)+jacobian' f as = transposeWith row t p+    where+        (p, t) = bind' f as+        row x as' = (primal x, tangent <$> as')+{-# INLINE jacobian' #-}++jacobianWith' :: (Traversable f, Traversable 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 = transposeWith row t p+    where+        (p, t) = bindWith' g' f as+        row x as' = (primal x, as')+        g' a ga = g a . tangent <$> ga+{-# INLINE jacobianWith' #-}++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 #-}++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+        (b, bs) = bind' f as+{-# INLINE grad' #-}++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 #-}++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)+{-# INLINE gradWith' #-}++-- | Compute the product of a vector with the Hessian using forward-on-forward-mode AD. +hessianProduct :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f (a, a) -> f a+hessianProduct f = duF $ grad $ decomposeMode . f . fmap composeMode++-- | Compute the gradient and hessian product using forward-on-forward-mode AD. +hessianProduct' :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f (a, a) -> f (a, a)+hessianProduct' f = duF' $ grad $ decomposeMode . f . fmap composeMode++-- * Experimental++-- data f :> a = a :< f (f :> a)+-- gradients :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (f :> a)

@@ -0,0 +1,156 @@+{-# LANGUAGE Rank2Types, TemplateHaskell, BangPatterns, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, UndecidableInstances, ScopedTypeVariables #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.AD.Mode.Reverse+-- Copyright   :  (c) Edward Kmett 2010+-- License     :  BSD3+-- Maintainer  :  ekmett@gmail.com+-- 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.+--+-----------------------------------------------------------------------------++module Numeric.AD.Mode.Reverse+    (+    -- * 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.Reverse++-- | The 'grad' function calculates the gradient of a non-scalar-to-scalar function with 'Reverse' AD in a single pass.+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' AD in a single pass.+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' 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 :: (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''+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' #-}++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 'd'' function calculates the value and derivative, as a+-- pair, of a scalar-to-scalar function.+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' #-}++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 #-}++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 :: (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 :: (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))+

@@ -0,0 +1,121 @@+{-# LANGUAGE Rank2Types, BangPatterns #-}+-----------------------------------------------------------------------------+-- |+-- Module      : Numeric.AD.Mode.Sparse+-- Copyright   : (c) Edward Kmett 2010+-- License     : BSD3+-- Maintainer  : ekmett@gmail.com+-- Stability   : experimental+-- Portability : GHC only+--+-- Higher order derivatives via a \"dual number tower\".+--+-----------------------------------------------------------------------------++module Numeric.AD.Mode.Sparse+    (+    -- * Sparse Gradients+      grad+    , grad'+    , gradWith+    , gradWith'+    , grads++    -- * Sparse Jacobians (synonyms)+    , jacobian+    , jacobian'+    , jacobianWith+    , jacobianWith'+    , jacobians++    -- * Sparse Hessians+    , hessian+    , hessian'++    , hessianF+    , hessianF'++    -- * Unsafe gradients+    , vgrad+    , vgrads++    -- * Exposed Types+    , Grad+    , Grads+    ) where++import Control.Comonad+import Control.Applicative ((<$>))+import Data.Traversable+import Control.Comonad.Cofree+import Numeric.AD.Types+import Numeric.AD.Internal.Sparse+import Numeric.AD.Internal.Combinators++second :: (a -> b) -> (c, a) -> (c, b)+second g (a,b) = (a, g b)+{-# INLINE second #-}++grad :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f a+grad f as = d as $ apply f as+{-# INLINE grad #-}++grad' :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f a)+grad' f as = d' as $ apply f as+{-# INLINE grad' #-}++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 = zipWithT g as $ grad f as+{-# INLINE gradWith #-}++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 = second (zipWithT g as) $ grad' f as+{-# INLINE gradWith' #-}++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 = d as <$> apply f as+{-# INLINE jacobian #-}++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 = d' as <$> apply f as+{-# INLINE jacobian' #-}++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 = zipWithT g as <$> jacobian f as+{-# INLINE jacobianWith #-}++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 = second (zipWithT g as) <$> jacobian' f as+{-# INLINE jacobianWith' #-}++grads :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> Cofree f a+grads f as = ds as $ apply f as+{-# INLINE grads #-}++jacobians :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (Cofree f a)+jacobians f as = ds as <$> apply f as+{-# INLINE jacobians #-}++d2 :: Functor f => Cofree f a -> f (f a)+d2 = headJet . tailJet . tailJet . jet+{-# INLINE d2 #-}++d2' :: Functor f => Cofree f a -> (a, f (a, f a))+d2' (a :< as) = (a, fmap (\(da :< das) -> (da, extract <$> das)) as)+{-# INLINE d2' #-}++hessian :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f (f a)+hessian f as = d2 $ grads f as+{-# INLINE hessian #-}++hessian' :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f (a, f a))+hessian' f as = d2' $ grads f as+{-# INLINE hessian' #-}++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 as = d2 <$> jacobians f as+{-# INLINE hessianF #-}++hessianF' :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f (a, f a))+hessianF' f as = d2' <$> jacobians f as+{-# INLINE hessianF' #-}

@@ -0,0 +1,123 @@+{-# LANGUAGE Rank2Types, BangPatterns #-}+-----------------------------------------------------------------------------+-- |+-- Module      : Numeric.AD.Mode.Tower+-- Copyright   : (c) Edward Kmett 2010+-- License     : BSD3+-- Maintainer  : ekmett@gmail.com+-- Stability   : experimental+-- Portability : GHC only+--+-- Higher order derivatives via a \"dual number tower\".+--+-----------------------------------------------------------------------------++module Numeric.AD.Mode.Tower+    (+    -- * Taylor Series+      taylor+    , taylor0+    -- * Maclaurin Series+    , maclaurin+    , maclaurin0+    -- * Derivatives+    , diff    -- first derivative of (a -> a)+    , diff'   -- answer and first derivative of (a -> a)+    , diffs   -- answer and all derivatives of (a -> a)+    , diffs0  -- zero padded derivatives of (a -> a)+    , diffsF  -- answer and all derivatives of (a -> f a)+    , diffs0F -- zero padded derivatives of (a -> f a)+    -- * Directional Derivatives+    , du      -- directional derivative of (a -> a)+    , du'     -- answer and directional derivative of (a -> a)+    , dus     -- answer and all directional derivatives of (a -> a)+    , dus0    -- answer and all zero padded directional derivatives of (a -> a)+    , duF     -- directional derivative of (a -> f a)+    , duF'    -- answer and directional derivative of (a -> f a)+    , dusF    -- answer and all directional derivatives of (a -> f a)+    , dus0F   -- answer and all zero padded directional derivatives of (a -> a)+    ) where++import Control.Applicative ((<$>))+import Numeric.AD.Types+import Numeric.AD.Internal.Tower++diffs :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]+diffs f a = getADTower $ apply f a+{-# INLINE diffs #-}++diffs0 :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]+diffs0 f a = zeroPad (diffs f a)+{-# INLINE diffs0 #-}++diffsF :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f [a]+diffsF f a = getADTower <$> apply f a+{-# INLINE diffsF #-}++diffs0F :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f [a]+diffs0F f a = (zeroPad . getADTower) <$> apply f a+{-# INLINE diffs0F #-}++taylor :: Fractional a => (forall s. Mode s => AD s a -> AD s a) -> a -> a -> [a]+taylor f x dx = go 1 1 (diffs f x)+    where+        go !n !acc (a:as) = a * acc : go (n + 1) (acc * dx / n) as+        go _ _ [] = []++taylor0 :: Fractional a => (forall s. Mode s => AD s a -> AD s a) -> a -> a -> [a]+taylor0 f x dx = zeroPad (taylor f x dx)+{-# INLINE taylor0 #-}++maclaurin :: Fractional a => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]+maclaurin f = taylor f 0+{-# INLINE maclaurin #-}++maclaurin0 :: Fractional a => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]+maclaurin0 f = taylor0 f 0+{-# INLINE maclaurin0 #-}++diff :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> a+diff f = d . diffs f+{-# INLINE diff #-}++diff' :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a)+diff' f = d' . diffs f+{-# INLINE diff' #-}++du :: (Functor f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f (a, a) -> a+du f = d . getADTower . f . fmap withD+{-# INLINE du #-}++du' :: (Functor f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f (a, a) -> (a, a)+du' f = d' . getADTower . f . fmap withD+{-# INLINE du' #-}++duF :: (Functor f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f (a, a) -> g a+duF f = fmap (d . getADTower) . f . fmap withD+{-# INLINE duF #-}++duF' :: (Functor f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f (a, a) -> g (a, a)+duF' f = fmap (d' . getADTower) . f . fmap withD+{-# INLINE duF' #-}++dus :: (Functor f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f [a] -> [a]+dus f = getADTower . f . fmap tower+{-# INLINE dus #-}++dus0 :: (Functor f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f [a] -> [a]+dus0 f = zeroPad . getADTower . f . fmap tower+{-# INLINE dus0 #-}++dusF :: (Functor f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f [a] -> g [a]+dusF f = fmap getADTower . f . fmap tower+{-# INLINE dusF #-}++dus0F :: (Functor f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f [a] -> g [a]+dus0F f = fmap getADTower . f . fmap tower+{-# INLINE dus0F #-}++-- TODO: higher order gradients+-- data f :> a = a :< f (f :> a)+-- gradients  :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f :> a+-- gradientsF, jacobians :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f :> a)+-- gradientsM :: (Traversable f, Monad m, Num a) => FF f m a -> f a -> m (f :> a)

@@ -0,0 +1,108 @@+{-# LANGUAGE Rank2Types, BangPatterns, ScopedTypeVariables #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.AD.Newton+-- Copyright   :  (c) Edward Kmett 2010+-- License     :  BSD3+-- Maintainer  :  ekmett@gmail.com+-- Stability   :  experimental+-- Portability :  GHC only+--+-----------------------------------------------------------------------------++module Numeric.AD.Newton+    (+    -- * Newton's Method (Forward AD)+      findZero+    , inverse+    , fixedPoint+    , extremum+    -- * Gradient Ascent/Descent (Reverse AD)+    , gradientDescent+    , gradientAscent+    ) where++import Prelude hiding (all)+import Data.Foldable (all)+import Data.Traversable (Traversable)+import Numeric.AD.Types+import Numeric.AD.Mode.Forward (diff, diff')+import Numeric.AD.Mode.Reverse (gradWith')+import Numeric.AD.Internal.Composition++-- | The 'findZero' function finds a zero of a scalar function using+-- Newton's method; its output is a stream of increasingly accurate+-- results.  (Modulo the usual caveats.)+--+-- Examples:+--+--  > take 10 $ findZero (\\x->x^2-4) 1  -- converge to 2.0+--+--  > module Data.Complex+--  > take 10 $ findZero ((+1).(^2)) (1 :+ 1)  -- converge to (0 :+ 1)@+--+findZero :: (Fractional a, Eq a) => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]+findZero f = go+    where+        go x = x : if y == 0 then [] else go (x - y/y')+            where+                (y,y') = diff' f x+{-# INLINE findZero #-}++-- | The 'inverseNewton' function inverts a scalar function using+-- Newton's method; its output is a stream of increasingly accurate+-- results.  (Modulo the usual caveats.)+--+-- Example:+--+-- > take 10 $ inverseNewton sqrt 1 (sqrt 10)  -- converges to 10+--+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+{-# INLINE inverse  #-}++-- | The 'fixedPoint' function find a fixedpoint of a scalar+-- function using Newton's method; its output is a stream of+-- increasingly accurate results.  (Modulo the usual caveats.)+--+-- > take 10 $ fixedPoint cos 1 -- converges to 0.7390851332151607+fixedPoint :: (Fractional a, Eq a) => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]+fixedPoint f = findZero (\x -> f x - x)+{-# INLINE fixedPoint #-}++-- | The 'extremum' function finds an extremum of a scalar+-- function using Newton's method; produces a stream of increasingly+-- accurate results.  (Modulo the usual caveats.)+--+-- > take 10 $ extremum cos 1 -- convert to 0+extremum :: (Fractional a, Eq a) => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]+extremum f = findZero (diff (decomposeMode . f . composeMode))+{-# INLINE extremum #-}++-- | The 'gradientDescent' function performs a multivariate+-- optimization, based on the naive-gradient-descent in the file+-- @stalingrad\/examples\/flow-tests\/pre-saddle-1a.vlad@ from the+-- VLAD compiler Stalingrad sources.  Its output is a stream of+-- increasingly accurate results.  (Modulo the usual caveats.)+--+-- It uses reverse mode automatic differentiation to compute the gradient.+gradientDescent :: (Traversable f, Fractional a, Ord a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> [f a]+gradientDescent f x0 = go x0 fx0 xgx0 0.1 (0 :: Int)+    where+        (fx0, xgx0) = gradWith' (,) f x0+        go x fx xgx !eta !i+            | eta == 0     = [] -- step size is 0+            | fx1 > fx     = go x fx xgx (eta/2) 0 -- we stepped too far+            | zeroGrad xgx = [] -- gradient is 0+            | otherwise    = x1 : if i == 10+                                  then go x1 fx1 xgx1 (eta*2) 0+                                  else go x1 fx1 xgx1 eta (i+1)+            where+                zeroGrad = all (\(_,g) -> g == 0)+                x1 = fmap (\(xi,gxi) -> xi - eta * gxi) xgx+                (fx1, xgx1) = gradWith' (,) f x1+{-# INLINE gradientDescent #-}++gradientAscent :: (Traversable f, Fractional a, Ord a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> [f a]+gradientAscent f = gradientDescent (negate . f)+{-# INLINE gradientAscent #-}

@@ -0,0 +1,49 @@+{-# LANGUAGE Rank2Types #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.AD.Types+-- Copyright   :  (c) Edward Kmett 2010+-- License     :  BSD3+-- Maintainer  :  ekmett@gmail.com+-- Stability   :  experimental+-- Portability :  GHC only+--+-----------------------------------------------------------------------------++module Numeric.AD.Types+    (+    -- * AD modes+      Mode(..)+    -- * AD variables+    , AD(..)+    -- * Jets+    , Jet(..)+    , headJet+    , tailJet+    , jet+    -- * Apply functions that use 'lift'+    , lowerUU, lowerUF, lowerFU, lowerFF+    ) where++import Numeric.AD.Internal.Identity+import Numeric.AD.Internal.Types+import Numeric.AD.Internal.Jet+import Numeric.AD.Internal.Classes++-- these exploit the 'magic' that is probed to avoid the need for Functor, etc.++lowerUU :: (forall s. Mode s => AD s a -> AD s a) -> a -> a+lowerUU f = unprobe . f . probe+{-# INLINE lowerUU #-}++lowerUF :: (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f a+lowerUF f = unprobed . f . probe+{-# INLINE lowerUF #-}++lowerFU :: (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> a+lowerFU f = unprobe . f . probed+{-# INLINE lowerFU #-}++lowerFF :: (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g a+lowerFF f = unprobed . f . probed+{-# INLINE lowerFF #-}

@@ -0,0 +1,29 @@+{-# LANGUAGE Rank2Types, TemplateHaskell, BangPatterns, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, UndecidableInstances, ScopedTypeVariables #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.AD.Variadic+-- Copyright   :  (c) Edward Kmett 2010-2012+-- License     :  BSD3+-- Maintainer  :  ekmett@gmail.com+-- Stability   :  experimental+-- Portability :  non-portable+--+-- Variadic combinators for variadic mixed-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-mode variadic gradient+      Grad , vgrad, vgrad'+    -- * Sparse forward mode variadic jet+    , Grads, vgrads+    ) where++import Numeric.AD.Variadic.Reverse+import Numeric.AD.Variadic.Sparse (Grads, vgrads)

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

@@ -0,0 +1,27 @@+{-# LANGUAGE Rank2Types, TemplateHaskell, BangPatterns, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, UndecidableInstances, ScopedTypeVariables #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.AD.Variadic.Sparse+-- Copyright   :  (c) Edward Kmett 2010-2012+-- License     :  BSD3+-- Maintainer  :  ekmett@gmail.com+-- Stability   :  experimental+-- Portability :  non-portable+--+-- Variadic combinators for sparse forward 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.Sparse+    (+    -- * Unsafe Variadic Gradient+      Grad , vgrad, vgrad'+    , Grads, vgrads+    ) where++import Numeric.AD.Internal.Sparse