diff --git a/Numeric/AD.hs b/Numeric/AD.hs deleted file mode 100644 --- a/Numeric/AD.hs +++ /dev/null @@ -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) diff --git a/Numeric/AD/Classes.hs b/Numeric/AD/Classes.hs deleted file mode 100644 --- a/Numeric/AD/Classes.hs +++ /dev/null @@ -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 diff --git a/Numeric/AD/Halley.hs b/Numeric/AD/Halley.hs deleted file mode 100644 --- a/Numeric/AD/Halley.hs +++ /dev/null @@ -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 #-} - diff --git a/Numeric/AD/Internal/Classes.hs b/Numeric/AD/Internal/Classes.hs deleted file mode 100644 --- a/Numeric/AD/Internal/Classes.hs +++ /dev/null @@ -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' diff --git a/Numeric/AD/Internal/Combinators.hs b/Numeric/AD/Internal/Combinators.hs deleted file mode 100644 --- a/Numeric/AD/Internal/Combinators.hs +++ /dev/null @@ -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) diff --git a/Numeric/AD/Internal/Composition.hs b/Numeric/AD/Internal/Composition.hs deleted file mode 100644 --- a/Numeric/AD/Internal/Composition.hs +++ /dev/null @@ -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 - diff --git a/Numeric/AD/Internal/Dense.hs b/Numeric/AD/Internal/Dense.hs deleted file mode 100644 --- a/Numeric/AD/Internal/Dense.hs +++ /dev/null @@ -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) diff --git a/Numeric/AD/Internal/Forward.hs b/Numeric/AD/Internal/Forward.hs deleted file mode 100644 --- a/Numeric/AD/Internal/Forward.hs +++ /dev/null @@ -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 - diff --git a/Numeric/AD/Internal/Identity.hs b/Numeric/AD/Internal/Identity.hs deleted file mode 100644 --- a/Numeric/AD/Internal/Identity.hs +++ /dev/null @@ -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 diff --git a/Numeric/AD/Internal/Reverse.hs b/Numeric/AD/Internal/Reverse.hs deleted file mode 100644 --- a/Numeric/AD/Internal/Reverse.hs +++ /dev/null @@ -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 - diff --git a/Numeric/AD/Internal/Sparse.hs b/Numeric/AD/Internal/Sparse.hs deleted file mode 100644 --- a/Numeric/AD/Internal/Sparse.hs +++ /dev/null @@ -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 #-} - diff --git a/Numeric/AD/Internal/Tensors.hs b/Numeric/AD/Internal/Tensors.hs deleted file mode 100644 --- a/Numeric/AD/Internal/Tensors.hs +++ /dev/null @@ -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 #-} diff --git a/Numeric/AD/Internal/Tower.hs b/Numeric/AD/Internal/Tower.hs deleted file mode 100644 --- a/Numeric/AD/Internal/Tower.hs +++ /dev/null @@ -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) diff --git a/Numeric/AD/Internal/Types.hs b/Numeric/AD/Internal/Types.hs deleted file mode 100644 --- a/Numeric/AD/Internal/Types.hs +++ /dev/null @@ -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 diff --git a/Numeric/AD/Mode/Directed.hs b/Numeric/AD/Mode/Directed.hs deleted file mode 100644 --- a/Numeric/AD/Mode/Directed.hs +++ /dev/null @@ -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' #-} - diff --git a/Numeric/AD/Mode/Forward.hs b/Numeric/AD/Mode/Forward.hs deleted file mode 100644 --- a/Numeric/AD/Mode/Forward.hs +++ /dev/null @@ -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) diff --git a/Numeric/AD/Mode/Mixed.hs b/Numeric/AD/Mode/Mixed.hs deleted file mode 100644 --- a/Numeric/AD/Mode/Mixed.hs +++ /dev/null @@ -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 , 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 diff --git a/Numeric/AD/Mode/Reverse.hs b/Numeric/AD/Mode/Reverse.hs deleted file mode 100644 --- a/Numeric/AD/Mode/Reverse.hs +++ /dev/null @@ -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)) - diff --git a/Numeric/AD/Mode/Sparse.hs b/Numeric/AD/Mode/Sparse.hs deleted file mode 100644 --- a/Numeric/AD/Mode/Sparse.hs +++ /dev/null @@ -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' #-} - diff --git a/Numeric/AD/Mode/Tower.hs b/Numeric/AD/Mode/Tower.hs deleted file mode 100644 --- a/Numeric/AD/Mode/Tower.hs +++ /dev/null @@ -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) diff --git a/Numeric/AD/Newton.hs b/Numeric/AD/Newton.hs deleted file mode 100644 --- a/Numeric/AD/Newton.hs +++ /dev/null @@ -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 #-} diff --git a/Numeric/AD/Types.hs b/Numeric/AD/Types.hs deleted file mode 100644 --- a/Numeric/AD/Types.hs +++ /dev/null @@ -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 #-} diff --git a/ad.cabal b/ad.cabal --- a/ad.cabal +++ b/ad.cabal @@ -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 diff --git a/src/Numeric/AD.hs b/src/Numeric/AD.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AD.hs @@ -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 , 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 diff --git a/src/Numeric/AD/Halley.hs b/src/Numeric/AD/Halley.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AD/Halley.hs @@ -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 #-} + diff --git a/src/Numeric/AD/Internal/Classes.hs b/src/Numeric/AD/Internal/Classes.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AD/Internal/Classes.hs @@ -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' diff --git a/src/Numeric/AD/Internal/Combinators.hs b/src/Numeric/AD/Internal/Combinators.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AD/Internal/Combinators.hs @@ -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) diff --git a/src/Numeric/AD/Internal/Composition.hs b/src/Numeric/AD/Internal/Composition.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AD/Internal/Composition.hs @@ -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 + diff --git a/src/Numeric/AD/Internal/Dense.hs b/src/Numeric/AD/Internal/Dense.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AD/Internal/Dense.hs @@ -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) diff --git a/src/Numeric/AD/Internal/Forward.hs b/src/Numeric/AD/Internal/Forward.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AD/Internal/Forward.hs @@ -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 + diff --git a/src/Numeric/AD/Internal/Identity.hs b/src/Numeric/AD/Internal/Identity.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AD/Internal/Identity.hs @@ -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 diff --git a/src/Numeric/AD/Internal/Jet.hs b/src/Numeric/AD/Internal/Jet.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AD/Internal/Jet.hs @@ -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 #-} diff --git a/src/Numeric/AD/Internal/Reverse.hs b/src/Numeric/AD/Internal/Reverse.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AD/Internal/Reverse.hs @@ -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 + diff --git a/src/Numeric/AD/Internal/Sparse.hs b/src/Numeric/AD/Internal/Sparse.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AD/Internal/Sparse.hs @@ -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 #-} + diff --git a/src/Numeric/AD/Internal/Tower.hs b/src/Numeric/AD/Internal/Tower.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AD/Internal/Tower.hs @@ -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) diff --git a/src/Numeric/AD/Internal/Types.hs b/src/Numeric/AD/Internal/Types.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AD/Internal/Types.hs @@ -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 diff --git a/src/Numeric/AD/Mode/Directed.hs b/src/Numeric/AD/Mode/Directed.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AD/Mode/Directed.hs @@ -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' #-} + diff --git a/src/Numeric/AD/Mode/Forward.hs b/src/Numeric/AD/Mode/Forward.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AD/Mode/Forward.hs @@ -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) diff --git a/src/Numeric/AD/Mode/Reverse.hs b/src/Numeric/AD/Mode/Reverse.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AD/Mode/Reverse.hs @@ -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)) + diff --git a/src/Numeric/AD/Mode/Sparse.hs b/src/Numeric/AD/Mode/Sparse.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AD/Mode/Sparse.hs @@ -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' #-} diff --git a/src/Numeric/AD/Mode/Tower.hs b/src/Numeric/AD/Mode/Tower.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AD/Mode/Tower.hs @@ -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) diff --git a/src/Numeric/AD/Newton.hs b/src/Numeric/AD/Newton.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AD/Newton.hs @@ -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 #-} diff --git a/src/Numeric/AD/Types.hs b/src/Numeric/AD/Types.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AD/Types.hs @@ -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 #-} diff --git a/src/Numeric/AD/Variadic.hs b/src/Numeric/AD/Variadic.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AD/Variadic.hs @@ -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) diff --git a/src/Numeric/AD/Variadic/Reverse.hs b/src/Numeric/AD/Variadic/Reverse.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AD/Variadic/Reverse.hs @@ -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 diff --git a/src/Numeric/AD/Variadic/Sparse.hs b/src/Numeric/AD/Variadic/Sparse.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AD/Variadic/Sparse.hs @@ -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