ad (empty) → 0.12

raw patch · 15 files changed

+1513/−0 lines, 15 filesdep +arraydep +basedep +containerssetup-changed

Dependencies added: array, base, containers, data-reify, template-haskell

Files

LICENSE +30/−0
Numeric/AD.hs +85/−0
Numeric/AD/Classes.hs +321/−0
Numeric/AD/Directed.hs +103/−0
Numeric/AD/Forward.hs +100/−0
Numeric/AD/Internal.hs +128/−0
Numeric/AD/Internal/Forward.hs +112/−0
Numeric/AD/Internal/Reverse.hs +202/−0
Numeric/AD/Internal/Tower.hs +101/−0
Numeric/AD/Newton.hs +66/−0
Numeric/AD/Reverse.hs +115/−0
Numeric/AD/Tower.hs +98/−0
Setup.hs +2/−0
TODO +13/−0
ad.cabal +37/−0

+ LICENSE view

@@ -0,0 +1,30 @@+Copyright (c) 2010, Edward Kmett+All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are+met:++    * Redistributions of source code must retain the above copyright+      notice, this list of conditions and the following disclaimer.++    * Redistributions in binary form must reproduce the above+      copyright notice, this list of conditions and the following+      disclaimer in the documentation and/or other materials provided+      with the distribution.++    * Neither the name of Edward Kmett nor the names of other+      contributors may be used to endorse or promote products derived+      from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

+ Numeric/AD.hs view

@@ -0,0 +1,85 @@+{-# LANGUAGE Rank2Types, TypeFamilies #-}+-----------------------------------------------------------------------------+-- |+-- 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.+-----------------------------------------------------------------------------++module Numeric.AD +    ( +    -- * Gradients+      grad, grad2++    -- * Jacobians+    , jacobian, jacobian2++    -- * Synonyms+    , diff+    , diff2+    , diffs+    , diffs0++    -- * Derivatives (Forward)+    , diffUU+    , diffUF++    , diff2UU+    , diff2UF++    -- * Derivatives (Reverse)+    , diffFU+    , diff2FU++    -- * Derivatives (Tower)+    , diffsUU+    , diffsUF++    , diffs0UU+    , diffs0UF++    -- * Taylor Series (Tower)+    , taylor+    , taylor0++    -- * Exposed Types+    , AD(..)+    , Mode(..)+    ) where++import Data.Traversable (Traversable)+import Data.Foldable (Foldable, foldr')+import Control.Applicative+import Numeric.AD.Classes  (Mode(..))+import Numeric.AD.Internal (AD(..), probe, unprobe)+import Numeric.AD.Forward  (diff, diffUU, diff2, diff2UU, diffUF, diff2UF)+import Numeric.AD.Tower    (diffsUU, diffs0UU , diffsUF, diffs0UF , diffs, diffs0, taylor, taylor0) +import Numeric.AD.Reverse  (diffFU, diff2FU, grad, grad2)++import qualified Numeric.AD.Forward as Forward+import qualified Numeric.AD.Reverse as Reverse++-- | 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+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 bs = snd <$> jacobian2 f bs+{-# INLINE jacobian #-}++-- | Calculate the answer and Jacobian of a non-scalar-to-non-scalar function, automatically choosing between forward- and reverse- mode AD based on the relative, number of inputs and outputs. If you need to support functions where the output is only a 'Functor', consider using 'jacobianT' from "Numeric.AD.Forward" or 'jacobian2' from "Numeric.AD.Reverse" directly.+jacobian2 :: (Traversable f, Traversable g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f a)+jacobian2 f bs | n == 0    = fmap (\x -> (unprobe x, bs)) as+               | n > m     = Reverse.jacobian2 f bs+               | otherwise = Forward.jacobian2 f bs+    where+        as = f (probe <$> bs)+        n = size bs+        m = size as+        size :: Foldable f => f a -> Int+        size = foldr' (\_ b -> 1 + b) 0 +{-# INLINE jacobian2 #-}

+ Numeric/AD/Classes.hs view

@@ -0,0 +1,321 @@+{-# LANGUAGE Rank2Types, TypeFamilies, FlexibleInstances, MultiParamTypeClasses, FlexibleContexts, FunctionalDependencies, UndecidableInstances, GeneralizedNewtypeDeriving, TemplateHaskell #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.AD.Classes+-- Copyright   :  (c) Edward Kmett 2010+-- License     :  BSD3+-- Maintainer  :  ekmett@gmail.com+-- Stability   :  experimental+-- Portability :  GHC only +--+-----------------------------------------------------------------------------++module Numeric.AD.Classes+    ( +    -- * AD modes+      Mode(..) +    , one+    -- * Automatically Deriving AD+    , Jacobian(..)+    , Primal(..)+    , deriveLifted+    , deriveNumeric+    , Lifted(..)+    ) where++import Control.Applicative+import Data.Char+import Language.Haskell.TH+-- import Text.Show++infixl 8 **!+infixl 7 *!, /!, ^*, *^, ^/+infixl 6 +!, -!, <+> +infix 4 ==!++-- | Higher-order versions of the stock numerical methods.+class Lifted t where+-- class Show1 t where+    showsPrec1 :: Show a => Int -> t a -> ShowS+--    show1 :: Show a => t a -> String+--    showList1 :: Show a => [t a] -> String -> String++-- class Eq1 t where+    (==!) :: (Num a, Eq a) => t a -> t a -> Bool+    -- (/=!) :: (Num a, Eq a) => t a -> t a -> Bool++-- class Eq1 => Ord1 t where+    compare1 :: (Num a, Ord a) => t a -> t a -> Ordering+    -- (<!) :: (Num a, Ord a) => t a -> t a -> Bool+    -- (>=!) :: (Num a, Ord a) => t a -> t a -> Bool+    -- (>!) :: (Num a, Ord a) => t a -> t a -> Bool+    -- (<=!) :: (Num a, Ord a) => t a -> t a -> Bool+    -- min1 :: (Num a, Ord a) => t a -> t a -> t a+    -- max1 :: (Num a, Ord a) => t a -> t a -> t a++-- class (Show1 t, Eq t) => Num1 t where+    fromInteger1 :: Num a => Integer -> t a+    (+!),(-!),(*!) :: Num a => t a -> t a -> t a+    negate1, abs1, signum1 :: Num a => t a -> t a ++-- class Num1 t => Fractional1 t where+    (/!) :: Fractional a => t a -> t a -> t a+    recip1 :: Fractional a => t a -> t a+    fromRational1 :: Fractional a => Rational -> t a ++-- class (Num1 t, Ord1 t) => Real1 t +    toRational1 :: Real a => t a -> Rational -- unsafe++-- class Fractional1 t => Floating1 t +    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++-- class (Real1 t, Fractional1 t) => RealFrac1 t where+    properFraction1 :: (RealFrac a, Integral b) => t a -> (b, t a)++    truncate1, round1, ceiling1, floor1 :: (RealFrac a, Integral b) => t a -> b++-- class (RealFrac1 t, Floating1 t) => RealFloat1 t where+    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++-- class Enum1 t where+    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]++-- class Bounded1 t where+    minBound1 :: (Num a, Bounded a) => t a +    maxBound1 :: (Num a, Bounded a) => t a++class Lifted t => Mode t where++    -- | 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 ++    -- | > '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' = 'lift' 1+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++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 #-}++on :: (a -> a -> c) -> (b -> a) -> b -> b -> c+on f g a b = f (g a) (g b)++discrete1 :: (Primal t, Num a) => (a -> c) -> t a -> c+discrete1 f x = f (primal x)++discrete2 :: (Primal t, Num a) => (a -> a -> c) -> t a -> t a -> c+discrete2 f x y = f (primal x) (primal y)++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)++-- | @'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 Type -> Q [Dec]+deriveLifted _t = [d| +    instance Lifted $_t where+        (==!)         = (==) `on` primal+        compare1      = compare `on` primal+        maxBound1     = lift maxBound+        minBound1     = lift minBound+        showsPrec1    = showsPrec+        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+        (/!)          = lift2 (/) $ \x y -> (recip1 y, x)+        recip1        = lift1 recip (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))+        (**!)     = lift2_ (**) (\z x y -> (y *! z /! x, z *! log1 x)) -- error at 0 ** n +        sin1      = lift1 sin cos1+        cos1      = lift1 cos $ \x -> negate1 (sin1 x)+        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+    ClassI (ClassD _ _ _ _ ds) <- reify ''Lifted+    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 Type -> Q Type -> Q [Dec]+deriveNumeric t' t = do+    members <- liftedMembers+    let keep n = nameBase n `elem` members+    xs <- lowerInstance keep (classP ''Num [varA]:) t t' `mapM` [''Enum, ''Eq, ''Ord, ''Bounded]+    ys <- lowerInstance keep id                     t t' `mapM` [''Show, ''Num, ''Fractional, ''Floating, ''RealFloat,''RealFrac, ''Real]+    return (xs ++ ys)++lowerInstance :: (Name -> Bool) -> ([Q Pred] -> [Q Pred]) -> Q Type -> Q Type -> Name -> Q Dec+lowerInstance p f t t' n = do+    ClassI (ClassD _ _ _ _ ds) <- reify n+    instanceD (cxt (f [classP ''Lifted [t], 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'

+ Numeric/AD/Directed.hs view

@@ -0,0 +1,103 @@+{-# LANGUAGE Rank2Types #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.AD.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.Directed+    ( +    -- * Derivatives+      diffUU+    , diff2UU+    -- * Common access patterns+    , diff+    , diff2+    -- * Jacobians+    , jacobian+    , jacobian2+    -- * Gradients+    , grad+    , grad2+    -- * Exposed Types+    , Direction(..)+    , Mode(..)+    , AD(..)+    ) where++import Prelude hiding (reverse)+import Numeric.AD.Classes+import Numeric.AD.Internal+import Data.Traversable (Traversable)+import qualified Numeric.AD.Reverse as R+import qualified Numeric.AD.Forward as F+import qualified Numeric.AD.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)++diffUU :: Num a => Direction -> (forall s. Mode s => AD s a -> AD s a) -> a -> a+diffUU Forward = F.diffUU+diffUU Reverse = R.diffUU+diffUU Tower = T.diffUU+diffUU Mixed = F.diffUU+{-# INLINE diffUU #-}++diff2UU :: Num a => Direction -> (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a) +diff2UU Forward = F.diff2UU+diff2UU Reverse = R.diff2UU+diff2UU Tower = T.diff2UU+diff2UU Mixed = F.diff2UU+{-# INLINE diff2UU #-}++diff :: Num a => Direction -> (forall s. Mode s => AD s a -> AD s a) -> a -> a+diff = diffUU+{-# INLINE diff #-}++diff2 :: Num a => Direction -> (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a) +diff2 = diff2UU+{-# INLINE diff2 #-}++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 = error "jacobian Tower: unimplemented"+jacobian Mixed = M.jacobian+{-# INLINE jacobian #-}++jacobian2 :: (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)+jacobian2 Forward = F.jacobian2+jacobian2 Reverse = R.jacobian2+jacobian2 Tower = error "jacobian2 Tower: unimplemented"+jacobian2 Mixed = M.jacobian2+{-# INLINE jacobian2 #-}++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   = error "grad Tower: unimplemented"+grad Mixed   = M.grad+{-# INLINE grad #-}++grad2 :: (Traversable f, Num a) => Direction -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f a)+grad2 Forward = F.grad2+grad2 Reverse = R.grad2+grad2 Tower   = error "grad2 Tower: unimplemented"+grad2 Mixed   = M.grad2+{-# INLINE grad2 #-}+

+ Numeric/AD/Forward.hs view

@@ -0,0 +1,100 @@+{-# LANGUAGE Rank2Types #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.AD.Forward+-- Copyright   :  (c) Edward Kmett 2010+-- License     :  BSD3+-- Maintainer  :  ekmett@gmail.com+-- Stability   :  experimental+-- Portability :  GHC only +--+-- Forward mode automatic differentiation+--+-----------------------------------------------------------------------------++module Numeric.AD.Forward+    ( +    -- * Gradient+      grad+    , grad2+    -- * Jacobian+    , jacobian+    , jacobian2+    , jacobianT+    -- * Derivatives+    , diffUU+    , diff2UU+    , diffUF+    , diff2UF+    -- * Synonyms+    , diff+    , diff2+    -- * Exposed Types+    , AD(..)+    , Mode(..)+    ) where++import Data.Traversable (Traversable)+import Control.Applicative+import Numeric.AD.Classes+import Numeric.AD.Internal+import Numeric.AD.Internal.Forward++-- | The 'diff2' function calculates the first derivative of scalar-to-scalar function by 'Forward' 'AD'+diff :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> a+diff = diffUU+{-# INLINE diff #-}++-- | The 'diff2' function calculates the result and first derivative of scalar-to-scalar function by 'Forward' 'AD'+diff2 :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a) +diff2 = diff2UU+{-# INLINE diff2 #-}++-- | The 'diffUU' function calculates the first derivative of a scalar-to-scalar function by 'Forward' 'AD'+diffUU :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> a+diffUU f a = tangent $ apply f a+{-# INLINE diffUU #-}++-- | The 'diff2UU' function calculates the result and first derivative of scalar-to-scalar function by 'Forward' 'AD'+diff2UU :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a) +diff2UU f a = unbundle $ apply f a+{-# INLINE diff2UU #-}++-- | The 'diffUF' function calculates the first derivative of scalar-to-nonscalar function by 'Forward' 'AD'+diffUF :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f a+diffUF f a = tangent <$> apply f a+{-# INLINE diffUF #-}++-- | The 'diff2UF' function calculates the result and first derivative of a scalar-to-non-scalar function by 'Forward' 'AD'+diff2UF :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f (a, a) +diff2UF f a = unbundle <$> apply f a+{-# INLINE diff2UF #-}++-- A fast, simple transposed forward jacobian+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)+-- jacobianT f as = fmap tangent <$> bind f as+{-# INLINE jacobianT #-}++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) = bind2 (fmap tangent . f) as +{-# INLINE jacobian #-}++jacobian2 :: (Traversable f, Traversable g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f a)+jacobian2 f as = transposeWith row t p+    where +        (p, t) = bind2 f as +        row x as' = (primal x, tangent <$> as') +{-# INLINE jacobian2 #-}++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 #-}++grad2 :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f a)+grad2 f as = (primal b, tangent <$> bs)+    where +        (b, bs) = bind2 f as+{-# INLINE grad2 #-}

+ Numeric/AD/Internal.hs view

@@ -0,0 +1,128 @@+{-# LANGUAGE GeneralizedNewtypeDeriving, TemplateHaskell, FlexibleContexts, FlexibleInstances #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.AD.Internal+-- Copyright   :  (c) Edward Kmett 2010+-- License     :  BSD3+-- Maintainer  :  ekmett@gmail.com+-- Stability   :  experimental+-- Portability :  GHC only +--+-----------------------------------------------------------------------------+module Numeric.AD.Internal+    ( zipWithT+    , AD(..)+    , Id(..)+    , probe+    , unprobe+    ) where++import Control.Applicative+import Language.Haskell.TH+import Numeric.AD.Classes+import Data.Monoid+import Data.Traversable (Traversable, mapAccumL)+import Data.Foldable (Foldable, toList)++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)++-- | '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 (Lifted, Mode, Primal)++let f = varT (mkName "f") in deriveNumeric (conT ''AD `appT` f) f++newtype Id a = Id a deriving+    (Eq, Ord, Show, Enum, Bounded, Num, Real, Fractional, Floating, RealFrac, RealFloat, Monoid)++probe :: a -> AD Id a+probe a = AD (Id a)++unprobe :: AD Id a -> a+unprobe (AD (Id a)) = a++instance Functor Id where+    fmap 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)++instance Primal Id where+    primal (Id a) = a

+ Numeric/AD/Internal/Forward.hs view

@@ -0,0 +1,112 @@+{-# LANGUAGE Rank2Types, TypeFamilies, DeriveDataTypeable, TemplateHaskell, UndecidableInstances, BangPatterns #-}+-----------------------------------------------------------------------------+-- |+-- 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+    , bind2+    , 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.Classes+import Numeric.AD.Internal++data Forward a = Forward a a deriving (Show, Data, Typeable)++tangent :: AD Forward a -> a+tangent (AD (Forward _ da)) = da+{-# INLINE tangent #-}++unbundle :: AD Forward a -> (a, a)+unbundle (AD (Forward a da)) = (a, da)+{-# 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++instance Lifted Forward => Mode Forward where+    lift a = Forward a 0+    Forward a da <+> Forward b db = Forward (a + b) (da + db)+    a *^ Forward b db = Forward (a * b) (a * db)+    Forward a da ^* b = Forward (a * b) (da * b)+    Forward a da ^/ b = Forward (a / b) (da / b)++instance Lifted Forward => Jacobian Forward where+    type D Forward = Id+    unary f (Id dadb) (Forward b db) = Forward (f b) (dadb * db)+    lift1 f df (Forward b db) = Forward (f b) (dadb * db)+        where +            Id dadb = df (Id b)+    lift1_ f df (Forward b db) = Forward a da+        where +            a = f b+            Id da = df (Id a) (Id b) ^* db++    binary f (Id dadb) (Id dadc) (Forward b db) (Forward c dc) = Forward (f b c) da+        where +            da = dadb * db + dc * dadc+    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 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 $ 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, bundle a $ if i == j then 1 else 0)++bind2 :: (Traversable f, Num a) => (f (AD Forward a) -> b) -> f a -> (b, f b)+bind2 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, bundle a $ if i == j then 1 else 0)+        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+

+ Numeric/AD/Internal/Reverse.hs view

@@ -0,0 +1,202 @@+{-# LANGUAGE Rank2Types, TypeFamilies, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, FlexibleContexts, TemplateHaskell, UndecidableInstances #-}+-----------------------------------------------------------------------------+-- |+-- 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+    , derivative2+    , Var(..)+    ) where++import Prelude hiding (mapM)+import Control.Applicative (Applicative(..),(<$>))+import Control.Monad.ST+import Control.Monad (forM_)+import Data.List (foldl')+import Data.Array.ST+import Data.Array+import Data.IntMap (IntMap, fromListWith, findWithDefault)+import Data.Graph (graphFromEdges', topSort, Vertex)+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 Numeric.AD.Classes+import Numeric.AD.Internal++-- | A @Tape@ records the information needed back propagate from the output to each input during 'Reverse' 'Mode' AD.+data Tape a t+    = Lift a+    | Var a {-# UNPACK #-} !Int+    | Binary a a a t t+    | Unary a a t+    deriving (Show)++-- | @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)++instance MuRef (Reverse a) where+    type DeRef (Reverse a) = Tape a++    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)+    (<+>)  = 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++instance Primal Reverse where+    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 (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 (Lift b)) (Reverse (Lift c)) = Reverse (Lift (f b c))+    binary f _         (Id dadc) (Reverse (Lift b)) c                  = Reverse (Unary (f b (primal c)) dadc c)+    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 (conT ''Reverse)+-- deriveNumeric  ''Reverse++derivative :: Num a => AD Reverse a -> a+derivative = sum . map snd . partials+{-# INLINE derivative #-}++derivative2 :: Num a => AD Reverse a -> (a, a)+derivative2 r = (primal r, derivative r)+{-# INLINE derivative2 #-}++-- | 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++-- | 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_ (topSort g) $+                backPropagate vmap ss+            return ss+        sbounds ((a,_):as) = foldl' (\(lo,hi) (b,_) -> (min lo b, max hi b)) (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 Var t a | t -> a where+    var     :: a -> Int -> t+    fromVar :: t -> Int++    bind :: Traversable f => f a -> (f t, (Int,Int))+    unbind :: Functor f => f t -> Array Int a -> f a+    unbindMap :: (Functor f, Num a) => f t -> IntMap a -> f a++    -- TODO: tweak bounds+    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 xs ys = fmap (\v -> ys ! fromVar v) xs+    unbindMap xs ys = fmap (\v -> findWithDefault 0 (fromVar v) ys) xs++instance Var (Reverse a) a where+    var a v = Reverse (Var a v)+    fromVar (Reverse (Var _ v)) = v+    fromVar _ = error "fromVar: not a Var"++instance Var (AD Reverse a) a where+    var a v = AD (var a v)+    fromVar (AD v) = fromVar v

+ Numeric/AD/Internal/Tower.hs view

@@ -0,0 +1,101 @@+{-# LANGUAGE Rank2Types, TypeFamilies, FlexibleContexts, UndecidableInstances, TemplateHaskell #-}+-----------------------------------------------------------------------------+-- |+-- 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+    , d+    , d2+    , tangents+    , bundle+    , apply+    , getADTower+    ) where++import Numeric.AD.Classes+import Numeric.AD.Internal+import Language.Haskell.TH++-- | @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 (Show)++-- Local combinators++zeroPad :: Num a => [a] -> [a]+zeroPad xs = xs ++ repeat 0 +{-# INLINE zeroPad #-}++d :: Num a => [a] -> a+d (_:da:_) = da+d _ = 0+{-# INLINE d #-}++d2 :: Num a => [a] -> (a, a) +d2 (a:da:_) = (a, da)+d2 (a:_)    = (a, 0)+d2 _        = (0, 0)+{-# INLINE d2 #-}++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 #-}++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 #-}++instance Primal Tower where+    primal (Tower (x:_)) = x+    primal _ = 0++instance Lifted Tower => Mode Tower where+    lift a = Tower [a]+    zero = Tower []++    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 (conT ''Tower)

+ Numeric/AD/Newton.hs view

@@ -0,0 +1,66 @@+{-# LANGUAGE Rank2Types, BangPatterns #-}+-----------------------------------------------------------------------------+-- |+-- 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 Descent (Reverse AD)+    , gradientDescent+    -- * Exposed Types+    , AD(..)+    , Mode(..)+    ) where++import Prelude hiding (all)+import Numeric.AD.Classes+import Numeric.AD.Internal+import Data.Foldable (all)+import Data.Traversable (Traversable)+import Numeric.AD.Forward (diff, diff2) +import Numeric.AD.Reverse (grad2) ++findZero :: Fractional a => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]+findZero f x0 = iterate (\x -> let (y,y') = diff2 f x in x - y/y') x0+{-# INLINE findZero #-}++inverse :: Fractional 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  #-}++fixedPoint :: Fractional a => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]+fixedPoint f = findZero (\x -> f x - x)+{-# INLINE fixedPoint #-}++extremum :: Fractional a => (forall t s. (Mode t, Mode s) => AD t (AD s a) -> AD t (AD s a)) -> a -> [a]+extremum f x0 = findZero (diff f) x0+{-# INLINE extremum #-}++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 gx0 0.1 (0 :: Int)+    where+        (fx0, gx0) = grad2 f x0+        go x fx gx !eta !i +            | eta == 0     = [] -- step size is 0+            | fx1 > fx     = go x fx gx (eta/2) 0 -- we stepped too far+            | all (==0) gx = [] -- gradient is 0+            | otherwise    = x1 : if i == 10+                                  then go x1 fx1 gx1 (eta*2) 0+                                  else go x1 fx1 gx1 eta (i+1)+            where+                -- should check gx = 0 here+                x1 = zipWithT (\xi gxi -> xi - eta * gxi) x gx+                (fx1, gx1) = grad2 f x1+{-# INLINE gradientDescent #-}

+ Numeric/AD/Reverse.hs view

@@ -0,0 +1,115 @@+{-# LANGUAGE Rank2Types, TemplateHaskell, BangPatterns #-}+-----------------------------------------------------------------------------+-- |+-- Module      :  Numeric.AD.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.Reverse+    ( +    -- * Gradient+      grad+    , grad2+    -- * Jacobian+    , jacobian+    , jacobian2+    -- * Derivatives+    , diffUU+    , diff2UU+    , diffFU+    , diff2FU+    , diffUF+    , diff2UF+    -- * Synonyms+    , diff+    , diff2+    -- * Exposed Types+    , AD(..)+    , Mode(..)+    ) where++import Control.Applicative ((<$>))+import Data.Traversable (Traversable)++import Numeric.AD.Classes+import Numeric.AD.Internal+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 'grad2' function calculates the result and gradient of a non-scalar-to-scalar function with 'Reverse' AD in a single pass.+grad2 :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f a)+grad2 f as = (primal r, unbind vs $ partialArray bds r)+    where (vs, bds) = bind as+          r = f vs +{-# INLINE grad2 #-}++-- | 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 'jacobian2' 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'+jacobian2 :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f a)+jacobian2 f as = row <$> f vs where +    (vs, bds) = bind as+    row a = (primal a, unbind vs (partialArray bds a))+{-# INLINE jacobian2 #-}++diffUU :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> a+diffUU f a = derivative $ f (var a 0)+{-# INLINE diffUU #-}++diffUF :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f a+diffUF f a = derivative <$> f (var a 0)+{-# INLINE diffUF #-}++-- | The 'diff2UU' function calculates the value and derivative, as a+-- pair, of a scalar-to-scalar function.+diff2UU :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a)+diff2UU f a = derivative2 $ f (var a 0)+{-# INLINE diff2UU #-}++diff2UF :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f (a, a)+diff2UF f a = derivative2 <$> f (var a 0)+{-# INLINE diff2UF #-}++diffFU :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f a+diffFU f as = unbind vs $ partialArray bds (f vs)+    where (vs, bds) = bind as+{-# INLINE diffFU #-}++diff2FU :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f a)+diff2FU f as = (primal result, unbind vs $ partialArray bds result)+    where (vs, bds) = bind as+          result = f vs+{-# INLINE diff2FU #-}++-- | The 'diff' function is a synonym for 'diffUU'.+diff :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> a+diff = diffUU +{-# INLINE diff #-}++-- | The 'diff2' function is a synonym for 'diff2UU'.+diff2 :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a)+diff2 = diff2UU+{-# INLINE diff2 #-}+

+ Numeric/AD/Tower.hs view

@@ -0,0 +1,98 @@+{-# LANGUAGE Rank2Types, BangPatterns #-}+-----------------------------------------------------------------------------+-- |+-- Module      : Numeric.AD.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.Tower+    ( +    -- * Taylor Series+      taylor, taylor0+    , maclaurin, maclaurin0+    -- * Derivatives+    , diffUU+    , diff2UU+    , diffsUU+    , diffs0UU+    , diffsUF+    , diffs0UF +    -- * Synonyms+    , diffs, diffs0+    , diff, diff2+    -- * Exposed Types+    , Mode(..)+    , AD(..)+    ) where++-- TODO: argminNaiveGradient++import Control.Applicative ((<$>))+import Numeric.AD.Classes+import Numeric.AD.Internal+import Numeric.AD.Internal.Tower++diffsUU :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]+diffsUU f a = getADTower $ apply f a +{-# INLINE diffsUU #-}++diffs0UU :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]+diffs0UU f a = zeroPad (diffsUU f a)+{-# INLINE diffs0UU #-}++diffs0UF :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f [a]+diffs0UF f a = (zeroPad . getADTower) <$> apply f a+{-# INLINE diffs0UF #-}++diffsUF :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f [a]+diffsUF f a = getADTower <$> apply f a+{-# INLINE diffsUF #-}++diffs :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]+diffs = diffsUU+{-# INLINE diffs #-}++diffs0 :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]+diffs0 = diffs0UU+{-# INLINE diffs0 #-}++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 #-}++diffUU :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> a+diffUU f a = d $ diffs f a +{-# INLINE diffUU #-}++diff2UU :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a)+diff2UU f a = d2 $ diffs f a +{-# INLINE diff2UU #-}++diff :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> a+diff = diffUU+{-# INLINE diff #-}++diff2 :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a)+diff2 = diff2UU+{-# INLINE diff2 #-}

+ Setup.hs view

@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain

+ TODO view

@@ -0,0 +1,13 @@+* Implement the diffMF, etc. functionality from Numeric.FAD++* Allow the type to vary within our AD data type container, in the same fashion as Numeric.FAD.++    Although, while Pearlmutter and Siskind provided the functionality to permit it in derivative combinator, they provided+    no combinators to convert, for instance, a Dual tag Float to a Dual tag Double, so that extra functionality cannot+    currently be leveraged.++* Provide the ability to use reverse mode locally on FAD inputs, i.e.++    reverseCheckpoint :: (forall s. AD s a -> AD s a) -> AD t a -> AD t a ++* Provide forward-on-reverse computation of Hessians

+ ad.cabal view

@@ -0,0 +1,37 @@+Name:         ad+Version:      0.12+License:      BSD3+License-File: LICENSE+Copyright:    Edward Kmett 2010+Author:       Edward Kmett 2010+Maintainer:   ekmett@gmail.com+Stability:    Experimental+Category:     Math+Homepage:     http://comonad.com/reader/+Synopsis:     Mixed-Mode Automatic Differentiation.+Description:  Forward, reverse, and higher-order automatic differentiation with a common API++Build-Type:   Simple+Build-Depends:       +    base >= 4 && < 5,+    data-reify >= 0.5 && < 0.6, +    containers >= 0.2 && < 0.4,+    template-haskell >= 2.4 && < 2.5,+    array >= 0.2 && < 0.4++Exposed-Modules:+    Numeric.AD+    Numeric.AD.Forward+    Numeric.AD.Reverse+    Numeric.AD.Tower+    Numeric.AD.Directed+    Numeric.AD.Newton+    Numeric.AD.Classes+    Numeric.AD.Internal+    Numeric.AD.Internal.Forward+    Numeric.AD.Internal.Reverse+    Numeric.AD.Internal.Tower++Extra-Source-Files: TODO++GHC-Options: -Wall