packages feed

ad 3.4 → 4.0

raw patch · 39 files changed

+2310/−2407 lines, 39 filesdep +addep +criteriondep +transformersdep ~arraydep ~basedep ~comonad

Dependencies added: ad, criterion, transformers

Dependency ranges changed: array, base, comonad, erf, free, mtl, reflection, tagged, template-haskell

Files

.ghci view
@@ -1,1 +1,1 @@-:set -isrc -idist/build/autogen -optP-include -optPdist/build/autogen/cabal_macros.h+:set -isrc -idist/build/autogen -optP-include -optPdist/build/autogen/cabal_macros.h -optP-include -optPinclude
.travis.yml view
@@ -12,7 +12,7 @@  script:   - $script-  - hlint src --cpp-define HLINT+  - hlint src --cpp-define HLINT --cpp-include include  notifications:   irc:@@ -23,4 +23,4 @@       - "\x0313ad\x03/\x0306%{branch}\x03 \x0314%{commit}\x03 %{build_url} %{message}"  env:-  - mode="--enable-tests" script="cabal test"+  - mode="--enable-tests" script="cabal test --show-details=always"
CHANGELOG.markdown view
@@ -1,6 +1,12 @@+4.0+---+* An overhaul permitting monomorphic modes was completed by @alang9.+* Add a `ForwardDouble` monomorphic mode+ 3.4 --- * Added support for `erf` and `inverf`, etc. from `Data.Number.Erf`.+* Split the infinitesimal and mode into two separate parameters to facilitate inlining and easier extension of the API.  3.3.1 -----
README.markdown view
@@ -23,7 +23,7 @@  You can compute derivatives of functions -    Prelude Numeric.AD> diff sin 0 {-# cos 0 #-}+    Prelude Numeric.AD> diff sin 0 {- cos 0 -}     1.0  Or both the answer and the derivative of a function:@@ -59,30 +59,32 @@     Prelude Numeric.AD> take 10 $ diffs sin 1     [0.8414709848078965,0.5403023058681398,-0.8414709848078965,-0.5403023058681398,0.8414709848078965,0.5403023058681398,-0.8414709848078965,-0.5403023058681398,0.8414709848078965,0.5403023058681398] -or if your function takes multiple inputs, you can use grads, which returns an 'f-branching stream' of derivatives. Somewhat more intuitive answers can be obtained by converting the stream into the-polymorphically recursive `Tensors` data type. With that we can look at a single 'layer' of the answer at a time:+or if your function takes multiple inputs, you can use grads, which returns an 'f-branching stream' of derivatives, that you can+inspect lazily. Somewhat more intuitive answers can be obtained by converting the stream into the polymorphically recursive +`Jet` data type. With that we can look at a single "layer" of the answer at a time:  The answer: -    Prelude Numeric.AD> headJet $ tensors $  grads (\[x,y] -> exp (x * y)) [1,2]+    Prelude Numeric.AD Numeric.AD.Types> headJet $ jet $  grads (\[x,y] -> exp (x * y)) [1,2]     7.38905609893065  The gradient: -    Prelude Numeric.AD> headJet $ tailJet $ tensors $  grads (\[x,y] -> exp (x * y)) [1,2]+    Prelude Numeric.AD Numeric.AD.Types> headJet $ tailJet $ jet $  grads (\[x,y] -> exp (x * y)) [1,2]     [14.7781121978613,7.38905609893065]  The hessian (n * n matrix of 2nd derivatives) -    Prelude Numeric.AD> headJet $ tailJet $ tailJet $ tensors $  grads (\[x,y] -> exp (x * y)) [1,2]+    Prelude Numeric.AD Numeric.AD.Types> headJet $ tailJet $ tailJet $ jet $  grads (\[x,y] -> exp (x * y)) [1,2]     [[29.5562243957226,22.16716829679195],[22.16716829679195,7.38905609893065]]  Or even higher order tensors of derivatives. -    Prelude Numeric.AD> headJet $ tailJet $ tailJet $ tailJet $ tensors $  grads (\[x,y] -> exp (x * y)) [1,2]+    Prelude Numeric.AD Numeric.AD.Types> headJet $ tailJet $ tailJet $ tailJet $ jet $  grads (\[x,y] -> exp (x * y)) [1,2]     [[[59.1124487914452,44.3343365935839],[44.3343365935839,14.7781121978613]],[[44.3343365935839,14.7781121978613],[14.7781121978613,7.38905609893065]]] -Note the redundant values caused by the various symmetries in the tensors. The 'ad' library is careful to compute each distinct derivative only once and to share the resulting thunks.+Note the redundant values caused by the various symmetries in the tensors. The `ad` library is careful to compute +each distinct derivative only once, lazily and to share the resulting computation.  Overview --------
ad.cabal view
@@ -1,8 +1,8 @@ name:         ad-version:      3.4+version:      4.0 license:      BSD3 license-File: LICENSE-copyright:    (c) Edward Kmett 2010-2013,+copyright:    (c) Edward Kmett 2010-2014,               (c) Barak Pearlmutter and Jeffrey Mark Siskind 2008-2009 author:       Edward Kmett maintainer:   ekmett@gmail.com@@ -81,6 +81,7 @@ library   extensions: CPP   hs-source-dirs: src+  include-dirs: include    other-extensions:     BangPatterns@@ -99,51 +100,48 @@     UndecidableInstances    build-depends:-    array            >= 0.2 && < 0.5,-    base             == 4.*,-    comonad          >= 3,-    containers       >= 0.2 && < 0.6,-    data-reify       >= 0.6 && < 0.7,-    erf              >= 2.0 && < 2.1,-    free             >= 3,-    mtl              >= 2,-    reflection       >= 1.1.6,-    tagged           >= 0.4.2.1,-    template-haskell >= 2.5 && < 2.9+    array            >= 0.2   && < 0.6,+    base             >= 4.5   && < 5,+    comonad          >= 4     && < 5,+    containers       >= 0.2   && < 0.6,+    data-reify       >= 0.6   && < 0.7,+    erf              >= 2.0   && < 2.1,+    free             >= 4.6.1 && < 5,+    mtl              >= 2     && < 2.2,+    reflection       >= 1.4   && < 2,+    tagged           >= 0.7   && < 1,+    template-haskell,+    transformers     >= 0.3   && < 0.4    exposed-modules:     Numeric.AD-    Numeric.AD.Types -    Numeric.AD.Newton     Numeric.AD.Halley+    Numeric.AD.Jacobian+    Numeric.AD.Jet+    Numeric.AD.Newton +    Numeric.AD.Mode     Numeric.AD.Mode.Directed     Numeric.AD.Mode.Forward+    Numeric.AD.Mode.Forward.Double     Numeric.AD.Mode.Kahn     Numeric.AD.Mode.Reverse     Numeric.AD.Mode.Tower     Numeric.AD.Mode.Sparse -    Numeric.AD.Variadic-    Numeric.AD.Variadic.Kahn-    Numeric.AD.Variadic.Sparse--    Numeric.AD.Internal.Classes-    Numeric.AD.Internal.Combinators+    Numeric.AD.Internal.Dense     Numeric.AD.Internal.Forward-    Numeric.AD.Internal.Tower+    Numeric.AD.Internal.Forward.Double+    Numeric.AD.Internal.Identity     Numeric.AD.Internal.Kahn+    Numeric.AD.Internal.On     Numeric.AD.Internal.Reverse-    Numeric.AD.Internal.Var     Numeric.AD.Internal.Sparse-    Numeric.AD.Internal.Dense-    Numeric.AD.Internal.Composition+    Numeric.AD.Internal.Tower    other-modules:-    Numeric.AD.Internal.Types-    Numeric.AD.Internal.Jet-    Numeric.AD.Internal.Identity+    Numeric.AD.Internal.Combinators    if flag(lib-Werror)     ghc-options: -Werror@@ -166,3 +164,10 @@   if impl(ghc<7.6)     ghc-options: -Werror   hs-source-dirs: tests++benchmark blackscholes+  type: exitcode-stdio-1.0+  main-is: BlackScholes.hs+  hs-source-dirs: bench+  build-depends: base, ad, erf, criterion+  ghc-options: -fspec-constr -fdicts-cheap -O2
+ bench/BlackScholes.hs view
@@ -0,0 +1,74 @@+{-# LANGUAGE RankNTypes #-}+import Criterion.Main+import Data.Number.Erf+import qualified Numeric.AD as Mixed+import qualified Numeric.AD.Mode.Forward as Forward+import qualified Numeric.AD.Mode.Kahn as Kahn+import qualified Numeric.AD.Mode.Reverse as Reverse+import qualified Numeric.AD.Mode.Sparse as Sparse++blackScholes :: (Erf a) => a -> a -> a -> a -> a -> (a, a)+blackScholes r s v t k = (put, call)+  where+    put = k * exp (negate r * t) - s + call+    call = normcdf (negate d2) * k * exp (negate r * t) - normcdf (negate d1) * s+    d1 = (log (s / k) + (r + v * v / 2) * t) / (v * sqrt t)+    d2 = d1 - v * t++bs :: Erf a => [a] -> (a, a)+bs [r', s', v', t', k'] = blackScholes r' s' v' t' k'++fromPair :: (t, t) -> [t]+fromPair (a, b) = [a, b]++runF :: Num a => (a -> a -> a -> a -> a -> b) -> Int -> [b]+runF f n =+    [ f r s v t k+    | r <- xs, s <- xs, v <- xs, t <- xs, k <- xs]+  where+    xs = map fromIntegral [1..n]++runFloat :: (Float -> Float -> Float -> Float -> Float -> b) -> Int -> [b]+runFloat = runF++runDouble :: (Double -> Double -> Double -> Double -> Double -> b) -> Int -> [b]+runDouble = runF++main = defaultMain+    [ bgroup "Forward"+        [ bench "greeks Double" $ nf (runDouble $ \r s v t k -> Forward.jacobian (fromPair . bs) [r, s, v, t, k]) 2+        , bench "greeks Float" $ nf (runFloat $ \r s v t k -> Forward.jacobian (fromPair . bs) [r, s, v, t, k]) 2+        ]+    , bgroup "Kahn"+        [ bench "greeks Double" $ nf (runDouble $ \r s v t k -> Kahn.jacobian (fromPair . bs) [r, s, v, t, k]) 2+        , bench "higherGreeks Double" $ nf (runDouble $ \r s v t k -> (Kahn.hessian (fst . bs) [r, s, v, t, k], Kahn.hessian (snd . bs) [r, s, v, t, k])) 2+        , bench "highererGreeks Double" $ nf (runDouble $ \r s v t k -> Kahn.hessianF (fromPair . bs) [r, s, v, t, k]) 2+        , bench "greeks Float" $ nf (runFloat $ \r s v t k -> Kahn.jacobian (fromPair . bs) [r, s, v, t, k]) 2+        , bench "higherGreeks Float" $ nf (runFloat $ \r s v t k -> (Kahn.hessian (fst . bs) [r, s, v, t, k], Kahn.hessian (snd . bs) [r, s, v, t, k])) 2+        , bench "highererGreeks Float" $ nf (runFloat $ \r s v t k -> Kahn.hessianF (fromPair . bs) [r, s, v, t, k]) 2+        ]+    , bgroup "Reverse"+        [ bench "greeks Double" $ nf (runDouble $ \r s v t k -> Reverse.jacobian (fromPair . bs) [r, s, v, t, k]) 2+        , bench "higherGreeks Double" $ nf (runDouble $ \r s v t k -> (Reverse.hessian (fst . bs) [r, s, v, t, k], Reverse.hessian (snd . bs) [r, s, v, t, k])) 2+        , bench "highererGreeks Double" $ nf (runDouble $ \r s v t k -> Reverse.hessianF (fromPair . bs) [r, s, v, t, k]) 2+        , bench "greeks Float" $ nf (runFloat $ \r s v t k -> Reverse.jacobian (fromPair . bs) [r, s, v, t, k]) 2+        , bench "higherGreeks Float" $ nf (runFloat $ \r s v t k -> (Reverse.hessian (fst . bs) [r, s, v, t, k], Reverse.hessian (snd . bs) [r, s, v, t, k])) 2+        , bench "highererGreeks Float" $ nf (runFloat $ \r s v t k -> Reverse.hessianF (fromPair . bs) [r, s, v, t, k]) 2+        ]+    , bgroup "Sparse"+        [ bench "greeks Double" $ nf (runDouble $ \r s v t k -> Sparse.jacobian (fromPair . bs) [r, s, v, t, k]) 2+        , bench "higherGreeks Double" $ nf (runDouble $ \r s v t k -> (Sparse.hessian (fst . bs) [r, s, v, t, k], Sparse.hessian (snd . bs) [r, s, v, t, k])) 2+        , bench "highererGreeks Double" $ nf (runDouble $ \r s v t k -> Sparse.hessianF (fromPair . bs) [r, s, v, t, k]) 2+        , bench "greeks Float" $ nf (runFloat $ \r s v t k -> Sparse.jacobian (fromPair . bs) [r, s, v, t, k]) 2+        , bench "higherGreeks Float" $ nf (runFloat $ \r s v t k -> (Sparse.hessian (fst . bs) [r, s, v, t, k], Sparse.hessian (snd . bs) [r, s, v, t, k])) 2+        , bench "highererGreeks Float" $ nf (runFloat $ \r s v t k -> Sparse.hessianF (fromPair . bs) [r, s, v, t, k]) 2+        ]+--    , bgroup "Mixed"+--        [ bench "greeks Double" $ nf (runDouble $ \r s v t k -> Mixed.jacobian (fromPair . bs) [r, s, v, t, k]) 2+--        , bench "higherGreeks Double" $ nf (runDouble $ \r s v t k -> (Mixed.hessian (fst . bs) [r, s, v, t, k], Mixed.hessian (snd . bs) [r, s, v, t, k])) 2+--        , bench "highererGreeks Double" $ nf (runDouble $ \r s v t k -> Mixed.hessianF (fromPair . bs) [r, s, v, t, k]) 2+--        , bench "greeks Float" $ nf (runFloat $ \r s v t k -> Mixed.jacobian (fromPair . bs) [r, s, v, t, k]) 2+--        , bench "higherGreeks Float" $ nf (runFloat $ \r s v t k -> (Mixed.hessian (fst . bs) [r, s, v, t, k], Mixed.hessian (snd . bs) [r, s, v, t, k])) 2+--        , bench "highererGreeks Float" $ nf (runFloat $ \r s v t k -> Mixed.hessianF (fromPair . bs) [r, s, v, t, k]) 2+--        ]+    ]
src/Numeric/AD.hs view
@@ -1,8 +1,10 @@-{-# LANGUAGE Rank2Types, TypeFamilies, PatternGuards #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE PatternGuards #-} ----------------------------------------------------------------------------- -- |--- Module      :  Numeric.AD--- Copyright   :  (c) Edward Kmett 2010+-- Copyright   :  (c) Edward Kmett 2010-2014 -- License     :  BSD3 -- Maintainer  :  ekmett@gmail.com -- Stability   :  experimental@@ -39,160 +41,161 @@ -----------------------------------------------------------------------------  module Numeric.AD-    (-    -- * Gradients (Reverse Mode)-      grad-    , grad'-    , gradWith-    , gradWith'+  ( -    -- * Higher Order Gradients (Sparse-on-Reverse)-    , grads+  -- * AD modes+    Mode(auto)+  , Scalar -    -- * Jacobians (Sparse or Reverse)-    , jacobian-    , jacobian'-    , jacobianWith-    , jacobianWith'+  -- * Gradients (Reverse Mode)+  , grad+  , grad'+  , gradWith+  , gradWith' -    -- * Higher Order Jacobian (Sparse-on-Reverse)-    , jacobians+  -- * Higher Order Gradients (Sparse-on-Reverse)+  , grads -    -- * Transposed Jacobians (Forward Mode)-    , jacobianT-    , jacobianWithT+  -- * Variadic Gradients (Sparse or Kahn)+  -- $vgrad+  , Grad , vgrad, vgrad'+  , Grads, vgrads -    -- * Hessian (Sparse-On-Reverse)-    , hessian-    , hessian'+  -- * Jacobians (Sparse or Reverse)+  , jacobian+  , jacobian'+  , jacobianWith+  , jacobianWith' -    -- * Hessian Tensors (Sparse or Sparse-On-Reverse)-    , hessianF-    -- * Hessian Tensors (Sparse)-    , hessianF'+  -- * Higher Order Jacobian (Sparse-on-Reverse)+  , jacobians -    -- * Hessian Vector Products (Forward-On-Reverse)-    , hessianProduct-    , hessianProduct'+  -- * Transposed Jacobians (Forward Mode)+  , jacobianT+  , jacobianWithT -    -- * Derivatives (Forward Mode)-    , diff-    , diffF+  -- * Hessian (Sparse-On-Reverse)+  , hessian+  , hessian' -    , diff'-    , diffF'+  -- * Hessian Tensors (Sparse or Sparse-On-Reverse)+  , hessianF -    -- * Derivatives (Tower)-    , diffs-    , diffsF+  -- * Hessian Tensors (Sparse)+  , hessianF' -    , diffs0-    , diffs0F+  -- * Hessian Vector Products (Forward-On-Reverse)+  , hessianProduct+  , hessianProduct' -    -- * Directional Derivatives (Forward Mode)-    , du-    , du'-    , duF-    , duF'+  -- * Derivatives (Forward Mode)+  , diff+  , diffF -    -- * Directional Derivatives (Tower)-    , dus-    , dus0-    , dusF-    , dus0F+  , diff'+  , diffF' -    -- * Taylor Series (Tower)-    , taylor-    , taylor0+  -- * Derivatives (Tower)+  , diffs+  , diffsF -    -- * Maclaurin Series (Tower)-    , maclaurin-    , maclaurin0+  , diffs0+  , diffs0F -    -- * Gradient Descent-    , gradientDescent-    , gradientAscent-    , conjugateGradientDescent-    , conjugateGradientAscent-    ) where+  -- * Directional Derivatives (Forward Mode)+  , du+  , du'+  , duF+  , duF' -import Data.Traversable (Traversable)-import Data.Foldable (Foldable, foldr')+  -- * Directional Derivatives (Tower)+  , dus+  , dus0+  , dusF+  , dus0F++  -- * Taylor Series (Tower)+  , taylor+  , taylor0++  -- * Maclaurin Series (Tower)+  , maclaurin+  , maclaurin0++  -- * Gradient Descent+  , gradientDescent+  , gradientAscent+  , conjugateGradientDescent+  , conjugateGradientAscent++  ) where+ import Control.Applicative+import Data.Functor.Compose+import Data.Traversable (Traversable)+import Data.Reflection (Reifies)+import Numeric.AD.Internal.Forward (Forward)+import Numeric.AD.Internal.Kahn (Grad, vgrad, vgrad')+import Numeric.AD.Internal.On+import Numeric.AD.Internal.Reverse (Reverse, Tape)+import Numeric.AD.Internal.Sparse (Sparse, Grads, vgrads) -import Numeric.AD.Types-import Numeric.AD.Internal.Composition-import Numeric.AD.Internal.Identity+import Numeric.AD.Mode  import Numeric.AD.Mode.Forward-    ( diff, diff', diffF, diffF'-    , du, du', duF, duF'-    , jacobianT, jacobianWithT )+  ( 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 )+  ( 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')+  ( 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')+  ( grads, jacobians, hessian', hessianF')  import Numeric.AD.Newton --- | 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.+-- | Calculate the Jacobian of a non-scalar-to-non-scalar function, automatically choosing between sparse and Reverse mode AD. ----- 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)+-- If you know that you have relatively many outputs per input, consider using 'Numeric.AD.Sparse.jacobian'.+--+-- >>> jacobian (\[x,y] -> [y,x,x+y,x*y,exp x * sin y]) [pi,1]+-- [[0.0,1.0],[1.0,0.0],[1.0,1.0],[1.0,3.141592653589793],[19.472221418841606,12.502969588876512]]+jacobian :: (Traversable f, Functor g, Num a) => (forall s. Reifies s Tape => f (Reverse a s) -> g (Reverse a s)) -> 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+-- | Calculate both the answer and Jacobian of a non-scalar-to-non-scalar function, using reverse-mode AD. ----- 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+-- If you have relatively many outputs per input, consider using 'Numeric.AD.Sparse.jacobian''.+jacobian' :: (Traversable f, Functor g, Num a) => (forall s. Reifies s Tape => f (Reverse a s) -> g (Reverse a s)) -> f a -> g (a, f a)+jacobian' = Reverse.jacobian' {-# 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.+-- | @'jacobianWith' g f@ calculates the Jacobian of a non-scalar-to-non-scalar function, using Reverse mode AD. -- -- 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)+-- If you know that you have relatively many outputs per input, consider using 'Numeric.AD.Sparse.jacobianWith'.+jacobianWith :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. Reifies s Tape => f (Reverse a s) -> g (Reverse a s)) -> 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.+-- | @'jacobianWith'' g f@ calculates the answer and Jacobian of a non-scalar-to-non-scalar function, using Reverse mode AD. -- -- 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+-- If you know that you have relatively many outputs per input, consider using 'Numeric.AD.Sparse.jacobianWith''.+jacobianWith' :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. Reifies s Tape => f (Reverse a s) -> g (Reverse a s)) -> f a -> g (a, f b)+jacobianWith' = Reverse.jacobianWith' {-# INLINE jacobianWith' #-}  -- | @'hessianProduct' f wv@ computes the product of the hessian @H@ of a non-scalar-to-scalar function @f@ at @w = 'fst' <$> wv@ with a vector @v = snd <$> wv@ using \"Pearlmutter\'s method\" from <http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.29.6143>, which states:@@ -201,23 +204,36 @@ -- -- 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 :: (Traversable f, Num a) => (forall s s'. Reifies s Tape => f (On (Reverse (Forward a s') s)) -> On (Reverse (Forward a s') s)) -> f (a, a) -> f a+hessianProduct f = duF (grad (off . f . fmap On))  -- | @'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))+hessianProduct' :: (Traversable f, Num a) => (forall s s'. Reifies s Tape => f (On (Reverse (Forward a s') s)) -> On (Reverse (Forward a s') s)) -> f (a, a) -> f (a, a)+hessianProduct' f = duF' (grad (off . f . fmap On))  -- | 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))+--+-- >>> hessian (\[x,y] -> x*y) [1,2]+-- [[0,1],[1,0]]+hessian :: (Traversable f, Num a) => (forall s s'. Reifies s Tape => f (On (Reverse (Sparse a s') s)) -> On (Reverse (Sparse a s') s)) -> f a -> f (f a)+hessian f = Sparse.jacobian (grad (off . f . fmap On)) --- | 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+-- | Compute the order 3 Hessian tensor on a non-scalar-to-non-scalar function using 'Sparse'-on-'Reverse'+--+-- >>> hessianF (\[x,y] -> [x*y,x+y,exp x*cos y]) [1,2]+-- [[[0.0,1.0],[1.0,0.0]],[[0.0,0.0],[0.0,0.0]],[[-1.1312043837568135,-2.4717266720048188],[-2.4717266720048188,1.1312043837568135]]]+hessianF :: (Traversable f, Functor g, Num a) => (forall s s'. Reifies s Tape => f (On (Reverse (Sparse a s') s)) -> g (On (Reverse (Sparse a s') s))) -> f a -> g (f (f a))+hessianF f as = getCompose $ Sparse.jacobian (Compose . Reverse.jacobian (fmap off . f . fmap On)) as++-- $vgrad+--+-- 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 'auto' calls you use when taking the gradient+-- of a function that takes gradients!
src/Numeric/AD/Halley.hs view
@@ -1,8 +1,8 @@-{-# LANGUAGE Rank2Types, ScopedTypeVariables #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE ScopedTypeVariables #-} ----------------------------------------------------------------------------- -- |--- Module      :  Numeric.AD.Halley--- Copyright   :  (c) Edward Kmett 2010+-- Copyright   :  (c) Edward Kmett 2010-2014 -- License     :  BSD3 -- Maintainer  :  ekmett@gmail.com -- Stability   :  experimental@@ -16,20 +16,25 @@ -----------------------------------------------------------------------------  module Numeric.AD.Halley-    (-    -- * Halley's Method (Tower AD)-      findZero-    , inverse-    , fixedPoint-    , extremum-    ) where+  (+  -- * Halley's Method (Tower AD)+    findZero+  , inverse+  , fixedPoint+  , extremum+  ) where  import Prelude hiding (all)-import Numeric.AD.Types+import Numeric.AD.Internal.Forward (Forward)+import Numeric.AD.Internal.On+import Numeric.AD.Internal.Tower (Tower)+import Numeric.AD.Mode import Numeric.AD.Mode.Tower (diffs0) import Numeric.AD.Mode.Forward (diff) -- , diff')-import Numeric.AD.Internal.Composition +-- $setup+-- >>> import Data.Complex+ -- | 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.) If the stream becomes constant@@ -40,10 +45,9 @@ -- >>> take 10 $ findZero (\x->x^2-4) 1 -- [1.0,1.8571428571428572,1.9997967892704736,1.9999999999994755,2.0] ----- >>> import Data.Complex -- >>> last $ take 10 $ findZero ((+1).(^2)) (1 :+ 1) -- 0.0 :+ 1.0-findZero :: (Fractional a, Eq a) => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]+findZero :: (Fractional a, Eq a) => (forall s. Tower a s -> Tower a s) -> a -> [a] findZero f = go where   go x = x : if x == xn then [] else go xn where     (y:y':y'':_) = diffs0 f x@@ -57,7 +61,7 @@ -- -- Note: the @take 10 $ inverse sqrt 1 (sqrt 10)@ example that works for Newton's method -- fails with Halley's method because the preconditions do not hold!-inverse :: (Fractional a, Eq a) => (forall s. Mode s => AD s a -> AD s a) -> a -> a -> [a]+inverse :: (Fractional a, Eq a) => (forall s. Tower a s -> Tower a s) -> a -> a -> [a] inverse f x0 y = findZero (\x -> f x - auto y) x0 {-# INLINE inverse  #-} @@ -70,10 +74,11 @@ -- -- >>> last $ take 10 $ fixedPoint cos 1 -- 0.7390851332151607-fixedPoint :: (Fractional a, Eq a) => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]+fixedPoint :: (Fractional a, Eq a) => (forall s. Tower a s -> Tower a s) -> 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.) If the stream becomes@@ -81,6 +86,6 @@ -- -- >>> take 10 $ extremum cos 1 -- [1.0,0.29616942658570555,4.59979519460002e-3,1.6220740159042513e-8,0.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))+extremum :: (Fractional a, Eq a) => (forall s s'. On (Forward (Tower a s') s) -> On (Forward (Tower a s') s)) -> a -> [a]+extremum f = findZero (diff (off . f . On)) {-# INLINE extremum #-}
− src/Numeric/AD/Internal/Classes.hs 
@@ -1,343 +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 Data.Function (on)-import Data.Number.Erf-import Language.Haskell.TH--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-    erf1            :: Erf a => t a -> t a-    erfc1           :: Erf a => t a -> t a-    normcdf1        :: Erf a => t a -> t a-    inverf1         :: InvErf a => t a -> t a-    inverfc1        :: InvErf a => t a -> t a-    invnormcdf1     :: InvErf a => t 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-    auto  :: 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 = auto a *! b-    a ^* b = a *! auto b--    a ^/ b = a ^* recip b--    zero = auto 0--one :: (Mode t, Num a) => t a-one = auto 1-{-# INLINE one #-}--negOne :: (Mode t, Num a) => t a-negOne = auto (-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 (auto . primal) would discard derivative--- information. The end user is protected from accidentally using this function--- by the universal quantification on the various combinators we expose.--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     = auto maxBound-        minBound1     = auto minBound-        showsPrec1 d  = showsPrec d . primal-        fromInteger1 0 = zero-        fromInteger1 n = auto (fromInteger n)-        (+!)          = (<+>) -- binary (+) one one-        (-!)          = binary (-) one negOne -- TODO: <-> ? as it is, this might be pretty bad for Tower-        (*!)          = lift2 (*) (\x y -> (y, x))-        negate1       = lift1 negate (const negOne)-        abs1          = lift1 abs signum1-        signum1       = lift1 signum (const zero)-        fromRational1 0 = zero-        fromRational1 r = auto (fromRational r)-        x /! y        = x *! recip1 y-        recip1        = lift1_ recip (const . negate1 . square1)-        pi1       = auto pi-        exp1      = lift1_ exp const-        log1      = lift1 log recip1-        logBase1 x y = log1 y /! log1 x-        sqrt1     = lift1_ sqrt (\z _ -> recip1 (auto 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      = lift1 tan $ recip1 . square1 . cos1-        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     = lift1 tanh $ recip1 . square1 . cosh1-        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               = auto . 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 = auto (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--        erf1 = lift1 erf $ \x -> (fromInteger1 2 /! sqrt1 pi1) *! exp1 (negate1 x *! x)-        erfc1 = lift1 erfc $ \x -> (fromInteger1 (-2) /! sqrt1 pi1) *! exp1 (negate1 x *! x)-        normcdf1 = lift1 normcdf $ \x -> (fromInteger1 (-1) /! sqrt1 pi1) *! exp1 (x *! x *! fromRational1 (- recip 2) /! sqrt1 (fromInteger1 2))--        inverf1 = lift1 inverfc $ \x -> recip1 $ (fromInteger1 2 /! sqrt1 pi1) *! exp1 (negate1 x *! x)-        inverfc1 = lift1 inverfc $ \x -> recip1 $ negate1 (fromInteger1 2 /! sqrt1 pi1) *! exp1 (negate1 x *! x)-        invnormcdf1 = lift1 invnormcdf $ \x -> recip1 $ (fromInteger1 (-1) /! sqrt1 pi1) *! exp1 (x *! x *! fromRational1 (- recip 2) /! sqrt1 (fromInteger1 2)) |]--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, ''Erf, ''InvErf]-    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'
src/Numeric/AD/Internal/Combinators.hs view
@@ -1,8 +1,8 @@ {-# OPTIONS_GHC -fno-warn-incomplete-patterns #-}+{-# OPTIONS_HADDOCK not-home #-} ----------------------------------------------------------------------------- -- |--- Module      :  Numeric.AD.Internal.Combinators--- Copyright   :  (c) Edward Kmett 2010+-- Copyright   :  (c) Edward Kmett 2010-2014 -- License     :  BSD3 -- Maintainer  :  ekmett@gmail.com -- Stability   :  experimental@@ -11,12 +11,16 @@ -- Combinators used internally by @Numeric.AD@ ----------------------------------------------------------------------------- module Numeric.AD.Internal.Combinators-    ( zipWithT-    , zipWithDefaultT-    ) where+  ( zipWithT+  , zipWithDefaultT+  , withPrimal+  , fromBy+  ) where  import Data.Traversable (Traversable, mapAccumL) import Data.Foldable (Foldable, toList)+import Numeric.AD.Mode+import Numeric.AD.Jacobian  -- | Zip a @'Foldable' f@ with a @'Traversable' g@ assuming @f@ has at least as many entries as @g@. zipWithT :: (Foldable f, Traversable g) => (a -> b -> c) -> f a -> g b -> g c@@ -25,3 +29,12 @@ -- | Zip a @'Foldable' f@ with a @'Traversable' g@ assuming @f@, using a default value after @f@ is exhausted. 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)++-- | Used internally to define various 'Enum' combinators.+withPrimal :: Jacobian t => t -> Scalar t -> t+withPrimal t a = unary (const a) 1 t+{-# INLINE withPrimal #-}++-- | Used internally to define various 'Enum' combinators.+fromBy :: Jacobian t => t -> t -> Int -> Scalar t -> t+fromBy a delta n x = binary (\_ _ -> x) 1 (fromIntegral n) a delta
− src/Numeric/AD/Internal/Composition.hs 
@@ -1,208 +0,0 @@-{-# LANGUAGE CPP, Rank2Types, TypeFamilies, MultiParamTypeClasses, FlexibleInstances, FlexibleContexts, UndecidableInstances, TypeOperators #-}--------------------------------------------------------------------------------- |--- 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--#ifndef MIN_VERSION_base-#define MIN_VERSION_base(x,y,z) 1-#endif--import Control.Applicative hiding ((<**>))-import Data.Data (Data(..), mkDataType, DataType, mkConstr, Constr, constrIndex, Fixity(..))-#if MIN_VERSION_base(4,4,0)-import Data.Typeable (Typeable1(..), Typeable(..), TyCon, mkTyCon3, mkTyConApp, typeOfDefault, gcast1)-#else-import Data.Typeable (Typeable1(..), Typeable(..), TyCon, mkTyCon, mkTyConApp, typeOfDefault, gcast1)-#endif-import Data.Foldable (Foldable(foldMap))-import Data.Traversable (Traversable(traverse))-import Numeric.AD.Internal.Classes-import Numeric.AD.Internal.Types--{-# ANN module "Hlint: ignore Eta reduce" #-}-{-# ANN module "Hlint: ignore Reduce duplication" #-}---- | 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-#if MIN_VERSION_base(4,4,0)-composeFunctorTyCon = mkTyCon3 "ad" "Numeric.AD.Internal.Composition" "ComposeFunctor"-#else-composeFunctorTyCon = mkTyCon "Numeric.AD.Internal.Composition.ComposeFunctor"-#endif--{-# 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-    auto = ComposeMode . auto . auto-    ComposeMode a <+> ComposeMode b = ComposeMode (a <+> b)-    a *^ ComposeMode b = ComposeMode (auto a *^ b)-    ComposeMode a ^* b = ComposeMode (a ^* auto b)-    ComposeMode a ^/ b = ComposeMode (a ^/ auto b)-    ComposeMode a <**> ComposeMode b = ComposeMode (a <**> b)--instance (Mode f, Mode g) => Lifted (ComposeMode f g) where-    showsPrec1 n (ComposeMode a) = showsPrec1 n a-    ComposeMode a ==! ComposeMode b  = a ==! b-    compare1 (ComposeMode a) (ComposeMode b) = compare1 a b-    fromInteger1 = ComposeMode . auto . 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 . auto . 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-    erf1 (ComposeMode a) = ComposeMode (erf1 a)-    erfc1 (ComposeMode a) = ComposeMode (erfc1 a)-    normcdf1 (ComposeMode a) = ComposeMode (normcdf1 a)-    inverf1 (ComposeMode a) = ComposeMode (inverf1 a)-    inverfc1 (ComposeMode a) = ComposeMode (inverfc1 a)-    invnormcdf1 (ComposeMode a) = ComposeMode (invnormcdf1 a)--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-#if MIN_VERSION_base(4,4,0)-composeModeTyCon = mkTyCon3 "ad" "Numeric.AD.Internal.Composition" "ComposeMode"-#else-composeModeTyCon = mkTyCon "Numeric.AD.Internal.Composition.ComposeMode"-#endif-{-# 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-
src/Numeric/AD/Internal/Dense.hs view
@@ -1,9 +1,17 @@-{-# LANGUAGE Rank2Types, TypeFamilies, FlexibleContexts, UndecidableInstances, TemplateHaskell, DeriveDataTypeable, BangPatterns #-}--- {-# OPTIONS_HADDOCK hide, prune #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE UndecidableInstances #-}+{-# OPTIONS_HADDOCK not-home #-}+ ----------------------------------------------------------------------------- -- |--- Module      : Numeric.AD.Internal.Dense--- Copyright   : (c) Edward Kmett 2010+-- Copyright   : (c) Edward Kmett 2010-2014 -- License     : BSD3 -- Maintainer  : ekmett@gmail.com -- Stability   : experimental@@ -25,162 +33,155 @@ -----------------------------------------------------------------------------  module Numeric.AD.Internal.Dense-    ( Dense(..)-    , ds-    , ds'-    , vars-    , apply-    ) where+  ( Dense(..)+  , ds+  , ds'+  , vars+  , apply+  ) where -import Language.Haskell.TH+import Control.Monad (join)+import Data.Functor import Data.Typeable () import Data.Traversable (Traversable, mapAccumL) import Data.Data ()-import Numeric.AD.Internal.Types+import Data.Number.Erf import Numeric.AD.Internal.Combinators-import Numeric.AD.Internal.Classes import Numeric.AD.Internal.Identity+import Numeric.AD.Jacobian+import Numeric.AD.Mode -data Dense f a-    = Lift !a-    | Dense !a (f a)-    | Zero+data Dense f a s+  = 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"+type instance Scalar (Dense f a s) = a -ds :: f a -> AD (Dense f) a -> f a-ds _ (AD (Dense _ da)) = da+instance Show a => Show (Dense f a s) where+  showsPrec d (Lift a)    = showsPrec d a+  showsPrec d (Dense a _) = showsPrec d a+  showsPrec _ Zero        = showString "0"++ds :: f a -> Dense f a s -> f a+ds _ (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)+ds' :: Num a => f a -> Dense f a s -> (a, f a)+ds' _ (Dense a da) = (a, da)+ds' z (Lift a) = (a, z)+ds' z 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)+vars :: (Traversable f, Num a) => f a -> f (Dense f a s)+vars as = snd $ mapAccumL outer (0 :: Int) as where+  outer !i a = (i + 1, 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 :: (Traversable f, Num a) => (f (Dense f a s) -> 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+primal :: Num a => Dense f a s -> a+primal Zero = 0+primal (Lift a) = a+primal (Dense a _) = a -instance (Traversable f, Lifted (Dense f)) => Mode (Dense f) where-    auto = Lift-    zero = Zero+instance (Num a, Traversable f) => Mode (Dense f a s) where+  auto = 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      = auto (0 ** primal y)-    _    <**> Zero   = auto 1-    x    <**> Lift y = lift1 (**y) (\z -> y *^ z ** Id (y - 1)) x-    x    <**> y      = lift2_ (**) (\z xi yi -> (yi *! z /! xi, z *! log1 xi)) x y+  _ *^ 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 -    _ *^ 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+(<+>) :: (Traversable f, Num a) => Dense f a s -> Dense f a s -> Dense f a s+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 -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)+(<**>) :: (Traversable f, Floating a) => Dense f a s -> Dense f a s -> Dense f a s+Zero <**> y      = auto (0 ** primal y)+_    <**> Zero   = auto 1+x    <**> Lift y = lift1 (**y) (\z -> y *^ z ** Id (y - 1)) x+x    <**> y      = lift2_ (**) (\z xi yi -> (yi * z / xi, z * log xi)) x y -    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)+instance (Traversable f, Num a) => Jacobian (Dense f a s) where+  type D (Dense f a s) = Id a s+  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 a (fmap (dadb *) db)-        where-            a = f b-            Id dadb = df (Id a) (Id b)+  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) -    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+  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) -    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+  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 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+  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 -let f = varT (mkName "f") in-    deriveLifted-        (classP ''Traversable [f]:)-        (conT ''Dense `appT` f)+  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++#define BODY1(x)    (Traversable f, x)+#define BODY2(x,y) (Traversable f, x, y)+#define HEAD Dense f a s+#include "instances.h"
src/Numeric/AD/Internal/Forward.hs view
@@ -1,9 +1,17 @@-{-# LANGUAGE Rank2Types, TypeFamilies, DeriveDataTypeable, TemplateHaskell, UndecidableInstances, BangPatterns #-}--- {-# OPTIONS_HADDOCK hide, prune #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE TemplateHaskell #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE UndecidableInstances #-}+{-# OPTIONS_HADDOCK not-home #-}+ ----------------------------------------------------------------------------- -- |--- Module      :  Numeric.AD.Internal.Forward--- Copyright   :  (c) Edward Kmett 2010+-- Copyright   :  (c) Edward Kmett 2010-2014 -- License     :  BSD3 -- Maintainer  :  ekmett@gmail.com -- Stability   :  experimental@@ -15,190 +23,188 @@ -----------------------------------------------------------------------------  module Numeric.AD.Internal.Forward-    ( Forward(..)-    , tangent-    , bundle-    , unbundle-    , apply-    , bind-    , bind'-    , bindWith-    , bindWith'-    , transposeWith-    ) where+  ( Forward(..)+  , primal+  , 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 Control.Monad (join)+import Control.Applicative hiding ((<**>)) import Data.Data-import Control.Applicative-import Numeric.AD.Internal.Types-import Numeric.AD.Internal.Classes+import Data.Foldable (Foldable, toList)+import Data.Number.Erf+import Data.Traversable (Traversable, mapAccumL)+import Numeric.AD.Internal.Combinators import Numeric.AD.Internal.Identity+import Numeric.AD.Jacobian+import Numeric.AD.Mode +#ifdef HLINT {-# ANN module "HLint: ignore Reduce duplication" #-}+#endif  -- | 'Forward' mode AD-data Forward a+data Forward a s   = Forward !a a   | Lift !a   | Zero   deriving (Show, Data, Typeable) +type instance Scalar (Forward a s) = a+ -- | Calculate the 'tangent' using forward mode AD.-tangent :: Num a => AD Forward a -> a-tangent (AD (Forward _ da)) = da+tangent :: Num a => Forward a s -> a+tangent (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)+unbundle :: Num a => Forward a s -> (a, a)+unbundle (Forward a da) = (a, da)+unbundle Zero = (0,0)+unbundle (Lift a) = (a, 0) {-# INLINE unbundle #-} -bundle :: a -> a -> AD Forward a-bundle a da = AD (Forward a da)+bundle :: a -> a -> Forward a s+bundle = Forward {-# INLINE bundle #-} -apply :: Num a => (AD Forward a -> b) -> a -> b+apply :: Num a => (Forward a s -> 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-    auto = Lift-    zero = Zero+primal :: Num a => Forward a s -> a+primal (Forward a _) = a+primal (Lift a) = a+primal Zero = 0 -    isKnownZero Zero = True-    isKnownZero _    = False+instance Num a => Mode (Forward a s) where+  auto = Lift+  zero = Zero -    isKnownConstant Forward{} = False-    isKnownConstant _ = True+  isKnownZero Zero = True+  isKnownZero _    = False -    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)+  isKnownConstant Forward{} = False+  isKnownConstant _ = True -    Zero <**> y      = auto (0 ** primal y)-    _    <**> Zero   = auto 1-    x    <**> Lift y = lift1 (**y) (\z -> y *^ z ** Id (y - 1)) x-    x    <**> y      = lift2_ (**) (\z xi yi -> (yi *! z /! xi, z *! log1 xi)) x y+  a *^ Forward b db = Forward (a * b) (a * db)+  a *^ Lift b = Lift (a * b)+  _ *^ Zero = Zero -    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+  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+(<+>) :: Num a => Forward a s -> Forward a s -> Forward a s+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) -instance Lifted Forward => Jacobian Forward where-    type D Forward = Id+(<**>) :: Floating a => Forward a s -> Forward a s -> Forward a s+Zero <**> y      = auto (0 ** primal y)+_    <**> Zero   = auto 1+x    <**> Lift y = lift1 (**y) (\z -> y *^ z ** Id (y - 1)) x+x    <**> y      = lift2_ (**) (\z xi yi -> (yi * z / xi, z * log xi)) x y +instance Num a => Jacobian (Forward a s) where+  type D (Forward a s) = Id a s -    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)+  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 (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+  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+  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 (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+  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+#define HEAD Forward a s+#include "instances.h" -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 auto a)+bind :: (Traversable f, Num a) => (f (Forward a s) -> 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 auto a) -bind' :: (Traversable f, Num a) => (f (AD Forward a) -> b) -> f a -> (b, f b)-bind' f as = dropIx $ mapAccumL outer (0 :: Int, b0) as-    where-        outer (!i, _) _ = let b = f $ snd $ mapAccumL (inner i) (0 :: Int) as in ((i + 1, b), b)-        inner !i !j a = (j + 1, if i == j then bundle a 1 else auto a)-        b0 = f (auto <$> as)-        dropIx ((_,b),bs) = (b,bs)+bind' :: (Traversable f, Num a) => (f (Forward a s) -> 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 auto a)+  b0 = f (auto <$> as)+  dropIx ((_,b),bs) = (b,bs) -bindWith :: (Traversable f, Num a) => (a -> b -> c) -> (f (AD Forward a) -> b) -> f a -> f c-bindWith g f as = snd $ mapAccumL outer (0 :: Int) as-    where-        outer !i a = (i + 1, g a $ f $ snd $ mapAccumL (inner i) 0 as)-        inner !i !j a = (j + 1, if i == j then bundle a 1 else auto a)+bindWith :: (Traversable f, Num a) => (a -> b -> c) -> (f (Forward a s) -> 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 auto a) -bindWith' :: (Traversable f, Num a) => (a -> b -> c) -> (f (AD Forward a) -> b) -> f a -> (b, f c)-bindWith' g f as = dropIx $ mapAccumL outer (0 :: Int, b0) as-    where-        outer (!i, _) a = let b = f $ snd $ mapAccumL (inner i) (0 :: Int) as in ((i + 1, b), g a b)-        inner !i !j a = (j + 1, if i == j then bundle a 1 else auto a)-        b0 = f (auto <$> as)-        dropIx ((_,b),bs) = (b,bs)+bindWith' :: (Traversable f, Num a) => (a -> b -> c) -> (f (Forward a s) -> 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 auto a)+  b0 = f (auto <$> as)+  dropIx ((_,b),bs) = (b,bs)  -- we can't transpose arbitrary traversables, since we can't construct one out of whole cloth, and the outer -- 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-+transposeWith f as = snd . mapAccumL go xss0 where+  go xss b = (tail <$> xss, f b (head <$> xss))+  xss0 = toList <$> as
+ src/Numeric/AD/Internal/Forward/Double.hs view
@@ -0,0 +1,227 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE ConstraintKinds #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE TemplateHaskell #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE UndecidableInstances #-}+{-# OPTIONS_HADDOCK not-home #-}++-----------------------------------------------------------------------------+---- |+---- Copyright   :  (c) Edward Kmett 2010-2014+---- 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.Double+  ( ForwardDouble(..)+  , bundle+  , unbundle+  , apply+  , bind+  , bind'+  , bindWith+  , bindWith'+  , transposeWith+  ) where++import Control.Applicative hiding ((<**>))+import Control.Monad (join)+import Data.Foldable (Foldable, toList)+import Data.Function (on)+import Data.Number.Erf+import Data.Traversable (Traversable, mapAccumL)+import Numeric.AD.Internal.Combinators+import Numeric.AD.Internal.Identity+import Numeric.AD.Jacobian+import Numeric.AD.Mode++data ForwardDouble a = ForwardDouble { primal, tangent :: {-# UNPACK #-} !Double }+  deriving (Read, Show)++type instance Scalar (ForwardDouble s) = Double++unbundle :: ForwardDouble s -> (Double, Double)+unbundle (ForwardDouble a da) = (a, da)+{-# INLINE unbundle #-}++bundle :: Double -> Double -> ForwardDouble s+bundle = ForwardDouble+{-# INLINE bundle #-}++apply :: (ForwardDouble s -> b) -> Double -> b+apply f a = f (bundle a 1)+{-# INLINE apply #-}++instance Mode (ForwardDouble s) where+  auto = flip ForwardDouble 0+  zero = ForwardDouble 0 0++  isKnownZero (ForwardDouble 0 0) = True+  isKnownZero _ = False++  isKnownConstant (ForwardDouble _ 0) = True+  isKnownConstant _ = False++  a *^ ForwardDouble b db = ForwardDouble (a * b) (a * db)++  ForwardDouble a da ^* b = ForwardDouble (a * b) (da * b)++  ForwardDouble a da ^/ b = ForwardDouble (a / b) (da / b)++(<+>) :: ForwardDouble s -> ForwardDouble s -> ForwardDouble s+ForwardDouble a da <+> ForwardDouble b db = ForwardDouble (a + b) (da + db)++instance Jacobian (ForwardDouble s) where+  type D (ForwardDouble s) = Id Double s++  unary f (Id dadb) (ForwardDouble b db) = ForwardDouble (f b) (dadb * db)++  lift1 f df (ForwardDouble b db) = ForwardDouble (f b) (dadb * db) where+    Id dadb = df (Id b)++  lift1_ f df (ForwardDouble b db) = ForwardDouble a da where+    a = f b+    Id da = df (Id a) (Id b) ^* db++  binary f (Id dadb) (Id dadc) (ForwardDouble b db) (ForwardDouble c dc) = ForwardDouble (f b c) $ dadb * db + dc * dadc++  lift2 f df (ForwardDouble b db) (ForwardDouble c dc) = ForwardDouble a da where+    a = f b c+    (Id dadb, Id dadc) = df (Id b) (Id c)+    da = dadb * db + dc * dadc++  lift2_ f df (ForwardDouble b db) (ForwardDouble c dc) = ForwardDouble a da where+    a = f b c+    (Id dadb, Id dadc) = df (Id a) (Id b) (Id c)+    da = dadb * db + dc * dadc++instance Eq (ForwardDouble s) where+  (==)          = on (==) primal++instance Ord (ForwardDouble s) where+  compare       = on compare primal++instance Num (ForwardDouble s) where+  fromInteger 0  = zero+  fromInteger n = auto (fromInteger n)+  (+)          = (<+>) -- binary (+) 1 1+  (-)          = binary (-) (auto 1) (auto (-1)) -- TODO: <-> ? as it is, this might be pretty bad for Tower+  (*)          = lift2 (*) (\x y -> (y, x))+  negate       = lift1 negate (const (auto (-1)))+  abs          = lift1 abs signum+  signum a     = lift1 signum (const zero) a++instance Fractional (ForwardDouble s) where+  fromRational 0 = zero+  fromRational r = auto (fromRational r)+  x / y        = x * recip y+  recip        = lift1_ recip (const . negate . join (*))++instance Floating (ForwardDouble s) where+  pi       = auto pi+  exp      = lift1_ exp const+  log      = lift1 log recip+  logBase x y = log y / log x+  sqrt     = lift1_ sqrt (\z _ -> recip (auto 2 * z))+  ForwardDouble 0 0 ** ForwardDouble a _ = ForwardDouble (0 ** a) 0+  _ ** ForwardDouble 0 0                 = ForwardDouble 1 0+  x ** ForwardDouble y 0 = lift1 (**y) (\z -> y *^ z ** Id (y - 1)) x+  x ** y                 = lift2_ (**) (\z xi yi -> (yi * z / xi, z * log xi)) x y+  sin      = lift1 sin cos+  cos      = lift1 cos $ negate . sin+  tan      = lift1 tan $ recip . join (*) . cos+  asin     = lift1 asin $ \x -> recip (sqrt (auto 1 - join (*) x))+  acos     = lift1 acos $ \x -> negate (recip (sqrt (1 - join (*) x)))+  atan     = lift1 atan $ \x -> recip (1 + join (*) x)+  sinh     = lift1 sinh cosh+  cosh     = lift1 cosh sinh+  tanh     = lift1 tanh $ recip . join (*) . cosh+  asinh    = lift1 asinh $ \x -> recip (sqrt (1 + join (*) x))+  acosh    = lift1 acosh $ \x -> recip (sqrt (join (*) x - 1))+  atanh    = lift1 atanh $ \x -> recip (1 - join (*) x)++instance Enum (ForwardDouble s) where+  succ                 = lift1 succ (const 1)+  pred                 = lift1 pred (const 1)+  toEnum               = auto . toEnum+  fromEnum             = fromEnum . primal+  enumFrom a           = withPrimal a <$> enumFrom (primal a)+  enumFromTo a b       = withPrimal a <$> enumFromTo (primal a) (primal b)+  enumFromThen a b     = zipWith (fromBy a delta) [0..] $ enumFromThen (primal a) (primal b) where delta = b - a+  enumFromThenTo a b c = zipWith (fromBy a delta) [0..] $ enumFromThenTo (primal a) (primal b) (primal c) where delta = b - a++instance Real (ForwardDouble s) where+  toRational      = toRational . primal++instance RealFloat (ForwardDouble s) where+  floatRadix      = floatRadix . primal+  floatDigits     = floatDigits . primal+  floatRange      = floatRange . primal+  decodeFloat     = decodeFloat . primal+  encodeFloat m e = auto (encodeFloat m e)+  isNaN           = isNaN . primal+  isInfinite      = isInfinite . primal+  isDenormalized  = isDenormalized . primal+  isNegativeZero  = isNegativeZero . primal+  isIEEE          = isIEEE . primal+  exponent = exponent+  scaleFloat n = unary (scaleFloat n) (scaleFloat n 1)+  significand x =  unary significand (scaleFloat (- floatDigits x) 1) x+  atan2 = lift2 atan2 $ \vx vy -> let r = recip (join (*) vx + join (*) vy) in (vy * r, negate vx * r)++instance RealFrac (ForwardDouble s) where+  properFraction a = (w, a `withPrimal` pb) where+    pa = primal a+    (w, pb) = properFraction pa+  truncate = truncate . primal+  round    = round . primal+  ceiling  = ceiling . primal+  floor    = floor . primal++instance Erf (ForwardDouble s) where+  erf = lift1 erf $ \x -> (2 / sqrt pi) * exp (negate x * x)+  erfc = lift1 erfc $ \x -> ((-2) / sqrt pi) * exp (negate x * x)+  normcdf = lift1 normcdf $ \x -> ((-1) / sqrt pi) * exp (x * x * fromRational (- recip 2) / sqrt 2)++instance InvErf (ForwardDouble s) where+  inverf = lift1 inverfc $ \x -> recip $ (2 / sqrt pi) * exp (negate x * x)+  inverfc = lift1 inverfc $ \x -> recip $ negate (2 / sqrt pi) * exp (negate x * x)+  invnormcdf = lift1 invnormcdf $ \x -> recip $ ((-1) / sqrt pi) * exp (x * x * fromRational (- recip 2) / sqrt 2)++bind :: (Traversable f) => (f (ForwardDouble s) -> b) -> f Double -> 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 auto a)++bind' :: (Traversable f) => (f (ForwardDouble s) -> b) -> f Double -> (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 auto a)+  b0 = f (auto <$> as)+  dropIx ((_,b),bs) = (b,bs)++bindWith :: (Traversable f) => (Double -> b -> c) -> (f (ForwardDouble s) -> b) -> f Double -> 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 auto a)++bindWith' :: (Traversable f) => (Double -> b -> c) -> (f (ForwardDouble s) -> b) -> f Double -> (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 auto a)+  b0 = f (auto <$> as)+  dropIx ((_,b),bs) = (b,bs)++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
src/Numeric/AD/Internal/Identity.hs view
@@ -1,9 +1,14 @@-{-# LANGUAGE GeneralizedNewtypeDeriving, FlexibleContexts, FlexibleInstances, MultiParamTypeClasses, DeriveDataTypeable #-}-{-# OPTIONS_HADDOCK hide #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE TypeFamilies #-}+{-# OPTIONS_HADDOCK not-home #-}+ ----------------------------------------------------------------------------- -- |--- Module      :  Numeric.AD.Internal.Identity--- Copyright   :  (c) Edward Kmett 2010+-- Copyright   :  (c) Edward Kmett 2010-2014 -- License     :  BSD3 -- Maintainer  :  ekmett@gmail.com -- Stability   :  experimental@@ -11,146 +16,43 @@ -- ----------------------------------------------------------------------------- module Numeric.AD.Internal.Identity-    ( Id(..)-    , probe-    , unprobe-    , probed-    , unprobed-    ) where+  ( Id(..)+  , probe+  , unprobe+  , probed+  , unprobed+  ) where -import Control.Applicative import Data.Data (Data)-import Data.Foldable (Foldable, foldMap) import Data.Monoid import Data.Number.Erf import Data.Typeable (Typeable)-import Data.Traversable (Traversable, traverse)-import Numeric.AD.Internal.Classes-import Numeric.AD.Internal.Types--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)+import Numeric.AD.Mode -instance Foldable Id where-    foldMap f (Id a) = f a+newtype Id a s = Id { runId :: a } deriving+  (Eq, Ord, Show, Enum, Bounded, Num, Real, Fractional, Floating, RealFrac, RealFloat, Monoid, Data, Typeable, Erf, InvErf) -instance Traversable Id where-    traverse f (Id a) = Id <$> f a+type instance Scalar (Id a s) = a -instance Applicative Id where-    pure = Id-    Id f <*> Id a = Id (f a)+probe :: a -> Id a s+probe = Id -instance Monad Id where-    return = Id-    Id a >>= f = f a+unprobe :: Id a s -> a+unprobe = runId -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-    erf1 = erf-    erfc1 = erfc-    normcdf1 = normcdf-    inverf1 = inverf-    inverfc1 = inverfc-    invnormcdf1 = invnormcdf+pid :: Functor f => f a -> f (Id a s)+pid = fmap probe -instance Mode Id where-    auto = 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)+unpid :: Functor f => f (Id a s) -> f a+unpid = fmap unprobe -instance Primal Id where-    primal (Id a) = a+probed :: Functor f => f a -> f (Id a s)+probed = pid -instance Erf a => Erf (Id a) where-  erf = Id . erf . runId-  erfc = Id . erfc . runId-  normcdf = Id . normcdf . runId+unprobed :: Functor f => f (Id a s) -> f a+unprobed = unpid -instance InvErf a => InvErf (Id a) where-  inverf = Id . inverf . runId-  inverfc = Id . inverfc . runId-  invnormcdf = Id . invnormcdf . runId+instance Num a => Mode (Id a s) where+  auto = Id+  Id a ^* b = Id (a * b)+  a *^ Id b = Id (a * b)
− src/Numeric/AD/Internal/Jet.hs 
@@ -1,97 +0,0 @@-{-# LANGUAGE CPP, TypeOperators, 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--#ifndef MIN_VERSION_base-#define MIN_VERSION_base(x,y,z) 1-#endif--import Control.Applicative-import Data.Foldable-import Data.Traversable-import Data.Monoid-#if MIN_VERSION_base(4,4,0)-import Data.Typeable (Typeable1(..), TyCon, mkTyCon3, mkTyConApp)-#else-import Data.Typeable (Typeable1(..), TyCon, mkTyCon, mkTyConApp)-#endif-import Control.Comonad.Cofree--infixl 3 :----- | A 'Jet' is a tower of all (higher order) partial derivatives of a function------ At each step, a @'Jet' f@ is wrapped in another layer worth of @f@.------ > a :- f a :- f (f a) :- f (f (f a)) :- ...-data Jet f a = a :- Jet f (f a)---- | Used to sidestep the need for UndecidableInstances.-newtype Showable = Showable (Int -> String -> String)--instance Show Showable where-  showsPrec d (Showable f) = f d--showable :: Show a => a -> Showable-showable a = Showable (`showsPrec` 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---- | Take the tail of a 'Jet'.-tailJet :: Jet f a -> Jet f (f a)-tailJet (_ :- as) = as-{-# INLINE tailJet #-}---- | Take the head of a 'Jet'.-headJet :: Jet f a -> a-headJet (a :- _) = a-{-# INLINE headJet #-}---- | Construct a 'Jet' by unzipping the layers of a 'Cofree' 'Comonad'.-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 MIN_VERSION_base(4,4,0)-jetTyCon = mkTyCon3 "ad" "Numeric.AD.Internal.Jet" "Jet"-#else-jetTyCon = mkTyCon "Numeric.AD.Internal.Jet.Jet"-#endif-{-# NOINLINE jetTyCon #-}
src/Numeric/AD/Internal/Kahn.hs view
@@ -1,10 +1,18 @@-{-# LANGUAGE Rank2Types, MultiParamTypeClasses, FlexibleInstances, UndecidableInstances, ScopedTypeVariables, TemplateHaskell, TypeFamilies, DeriveDataTypeable, FunctionalDependencies #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE UndecidableInstances #-}+{-# OPTIONS_GHC -fno-full-laziness #-}+{-# OPTIONS_HADDOCK not-home #-} --- {-# OPTIONS_HADDOCK hide, prune #-} ----------------------------------------------------------------------------- -- |--- Module      :  Numeric.AD.Internal.Kahn--- Copyright   :  (c) Edward Kmett 2010+-- Copyright   :  (c) Edward Kmett 2010-2014 -- License     :  BSD3 -- Maintainer  :  ekmett@gmail.com -- Stability   :  experimental@@ -21,128 +29,143 @@ -----------------------------------------------------------------------------  module Numeric.AD.Internal.Kahn-    ( Kahn(..)-    , Tape(..)-    , partials-    , partialArray-    , partialMap-    , derivative-    , derivative'-    , vgrad, vgrad'-    , Grad(..)-    ) where+  ( Kahn(..)+  , Tape(..)+  , partials+  , partialArray+  , partialMap+  , derivative+  , derivative'+  , vgrad, vgrad'+  , Grad(..)+  , bind+  , unbind+  , unbindMap+  , unbindWith+  , unbindMapWithDefault+  , primal+  , var+  , varId+  ) where  import Prelude hiding (mapM) import Control.Applicative (Applicative(..),(<$>)) import Control.Monad.ST-import Control.Monad (forM_)+import Control.Monad hiding (mapM)+import Control.Monad.Trans.State import Data.List (foldl') import Data.Array.ST import Data.Array-import Data.IntMap (IntMap, fromListWith)+import Data.IntMap (IntMap, fromListWith, findWithDefault) import Data.Graph (Vertex, transposeG, Graph)+import Data.Number.Erf import Data.Reify (reifyGraph, MuRef(..)) import qualified Data.Reify.Graph as Reified import System.IO.Unsafe (unsafePerformIO)-import Language.Haskell.TH import Data.Data (Data)+import Data.Traversable (Traversable, mapM) import Data.Typeable (Typeable)-import Numeric.AD.Internal.Types-import Numeric.AD.Internal.Classes+import Numeric.AD.Internal.Combinators import Numeric.AD.Internal.Identity-import Numeric.AD.Internal.Var+import Numeric.AD.Jacobian+import Numeric.AD.Mode  -- | 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)+  = Zero+  | Lift !a+  | Var !a {-# UNPACK #-} !Int+  | Binary !a a a t t+  | Unary !a a t+  deriving (Show, Data, Typeable)  -- | @Kahn@ is a 'Mode' using reverse-mode automatic differentiation that provides fast 'diffFU', 'diff2FU', 'grad', 'grad2' and a fast 'jacobian' when you have a significantly smaller number of outputs than inputs.-newtype Kahn a = Kahn (Tape a (Kahn a)) deriving (Show, Typeable)+newtype Kahn a s = Kahn (Tape a (Kahn a s)) deriving (Show, Typeable) +type instance Scalar (Kahn a s) = a+ -- deriving instance (Data (Tape a (Kahn a)) => Data (Kahn a) -instance MuRef (Kahn a) where-    type DeRef (Kahn a) = Tape a+instance MuRef (Kahn a s) where+  type DeRef (Kahn a s) = Tape a -    mapDeRef _ (Kahn Zero) = pure Zero-    mapDeRef _ (Kahn (Lift a)) = pure (Lift a)-    mapDeRef _ (Kahn (Var a v)) = pure (Var a v)-    mapDeRef f (Kahn (Binary a dadb dadc b c)) = Binary a dadb dadc <$> f b <*> f c-    mapDeRef f (Kahn (Unary a dadb b)) = Unary a dadb <$> f b+  mapDeRef _ (Kahn Zero) = pure Zero+  mapDeRef _ (Kahn (Lift a)) = pure (Lift a)+  mapDeRef _ (Kahn (Var a v)) = pure (Var a v)+  mapDeRef f (Kahn (Binary a dadb dadc b c)) = Binary a dadb dadc <$> f b <*> f c+  mapDeRef f (Kahn (Unary a dadb b)) = Unary a dadb <$> f b -instance Lifted Kahn => Mode Kahn where-    isKnownZero (Kahn Zero) = True-    isKnownZero _    = False+instance Num a => Mode (Kahn a s) where+  isKnownZero (Kahn Zero) = True+  isKnownZero _    = False -    isKnownConstant (Kahn Zero) = True-    isKnownConstant (Kahn (Lift _)) = True-    isKnownConstant _ = False+  isKnownConstant (Kahn Zero) = True+  isKnownConstant (Kahn (Lift _)) = True+  isKnownConstant _ = False -    auto a = Kahn (Lift a)-    zero   = Kahn Zero-    (<+>)  = binary (+) one one-    a *^ b = lift1 (a *) (\_ -> auto a) b-    a ^* b = lift1 (* b) (\_ -> auto b) a-    a ^/ b = lift1 (/ b) (\_ -> auto (recip b)) a+  auto a = Kahn (Lift a)+  zero   = Kahn Zero+  a *^ b = lift1 (a *) (\_ -> auto a) b+  a ^* b = lift1 (* b) (\_ -> auto b) a+  a ^/ b = lift1 (/ b) (\_ -> auto (recip b)) a -    Kahn Zero <**> y                = auto (0 ** primal y)-    _            <**> Kahn Zero     = auto 1-    x            <**> Kahn (Lift y) = lift1 (**y) (\z -> y *^ z ** Id (y-1)) x-    x            <**> y                = lift2_ (**) (\z xi yi -> (yi *! z /! xi, z *! log1 xi)) x y+(<+>) :: Num a => Kahn a s -> Kahn a s -> Kahn a s+(<+>)  = binary (+) 1 1 -instance Primal Kahn where-    primal (Kahn Zero) = 0-    primal (Kahn (Lift a)) = a-    primal (Kahn (Var a _)) = a-    primal (Kahn (Binary a _ _ _ _)) = a-    primal (Kahn (Unary a _ _)) = a+(<**>) :: Floating a => Kahn a s -> Kahn a s -> Kahn a s+Kahn Zero <**> y             = auto (0 ** primal y)+_         <**> Kahn Zero     = auto 1+x         <**> Kahn (Lift y) = lift1 (**y) (\z -> y *^ z ** Id (y-1)) x+x         <**> y             = lift2_ (**) (\z xi yi -> (yi * z / xi, z * xi)) x y -instance Lifted Kahn => Jacobian Kahn where-    type D Kahn = Id+primal :: Num a => Kahn a s -> a+primal (Kahn Zero) = 0+primal (Kahn (Lift a)) = a+primal (Kahn (Var a _)) = a+primal (Kahn (Binary a _ _ _ _)) = a+primal (Kahn (Unary a _ _)) = a -    unary f _         (Kahn Zero)     = Kahn (Lift (f 0))-    unary f _         (Kahn (Lift a)) = Kahn (Lift (f a))-    unary f (Id dadb) b                  = Kahn (Unary (f (primal b)) dadb b)+instance Num a => Jacobian (Kahn a s) where+  type D (Kahn a s) = Id a s -    lift1 f df b = unary f (df (Id pb)) b-        where pb = primal b+  unary f _         (Kahn Zero)     = Kahn (Lift (f 0))+  unary f _         (Kahn (Lift a)) = Kahn (Lift (f a))+  unary f (Id dadb) b                  = Kahn (Unary (f (primal b)) dadb b) -    lift1_ f df b = unary (const a) (df (Id a) (Id pb)) b-        where pb = primal b-              a = f pb+  lift1 f df b = unary f (df (Id pb)) b where+    pb = primal b -    binary f _         _         (Kahn Zero)     (Kahn Zero)     = Kahn (Lift (f 0 0))-    binary f _         _         (Kahn Zero)     (Kahn (Lift c)) = Kahn (Lift (f 0 c))-    binary f _         _         (Kahn (Lift b)) (Kahn Zero)     = Kahn (Lift (f b 0))-    binary f _         _         (Kahn (Lift b)) (Kahn (Lift c)) = Kahn (Lift (f b c))-    binary f _         (Id dadc) (Kahn Zero)     c                  = Kahn (Unary (f 0 (primal c)) dadc c)-    binary f _         (Id dadc) (Kahn (Lift b)) c                  = Kahn (Unary (f b (primal c)) dadc c)-    binary f (Id dadb) _         b                  (Kahn Zero)     = Kahn (Unary (f (primal b) 0) dadb b)-    binary f (Id dadb) _         b                  (Kahn (Lift c)) = Kahn (Unary (f (primal b) c) dadb b)-    binary f (Id dadb) (Id dadc) b                  c                  = Kahn (Binary (f (primal b) (primal c)) dadb dadc b c)+  lift1_ f df b = unary (const a) (df (Id a) (Id pb)) b where+    pb = primal b+    a = f pb -    lift2 f df b c = binary f dadb dadc b c-        where (dadb, dadc) = df (Id (primal b)) (Id (primal c))+  binary f _         _         (Kahn Zero)     (Kahn Zero)     = Kahn (Lift (f 0 0))+  binary f _         _         (Kahn Zero)     (Kahn (Lift c)) = Kahn (Lift (f 0 c))+  binary f _         _         (Kahn (Lift b)) (Kahn Zero)     = Kahn (Lift (f b 0))+  binary f _         _         (Kahn (Lift b)) (Kahn (Lift c)) = Kahn (Lift (f b c))+  binary f _         (Id dadc) (Kahn Zero)     c                  = Kahn (Unary (f 0 (primal c)) dadc c)+  binary f _         (Id dadc) (Kahn (Lift b)) c                  = Kahn (Unary (f b (primal c)) dadc c)+  binary f (Id dadb) _         b                  (Kahn Zero)     = Kahn (Unary (f (primal b) 0) dadb b)+  binary f (Id dadb) _         b                  (Kahn (Lift c)) = Kahn (Unary (f (primal b) c) dadb b)+  binary f (Id dadb) (Id dadc) b                  c                  = Kahn (Binary (f (primal b) (primal c)) dadb dadc b c) -    lift2_ f df b c = binary (\_ _ -> 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)+  lift2 f df b c = binary f dadb dadc b c where+    (dadb, dadc) = df (Id (primal b)) (Id (primal c)) -deriveLifted id (conT ''Kahn)+  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) -derivative :: Num a => AD Kahn a -> a+#define HEAD Kahn a s+#include <instances.h>++derivative :: Num a => Kahn a s -> a derivative = sum . map snd . partials {-# INLINE derivative #-} -derivative' :: Num a => AD Kahn a -> (a, a)+derivative' :: Num a => Kahn a s -> (a, a) derivative' r = (primal r, derivative r) {-# INLINE derivative' #-} @@ -150,15 +173,15 @@ backPropagate :: Num a => (Vertex -> (Tape a Int, Int, [Int])) -> STArray s Int a -> Vertex -> ST s () backPropagate vmap ss v = case node of   Unary _ g b -> do-      da <- readArray ss i-      db <- readArray ss b-      writeArray ss b (db + g*da)+    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)+    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@@ -166,99 +189,106 @@  topSortAcyclic :: Graph -> [Vertex] topSortAcyclic g = reverse $ runST $ do-    del <- newArray (bounds g) False :: ST s (STUArray s Int Bool)-    let tg = transposeG g-        starters = [ n | (n, []) <- assocs tg ]-        loop [] rs = return rs-        loop (n:ns) rs = do-            writeArray del n True-            let add [] = return ns-                add (m:ms) = do-                    b <- ok (tg!m)-                    ms' <- add ms-                    return $ if b then m : ms' else ms'-                ok [] = return True-                ok (x:xs) = do b <- readArray del x; if b then ok xs else return False-            ns' <- add (g!n)-            loop ns' (n : rs)-    loop starters []+  del <- newArray (bounds g) False :: ST s (STUArray s Int Bool)+  let tg = transposeG g+      starters = [ n | (n, []) <- assocs tg ]+      loop [] rs = return rs+      loop (n:ns) rs = do+        writeArray del n True+        let add [] = return ns+            add (m:ms) = do+              b <- ok (tg!m)+              ms' <- add ms+              return $ if b then m : ms' else ms'+            ok [] = return True+            ok (x:xs) = do b <- readArray del x; if b then ok xs else return False+        ns' <- add (g!n)+        loop ns' (n : rs)+  loop starters []  -- | This returns a list of contributions to the partials. -- The variable ids returned in the list are likely /not/ unique!-{-# SPECIALIZE partials :: AD Kahn Double -> [(Int, Double)] #-}-partials :: forall a . Num a => AD Kahn a -> [(Int, a)]-partials (AD tape) = [ let v = sensitivities ! ix in seq v (ident, v) | (ix, Var _ ident) <- xs ]-    where-        Reified.Graph xs start = unsafePerformIO $ reifyGraph tape-        g = array xsBounds [ (i, successors t) | (i, t) <- xs ]-        vertexMap = array xsBounds xs-        vmap i = (vertexMap ! i, i, [])-        xsBounds = sbounds xs+{-# SPECIALIZE partials :: Kahn Double s -> [(Int, Double)] #-}+partials :: forall s a . Num a => Kahn a s -> [(Int, a)]+partials tape = [ let v = sensitivities ! ix in seq v (ident, v) | (ix, Var _ ident) <- xs ] where+  Reified.Graph xs start = unsafePerformIO $ reifyGraph tape+  g = array xsBounds [ (i, successors t) | (i, t) <- xs ]+  vertexMap = array xsBounds xs+  vmap i = (vertexMap ! i, i, [])+  xsBounds = sbounds xs -        sensitivities = runSTArray $ do-            ss <- newArray xsBounds 0-            writeArray ss start 1-            forM_ (topSortAcyclic g) $-                backPropagate vmap ss-            return ss+  sensitivities = runSTArray $ do+    ss <- newArray xsBounds 0+    writeArray ss start 1+    forM_ (topSortAcyclic g) $+      backPropagate vmap ss+    return ss -        sbounds ((a,_):as) = foldl' (\(lo,hi) (b,_) -> let lo' = min lo b; hi' = max hi b in lo' `seq` hi' `seq` (lo', hi')) (a,a) as-        sbounds _ = undefined -- the graph can't be empty, it contains the output node!+  sbounds ((a,_):as) = foldl' (\(lo,hi) (b,_) -> let lo' = min lo b; hi' = max hi b in lo' `seq` hi' `seq` (lo', hi')) (a,a) as+  sbounds _ = undefined -- the graph can't be empty, it contains the output node! -        successors :: Tape a t -> [t]-        successors (Unary _ _ b) = [b]-        successors (Binary _ _ _ b c) = [b,c]-        successors _ = []+  successors :: Tape a t -> [t]+  successors (Unary _ _ b) = [b]+  successors (Binary _ _ _ b c) = [b,c]+  successors _ = []  -- | Return an 'Array' of 'partials' given bounds for the variable IDs.-partialArray :: Num a => (Int, Int) -> AD Kahn a -> Array Int a+partialArray :: Num a => (Int, Int) -> Kahn a s -> Array Int a partialArray vbounds tape = accumArray (+) 0 vbounds (partials tape) {-# INLINE partialArray #-}  -- | Return an 'IntMap' of sparse partials-partialMap :: Num a => AD Kahn a -> IntMap a+partialMap :: Num a => Kahn a s -> IntMap a partialMap = fromListWith (+) . partials {-# INLINE partialMap #-} --- A simple fresh variable supply monad-newtype S a = S { runS :: Int -> (a,Int) }-instance Monad S where-    return a = S (\s -> (a,s))-    S g >>= f = S (\s -> let (a,s') = g s in runS (f a) s')--instance Var Kahn where-    var a v = Kahn (Var a v)-    varId (Kahn (Var _ v)) = v-    varId _ = error "varId: not a Var"- class Num a => Grad i o o' a | i -> a o o', o -> a i o', o' -> a i o where-    pack :: i -> [AD Kahn a] -> AD Kahn a-    unpack :: ([a] -> [a]) -> o-    unpack' :: ([a] -> (a, [a])) -> o'+  pack :: i -> [Kahn a ()] -> Kahn a ()+  unpack :: ([a] -> [a]) -> o+  unpack' :: ([a] -> (a, [a])) -> o' -instance Num a => Grad (AD Kahn a) [a] (a, [a]) a where-    pack i _ = i-    unpack f = f []-    unpack' f = f []+instance Num a => Grad (Kahn a ()) [a] (a, [a]) a where+  pack i _ = i+  unpack f = f []+  unpack' f = f [] -instance Grad i o o' a => Grad (AD Kahn a -> i) (a -> o) (a -> o') a where-    pack f (a:as) = pack (f a) as-    pack _ [] = error "Grad.pack: logic error"-    unpack f a = unpack (f . (a:))-    unpack' f a = unpack' (f . (a:))+instance Grad i o o' a => Grad (Kahn a () -> i) (a -> o) (a -> o') a where+  pack f (a:as) = pack (f a) as+  pack _ [] = error "Grad.pack: logic error"+  unpack f a = unpack (f . (a:))+  unpack' f a = unpack' (f . (a:))  vgrad :: Grad i o o' a => i -> o-vgrad i = unpack (unsafeGrad (pack i))-    where-        unsafeGrad f as = unbind vs (partialArray bds $ f vs)-            where-                (vs,bds) = bind as+vgrad 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+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 +var :: a -> Int -> Kahn a s+var a v = Kahn (Var a v)++varId :: Kahn a s -> Int+varId (Kahn (Var _ v)) = v+varId _ = error "varId: not a Var"++bind :: Traversable f => f a -> (f (Kahn a s), (Int,Int))+bind xs = (r,(0,hi)) where+  (r,hi) = runState (mapM freshVar xs) 0+  freshVar a = state $ \s -> let s' = s + 1 in s' `seq` (var a s, s')++unbind :: Functor f => f (Kahn a s) -> Array Int a -> f a+unbind xs ys = fmap (\v -> ys ! varId v) xs++unbindWith :: (Functor f, Num a) => (a -> b -> c) -> f (Kahn a s) -> Array Int b -> f c+unbindWith f xs ys = fmap (\v -> f (primal v) (ys ! varId v)) xs++unbindMap :: (Functor f, Num a) => f (Kahn a s) -> IntMap a -> f a+unbindMap xs ys = fmap (\v -> findWithDefault 0 (varId v) ys) xs++unbindMapWithDefault :: (Functor f, Num a) => b -> (a -> b -> c) -> f (Kahn a s) -> IntMap b -> f c+unbindMapWithDefault z f xs ys = fmap (\v -> f (primal v) $ findWithDefault z (varId v) ys) xs
+ src/Numeric/AD/Internal/On.hs view
@@ -0,0 +1,49 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE StandaloneDeriving #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-# OPTIONS_HADDOCK not-home #-}++-----------------------------------------------------------------------------+-- |+-- Copyright   :  (c) Edward Kmett 2010-2014+-- License     :  BSD3+-- Maintainer  :  ekmett@gmail.com+-- Stability   :  experimental+-- Portability :  GHC only+--+-----------------------------------------------------------------------------++module Numeric.AD.Internal.On+  ( On(..)+  ) where++import Data.Number.Erf+import Data.Data+import Numeric.AD.Mode++#ifdef HLINT+#endif++------------------------------------------------------------------------------+-- On+------------------------------------------------------------------------------++-- | The composition of two AD modes is an AD mode in its own right+newtype On t = On { off :: t } deriving+  ( Eq, Enum, Ord, Bounded+  , Num, Real, Fractional+  , RealFrac, Floating, Erf+  , InvErf, RealFloat, Typeable+  )++type instance Scalar (On t) = Scalar (Scalar t)++instance (Mode t, Mode (Scalar t)) => Mode (On t) where+  auto = On . auto . auto+  a *^ On b = On (auto a *^ b)+  On a ^* b = On (a ^* auto b)
src/Numeric/AD/Internal/Reverse.hs view
@@ -1,9 +1,20 @@-{-# LANGUAGE CPP, Rank2Types, MultiParamTypeClasses, FlexibleInstances, UndecidableInstances, ScopedTypeVariables, TemplateHaskell, GADTs, TypeFamilies, DeriveDataTypeable, FlexibleContexts #-}--- {-# OPTIONS_HADDOCK hide, prune #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TemplateHaskell #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE UndecidableInstances #-}+{-# OPTIONS_GHC -fno-full-laziness #-}+{-# OPTIONS_HADDOCK not-home #-}+ ----------------------------------------------------------------------------- -- |--- Module      :  Numeric.AD.Internal.Reverse--- Copyright   :  (c) Edward Kmett 2012+-- Copyright   :  (c) Edward Kmett 2012-2014 -- License     :  BSD3 -- Maintainer  :  ekmett@gmail.com -- Stability   :  experimental@@ -21,37 +32,51 @@ -----------------------------------------------------------------------------  module Numeric.AD.Internal.Reverse-    ( Reverse(..)-    , Tape(..)-    , Head(..)-    , Cells(..)-    , reifyTape-    , partials-    , partialArrayOf-    , partialMapOf-    , derivativeOf-    , derivativeOf'-    ) where+  ( Reverse(..)+  , Tape(..)+  , Head(..)+  , Cells(..)+  , reifyTape+  , partials+  , partialArrayOf+  , partialMapOf+  , derivativeOf+  , derivativeOf'+  , bind+  , unbind+  , unbindMap+  , unbindWith+  , unbindMapWithDefault+  , var+  , varId+  , primal+  ) where +import Data.Functor+import Control.Monad hiding (mapM) import Control.Monad.ST+import Control.Monad.Trans.State import Data.Array.ST import Data.Array import Data.Array.Unsafe as Unsafe import Data.IORef-import Data.IntMap (IntMap, fromDistinctAscList)+import Data.IntMap (IntMap, fromDistinctAscList, findWithDefault)+import Data.Number.Erf import Data.Proxy import Data.Reflection+import Data.Traversable (Traversable, mapM) import Data.Typeable-import Language.Haskell.TH hiding (reify)-import Numeric.AD.Internal.Types-import Numeric.AD.Internal.Classes+import Numeric.AD.Internal.Combinators import Numeric.AD.Internal.Identity-import Numeric.AD.Internal.Var+import Numeric.AD.Jacobian+import Numeric.AD.Mode import Prelude hiding (mapM) import System.IO.Unsafe (unsafePerformIO) import Unsafe.Coerce +#ifdef HLINT {-# ANN module "HLint: ignore Reduce duplication" #-}+#endif  -- evil untyped tape #ifndef HLINT@@ -89,25 +114,27 @@  -- | This is used to create a new entry on the chain given a unary function, its derivative with respect to its input, -- the variable ID of its input, and the value of its input. Used by 'unary' and 'binary' internally.-unarily :: forall s a. Reifies s Tape => (a -> a) -> a -> Int -> a -> Reverse s a+unarily :: forall s a. Reifies s Tape => (a -> a) -> a -> Int -> a -> Reverse a s unarily f di i b = Reverse (unsafePerformIO (modifyTape (Proxy :: Proxy s) (un i di))) $! f b {-# INLINE unarily #-}  -- | This is used to create a new entry on the chain given a binary function, its derivatives with respect to its inputs, -- their variable IDs and values. Used by 'binary' internally.-binarily :: forall s a. Reifies s Tape => (a -> a -> a) -> a -> a -> Int -> a -> Int -> a -> Reverse s a+binarily :: forall s a. Reifies s Tape => (a -> a -> a) -> a -> a -> Int -> a -> Int -> a -> Reverse a s binarily f di dj i b j c = Reverse (unsafePerformIO (modifyTape (Proxy :: Proxy s) (bin i j di dj))) $! f b c {-# INLINE binarily #-}  #ifndef HLINT-data Reverse s a where-  Zero :: Reverse s a-  Lift :: a -> Reverse s a-  Reverse :: {-# UNPACK #-} !Int -> a -> Reverse s a+data Reverse a s where+  Zero :: Reverse a s+  Lift :: a -> Reverse a s+  Reverse :: {-# UNPACK #-} !Int -> a -> Reverse a s   deriving (Show, Typeable) #endif -instance (Reifies s Tape, Lifted (Reverse s)) => Mode (Reverse s) where+type instance Scalar (Reverse a s) = a++instance (Num a, Reifies s Tape) => Mode (Reverse a s) where   isKnownZero Zero = True   isKnownZero _    = False @@ -116,66 +143,70 @@    auto = Lift   zero = Zero-  (<+>)  = binary (+) one one   a *^ b = lift1 (a *) (\_ -> auto a) b   a ^* b = lift1 (* b) (\_ -> auto b) a   a ^/ b = lift1 (/ b) (\_ -> auto (recip b)) a -  Zero <**> y      = auto (0 ** primal y)-  _    <**> Zero   = auto 1-  x    <**> Lift y = lift1 (**y) (\z -> y *^ z ** Id (y - 1)) x-  x    <**> y      = lift2_ (**) (\z xi yi -> (yi *! z /! xi, z *! log1 xi)) x y+(<+>) :: (Reifies s Tape, Num a) => Reverse a s -> Reverse a s -> Reverse a s+(<+>)  = binary (+) 1 1 -instance Primal (Reverse s) where-    primal Zero = 0-    primal (Lift a) = a-    primal (Reverse _ a) = a+(<**>) :: (Reifies s Tape, Floating a) => Reverse a s -> Reverse a s -> Reverse a s+Zero <**> y      = auto (0 ** primal y)+_    <**> Zero   = auto 1+x    <**> Lift y = lift1 (**y) (\z -> y *^ z ** Id (y - 1)) x+x    <**> y      = lift2_ (**) (\z xi yi -> (yi * z / xi, z * log xi)) x y -instance (Reifies s Tape, Lifted (Reverse s)) => Jacobian (Reverse s) where-    type D (Reverse s) = Id+primal :: Num a => Reverse a s -> a+primal Zero = 0+primal (Lift a) = a+primal (Reverse _ a) = a -    unary f _         (Zero)   = Lift (f 0)-    unary f _         (Lift a) = Lift (f a)-    unary f (Id dadi) (Reverse i b) = unarily f dadi i b+instance (Reifies s Tape, Num a) => Jacobian (Reverse a s) where+  type D (Reverse a s) = Id a s -    lift1 f df b = unary f (df (Id pb)) b-        where pb = primal b+  unary f _         (Zero)   = Lift (f 0)+  unary f _         (Lift a) = Lift (f a)+  unary f (Id dadi) (Reverse i b) = unarily f dadi i b -    lift1_ f df b = unary (const a) (df (Id a) (Id pb)) b-        where pb = primal b-              a = f pb+  lift1 f df b = unary f (df (Id pb)) b where+    pb = primal 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)+  lift1_ f df b = unary (const a) (df (Id a) (Id pb)) b where+    pb = primal b+    a = f pb -    binary f _         (Id dadc) Zero        (Reverse i c) = unarily (f 0) dadc i c-    binary f _         (Id dadc) (Lift b)    (Reverse i c) = unarily (f b) dadc i c-    binary f (Id dadb) _         (Reverse i b) Zero        = unarily (`f` 0) dadb i b-    binary f (Id dadb) _         (Reverse i b) (Lift c)    = unarily (`f` c) dadb i b-    binary f (Id dadb) (Id dadc) (Reverse i b) (Reverse j c) = binarily f dadb dadc i b j c+  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) -    lift2 f df b c = binary f dadb dadc b c-        where (dadb, dadc) = df (Id (primal b)) (Id (primal c))+  binary f _         (Id dadc) Zero        (Reverse i c) = unarily (f 0) dadc i c+  binary f _         (Id dadc) (Lift b)    (Reverse i c) = unarily (f b) dadc i c+  binary f (Id dadb) _         (Reverse i b) Zero        = unarily (`f` 0) dadb i b+  binary f (Id dadb) _         (Reverse i b) (Lift c)    = unarily (`f` c) dadb i b+  binary f (Id dadb) (Id dadc) (Reverse i b) (Reverse j c) = binarily f dadb dadc i b j c -    lift2_ f df b c = binary (\_ _ -> 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)+  lift2 f df b c = binary f dadb dadc b c where+    (dadb, dadc) = df (Id (primal b)) (Id (primal c)) -let s = varT (mkName "s") in-  deriveLifted (classP ''Reifies [s, conT ''Tape] :) (conT ''Reverse `appT` s)+  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) --- | Helper that extracts the derivative of a chain when the chain was constructed with one variable.-derivativeOf :: (Reifies s Tape, Num a) => Proxy s -> AD (Reverse s) a -> a+#define BODY1(x) (Reifies s Tape,x)+#define BODY2(x,y) (Reifies s Tape,x,y)+#define HEAD Reverse a s+#include "instances.h"++-- | Helper that extracts the derivative of a chain when the chain was constructed with 1 variable.+derivativeOf :: (Reifies s Tape, Num a) => Proxy s -> Reverse a s -> a derivativeOf _ = sum . partials {-# INLINE derivativeOf #-} --- | Helper that extracts both the primal and derivative of a chain when the chain was constructed with one variable.-derivativeOf' :: (Reifies s Tape, Num a) => Proxy s -> AD (Reverse s) a -> (a, a)+-- | Helper that extracts both the primal and derivative of a chain when the chain was constructed with 1 variable.+derivativeOf' :: (Reifies s Tape, Num a) => Proxy s -> Reverse a s -> (a, a) derivativeOf' p r = (primal r, derivativeOf p r) {-# INLINE derivativeOf' #-} @@ -196,27 +227,27 @@   (backPropagate $! k - 1) xs ss  -- | Extract the partials from the current chain for a given AD variable.-{-# SPECIALIZE partials :: Reifies s Tape => AD (Reverse s) Double -> [Double] #-}-partials :: forall s a. (Reifies s Tape, Num a) => AD (Reverse s) a -> [a]-partials (AD Zero)        = []-partials (AD (Lift _))    = []-partials (AD (Reverse k _)) = map (sensitivities !) [0..vs] where-   Head n t = unsafePerformIO $ readIORef (getTape (reflect (Proxy :: Proxy s)))-   tk = dropCells (n - k) t-   (vs,sensitivities) = runST $ do-     ss <- newArray (0, k) 0-     writeArray ss k 1-     v <- backPropagate k tk ss-     as <- Unsafe.unsafeFreeze ss-     return (v, as)+{-# SPECIALIZE partials :: Reifies s Tape => Reverse Double s -> [Double] #-}+partials :: forall s a. (Reifies s Tape, Num a) => Reverse a s -> [a]+partials Zero        = []+partials (Lift _)    = []+partials (Reverse k _) = map (sensitivities !) [0..vs] where+  Head n t = unsafePerformIO $ readIORef (getTape (reflect (Proxy :: Proxy s)))+  tk = dropCells (n - k) t+  (vs,sensitivities) = runST $ do+    ss <- newArray (0, k) 0+    writeArray ss k 1+    v <- backPropagate k tk ss+    as <- Unsafe.unsafeFreeze ss+    return (v, as)  -- | Return an 'Array' of 'partials' given bounds for the variable IDs.-partialArrayOf :: (Reifies s Tape, Num a) => Proxy s -> (Int, Int) -> AD (Reverse s) a -> Array Int a+partialArrayOf :: (Reifies s Tape, Num a) => Proxy s -> (Int, Int) -> Reverse a s -> Array Int a partialArrayOf _ vbounds = accumArray (+) 0 vbounds . zip [0..] . partials {-# INLINE partialArrayOf #-}  -- | Return an 'IntMap' of sparse partials-partialMapOf :: (Reifies s Tape, Num a) => Proxy s -> AD (Reverse s) a -> IntMap a+partialMapOf :: (Reifies s Tape, Num a) => Proxy s -> Reverse a s -> IntMap a partialMapOf _ = fromDistinctAscList . zip [0..] . partials {-# INLINE partialMapOf #-} @@ -227,7 +258,26 @@   return (reify (Tape h) k) {-# NOINLINE reifyTape #-} -instance Var (Reverse s) where-    var a v = Reverse v a-    varId (Reverse v _) = v-    varId _ = error "varId: not a Var"+var :: a -> Int -> Reverse a s+var a v = Reverse v a++varId :: Reverse a s -> Int+varId (Reverse v _) = v+varId _ = error "varId: not a Var"++bind :: Traversable f => f a -> (f (Reverse a s), (Int,Int))+bind xs = (r,(0,hi)) where+  (r,hi) = runState (mapM freshVar xs) 0+  freshVar a = state $ \s -> let s' = s + 1 in s' `seq` (var a s, s')++unbind :: Functor f => f (Reverse a s) -> Array Int a -> f a+unbind xs ys = fmap (\v -> ys ! varId v) xs++unbindWith :: (Functor f, Num a) => (a -> b -> c) -> f (Reverse a s) -> Array Int b -> f c+unbindWith f xs ys = fmap (\v -> f (primal v) (ys ! varId v)) xs++unbindMap :: (Functor f, Num a) => f (Reverse a s) -> IntMap a -> f a+unbindMap xs ys = fmap (\v -> findWithDefault 0 (varId v) ys) xs++unbindMapWithDefault :: (Functor f, Num a) => b -> (a -> b -> c) -> f (Reverse a s) -> IntMap b -> f c+unbindMapWithDefault z f xs ys = fmap (\v -> f (primal v) $ findWithDefault z (varId v) ys) xs
src/Numeric/AD/Internal/Sparse.hs view
@@ -1,35 +1,48 @@-{-# LANGUAGE BangPatterns, TemplateHaskell, TypeFamilies, TypeOperators, FlexibleContexts, UndecidableInstances, DeriveDataTypeable, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE UndecidableInstances #-} {-# OPTIONS_GHC -fno-warn-name-shadowing #-}+{-# OPTIONS_HADDOCK not-home #-}+ module Numeric.AD.Internal.Sparse-    ( Index(..)-    , emptyIndex-    , addToIndex-    , indices-    , Sparse(..)-    , apply-    , vars-    , d, d', ds-    , skeleton-    , spartial-    , partial-    , vgrad-    , vgrad'-    , vgrads-    , Grad(..)-    , Grads(..)-    ) where+  ( 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 Control.Monad (join) import Data.Data-import Data.Typeable ()-import qualified Data.IntMap as IntMap import Data.IntMap (IntMap, mapWithKey, unionWith, findWithDefault, toAscList, singleton, insertWith, lookup)+import qualified Data.IntMap as IntMap+import Data.Number.Erf import Data.Traversable-import Language.Haskell.TH+import Data.Typeable ()+import Numeric.AD.Internal.Combinators+import Numeric.AD.Jacobian+import Numeric.AD.Mode  newtype Index = Index (IntMap Int) @@ -50,32 +63,33 @@ -- 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))+data Sparse a s+  = Sparse !a (IntMap (Sparse a s))   | Zero   deriving (Show, Data, Typeable) +type instance Scalar (Sparse a s) = a+ -- | 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 :: Num a => Sparse a s -> Int -> Sparse a s -> Sparse a s 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))+  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 $ auto 1)+vars :: (Traversable f, Num a) => f a -> f (Sparse a s)+vars = snd . mapAccumL var 0 where+  var !n a = (n + 1, Sparse a $ singleton n $ auto 1) {-# INLINE vars #-} -apply :: (Traversable f, Num a) => (f (AD Sparse a) -> b) -> f a -> b+apply :: (Traversable f, Num a) => (f (Sparse a s) -> b) -> f a -> b apply f = f . vars {-# INLINE apply #-} @@ -83,24 +97,23 @@ 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+d :: (Traversable f, Num a) => f b -> Sparse a s -> f a+d fs (Zero) = 0 <$ fs+d fs (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)+d' :: (Traversable f, Num a) => f a -> Sparse a s -> (a, f a)+d' fs Zero = (0, 0 <$ fs)+d' fs (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+ds :: (Traversable f, Num a) => f b -> Sparse a s -> Cofree f a+ds fs Zero = r where r = 0 :< (r <$ fs)+ds fs (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 #-}  {-@@ -129,128 +142,131 @@ {-# INLINE vds #-} -} -partial :: Num a => [Int] -> Sparse a -> a+partial :: Num a => [Int] -> Sparse a s -> a partial []     (Sparse a _)  = a partial (n:ns) (Sparse _ da) = partial ns $ findWithDefault (auto 0) n da partial _      Zero          = 0 {-# INLINE partial #-} -spartial :: Num a => [Int] -> Sparse a -> Maybe a+spartial :: Num a => [Int] -> Sparse a s -> Maybe a spartial [] (Sparse a _) = Just a spartial (n:ns) (Sparse _ da) = do-    a' <- lookup n da-    spartial ns a'+  a' <- lookup n da+  spartial ns a' spartial _  Zero         = Nothing {-# INLINE spartial #-} -instance Primal Sparse where-    primal (Sparse a _) = a-    primal Zero = 0+primal :: Num a => Sparse a s -> a+primal (Sparse a _) = a+primal Zero = 0 -instance Lifted Sparse => Mode Sparse where-    auto a = Sparse a IntMap.empty-    zero = Zero-    Zero <**> y    = auto (0 ** primal y)-    _    <**> Zero = auto 1-    x    <**> y@(Sparse b bs)-      | IntMap.null bs = lift1 (**b) (\z -> b *^ z <**> Sparse (b-1) IntMap.empty) x-      | otherwise      = lift2_ (**) (\z xi yi -> (yi *! z /! xi, z *! log1 xi)) x y-    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+(<**>) :: Floating a => Sparse a s -> Sparse a s -> Sparse a s+Zero <**> y    = auto (0 ** primal y)+_    <**> Zero = auto 1+x    <**> y@(Sparse b bs)+  | IntMap.null bs = lift1 (**b) (\z -> b *^ z <**> Sparse (b-1) IntMap.empty) x+  | otherwise      = lift2_ (**) (\z xi yi -> (yi * z / xi, z * log xi)) x y -instance Lifted Sparse => Jacobian Sparse where-    type D Sparse = Sparse-    unary f _ Zero = auto (f 0)-    unary f dadb (Sparse pb bs) = Sparse (f pb) $ mapWithKey (times dadb) bs+instance Num a => Mode (Sparse a s) where+  auto a = Sparse a IntMap.empty+  zero = Zero -    lift1 f _ Zero = auto (f 0)-    lift1 f df b@(Sparse pb bs) = Sparse (f pb) $ mapWithKey (times (df b)) 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 -    lift1_ f _  Zero = auto (f 0)-    lift1_ f df b@(Sparse pb bs) = a where-        a = Sparse (f pb) $ mapWithKey (times (df a b)) bs+infixr 6 <+> -    binary f _    _    Zero           Zero           = auto (f 0 0)-    binary f _    dadc Zero           (Sparse pc dc) = Sparse (f 0  pc) $ mapWithKey (times dadc) dc-    binary f dadb _    (Sparse pb db) Zero           = Sparse (f pb 0 ) $ mapWithKey (times dadb) db-    binary f dadb dadc (Sparse pb db) (Sparse pc dc) = Sparse (f pb pc) $-        unionWith (<+>)-            (mapWithKey (times dadb) db)-            (mapWithKey (times dadc) dc)+(<+>) :: Num a => Sparse a s -> Sparse a s -> Sparse a s+Zero <+> a = a+a <+> Zero = a+Sparse a as <+> Sparse b bs = Sparse (a + b) $ unionWith (<+>) as bs -    lift2 f _  Zero             Zero = auto (f 0 0)-    lift2 f df Zero c@(Sparse pc dc) = Sparse (f 0 pc) $ mapWithKey (times dadc) dc where dadc = snd (df zero c)-    lift2 f df b@(Sparse pb db) Zero = Sparse (f pb 0) $ mapWithKey (times dadb) db where dadb = fst (df b zero)-    lift2 f df b@(Sparse pb db) c@(Sparse pc dc) = Sparse (f pb pc) da where-        (dadb, dadc) = df b c-        da = unionWith (<+>)-            (mapWithKey (times dadb) db)-            (mapWithKey (times dadc) dc)+instance Num a => Jacobian (Sparse a s) where+  type D (Sparse a s) = Sparse a s+  unary f _ Zero = auto (f 0)+  unary f dadb (Sparse pb bs) = Sparse (f pb) $ mapWithKey (times dadb) bs -    lift2_ f _  Zero             Zero = auto (f 0 0)-    lift2_ f df b@(Sparse pb db) Zero = a where a = Sparse (f pb 0) (mapWithKey (times (fst (df a b zero))) db)-    lift2_ f df Zero c@(Sparse pc dc) = a where a = Sparse (f 0 pc) (mapWithKey (times (snd (df a zero c))) dc)-    lift2_ f df b@(Sparse pb db) c@(Sparse pc dc) = a where-        (dadb, dadc) = df a b c-        a = Sparse (f pb pc) da-        da = unionWith (<+>)-            (mapWithKey (times dadb) db)-            (mapWithKey (times dadc) dc)+  lift1 f _ Zero = auto (f 0)+  lift1 f df b@(Sparse pb bs) = Sparse (f pb) $ mapWithKey (times (df b)) bs -deriveLifted id $ conT ''Sparse+  lift1_ f _  Zero = auto (f 0)+  lift1_ f df b@(Sparse pb bs) = a where+    a = Sparse (f pb) $ mapWithKey (times (df a b)) bs +  binary f _    _    Zero           Zero           = auto (f 0 0)+  binary f _    dadc Zero           (Sparse pc dc) = Sparse (f 0  pc) $ mapWithKey (times dadc) dc+  binary f dadb _    (Sparse pb db) Zero           = Sparse (f pb 0 ) $ mapWithKey (times dadb) db+  binary f dadb dadc (Sparse pb db) (Sparse pc dc) = Sparse (f pb pc) $+    unionWith (<+>)+      (mapWithKey (times dadb) db)+      (mapWithKey (times dadc) dc) +  lift2 f _  Zero             Zero = auto (f 0 0)+  lift2 f df Zero c@(Sparse pc dc) = Sparse (f 0 pc) $ mapWithKey (times dadc) dc where dadc = snd (df zero c)+  lift2 f df b@(Sparse pb db) Zero = Sparse (f pb 0) $ mapWithKey (times dadb) db where dadb = fst (df b zero)+  lift2 f df b@(Sparse pb db) c@(Sparse pc dc) = Sparse (f pb pc) da where+    (dadb, dadc) = df b c+    da = unionWith (<+>)+      (mapWithKey (times dadb) db)+      (mapWithKey (times dadc) dc)++  lift2_ f _  Zero             Zero = auto (f 0 0)+  lift2_ f df b@(Sparse pb db) Zero = a where a = Sparse (f pb 0) (mapWithKey (times (fst (df a b zero))) db)+  lift2_ f df Zero c@(Sparse pc dc) = a where a = Sparse (f 0 pc) (mapWithKey (times (snd (df a zero c))) dc)+  lift2_ f df b@(Sparse pb db) c@(Sparse pc dc) = a where+    (dadb, dadc) = df a b c+    a = Sparse (f pb pc) da+    da = unionWith (<+>)+      (mapWithKey (times dadb) db)+      (mapWithKey (times dadc) dc)++#define HEAD Sparse a s+#include "instances.h"+ 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'+  pack :: i -> [Sparse a ()] -> 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 Num a => Grad (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:))+instance Grad i o o' a => Grad (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+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+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+  packs :: i -> [Sparse a ()] -> 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 Num a => Grads (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:))+instance Grads i o a => Grads (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+vgrads i = unpacks (unsafeGrads (packs i)) where+  unsafeGrads f as = ds as $ apply f as {-# INLINE vgrads #-}-
src/Numeric/AD/Internal/Tower.hs view
@@ -1,10 +1,18 @@-{-# LANGUAGE Rank2Types, TypeFamilies, FlexibleContexts, UndecidableInstances, TemplateHaskell, DeriveDataTypeable #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE DeriveDataTypeable #-} {-# OPTIONS_GHC -fno-warn-name-shadowing #-}--- {-# OPTIONS_HADDOCK hide, prune #-}+{-# OPTIONS_HADDOCK not-home #-}+ ----------------------------------------------------------------------------- -- |--- Module      : Numeric.AD.Tower.Internal--- Copyright   : (c) Edward Kmett 2010+-- Copyright   : (c) Edward Kmett 2010-2014 -- License     : BSD3 -- Maintainer  : ekmett@gmail.com -- Stability   : experimental@@ -13,35 +21,39 @@ -----------------------------------------------------------------------------  module Numeric.AD.Internal.Tower-    ( Tower(..)-    , zeroPad-    , zeroPadF-    , transposePadF-    , d-    , d'-    , withD-    , tangents-    , bundle-    , apply-    , getADTower-    , tower-    ) where+  ( Tower(..)+  , zeroPad+  , zeroPadF+  , transposePadF+  , d+  , d'+  , withD+  , tangents+  , bundle+  , apply+  , getADTower+  , tower+  ) where  import Prelude hiding (all) import Control.Applicative hiding ((<**>))+import Control.Monad (join) import Data.Foldable import Data.Data (Data)+import Data.Number.Erf import Data.Typeable (Typeable)-import Language.Haskell.TH-import Numeric.AD.Internal.Types-import Numeric.AD.Internal.Classes+import Numeric.AD.Internal.Combinators+import Numeric.AD.Jacobian+import Numeric.AD.Mode  -- | @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)+newtype Tower a s = Tower { getTower :: [a] } deriving (Data, Typeable) -instance Show a => Show (Tower a) where-    showsPrec n (Tower as) = showParen (n > 10) $ showString "Tower " . showList as+type instance Scalar (Tower a s) = a +instance Show a => Show (Tower a s) where+  showsPrec n (Tower as) = showParen (n > 10) $ showString "Tower " . showList as+ -- Local combinators  zeroPad :: Num a => [a] -> [a]@@ -55,13 +67,13 @@  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+  | 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@@ -74,67 +86,84 @@ d' _        = (0, 0) {-# INLINE d' #-} -tangents :: Tower a -> Tower a+tangents :: Tower a s -> Tower a s tangents (Tower []) = Tower [] tangents (Tower (_:xs)) = Tower xs {-# INLINE tangents #-} -bundle :: a -> Tower a -> Tower a+truncated :: Tower a s -> Bool+truncated (Tower []) = True+truncated _ = False+{-# INLINE truncated #-}++bundle :: a -> Tower a s -> Tower a s bundle a (Tower as) = Tower (a:as) {-# INLINE bundle #-} -withD :: (a, a) -> AD Tower a-withD (a, da) = AD (Tower [a,da])+withD :: (a, a) -> Tower a s+withD (a, da) = Tower [a,da] {-# INLINE withD #-} -apply :: Num a => (AD Tower a -> b) -> a -> b-apply f a = f (AD (Tower [a,1]))+apply :: Num a => (Tower a s -> b) -> a -> b+apply f a = f (Tower [a,1]) {-# INLINE apply #-} -getADTower :: AD Tower a -> [a]-getADTower (AD t) = getTower t+getADTower :: Tower a s -> [a]+getADTower = getTower {-# INLINE getADTower #-} -tower :: [a] -> AD Tower a-tower as = AD (Tower as)+tower :: [a] -> Tower a s+tower = Tower -instance Primal Tower where-    primal (Tower (x:_)) = x-    primal _ = 0+primal :: Num a => Tower a s -> a+primal (Tower (x:_)) = x+primal _ = 0 -instance Lifted Tower => Mode Tower where-    auto a = Tower [a]-    zero = Tower []-    Tower [] <**> y         = auto (0 ** primal y)-    _        <**> Tower []  = auto 1-    x        <**> Tower [y] = lift1 (**y) (\z -> y *^ z <**> Tower [y-1]) x-    x        <**> y         = lift2_ (**) (\z xi yi -> (yi *! z /! xi, z *! log1 xi)) x y+instance Num a => Mode (Tower a s) where+  auto 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) -    a *^ Tower bs = Tower (map (a*) bs)-    Tower as ^* b = Tower (map (*b) as)-    Tower as ^/ b = Tower (map (/b) as)+infixr 6 <+> -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)+(<+>) :: forall a s. Num a => Tower a s -> Tower a s -> Tower a s+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 :: Tower a s) -    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+instance Num a => Jacobian (Tower a s) where+  type D (Tower a s) = Tower a s+  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) -deriveLifted id (conT ''Tower)+  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)) tana where+     (dadb, dadc) = df b c+     tanb = tangents b+     tanc = tangents c+     tana = case (truncated tanb, truncated tanc) of+       (False, False) -> tanb * dadb + tanc * dadc+       (True, False) -> tanc * dadc+       (False, True) -> tanb * dadb+       (True, True) -> zero+  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++(<**>) :: Floating a => Tower a s -> Tower a s -> Tower a s+Tower [] <**> y         = auto (0 ** primal y)+_        <**> Tower []  = auto 1+x        <**> Tower [y] = lift1 (**y) (\z -> y *^ z <**> Tower [y-1]) x+x        <**> y         = lift2_ (**) (\z xi yi -> (yi * z / xi, z * log xi)) x y++#define HEAD Tower a s+#include <instances.h>
− src/Numeric/AD/Internal/Types.hs 
@@ -1,74 +0,0 @@-{-# LANGUAGE CPP #-}-{-# 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--#ifndef MIN_VERSION_base-#define MIN_VERSION_base (x,y,z) 1-#endif--import Data.Data (Data(..), mkDataType, DataType, mkConstr, Constr, constrIndex, Fixity(..))-#if MIN_VERSION_base(4,4,0)-import Data.Typeable (Typeable1(..), Typeable(..), TyCon, mkTyCon3, mkTyConApp, gcast1)-#else-import Data.Typeable (Typeable1(..), Typeable(..), TyCon, mkTyCon, mkTyConApp, gcast1)-#endif-import Language.Haskell.TH-import Numeric.AD.Internal.Classes--{-# ANN module "HLint: ignore Eta reduce" #-}---- | '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-#if MIN_VERSION_base(4,4,0)-adTyCon = mkTyCon3 "ad" "Numeric.AD.Internal.Types" "AD"-#else-adTyCon = mkTyCon "Numeric.AD.Internal.Types.AD"-#endif-{-# 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
− src/Numeric/AD/Internal/Var.hs 
@@ -1,74 +0,0 @@--- {-# OPTIONS_HADDOCK hide, prune #-}--------------------------------------------------------------------------------- |--- Module      :  Numeric.AD.Internal.Var--- Copyright   :  (c) Edward Kmett 2012--- License     :  BSD3--- Maintainer  :  ekmett@gmail.com--- Stability   :  experimental--- Portability :  GHC only------ Variables used for reverse-mode automatic differentiation.--------------------------------------------------------------------------------module Numeric.AD.Internal.Var-    ( Var(..)-    , bind-    , unbind-    , unbindMap-    , unbindWith-    , unbindMapWithDefault-    , Variable(..)-    , vary-    ) where--import Prelude hiding (mapM)-import Data.Array-import Data.IntMap (IntMap, findWithDefault)-import Data.Traversable (Traversable, mapM)-import Numeric.AD.Internal.Types-import Numeric.AD.Internal.Classes---- | 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 f => Var (AD f) where-    var a v = AD (var a v)-    varId (AD v) = varId v---- 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')--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--data Variable a = Variable a {-# UNPACK #-} !Int--instance Var Variable where-  var = Variable-  varId (Variable _ i) = i--instance Primal Variable where-  primal (Variable a _) = a--vary :: Var f => Variable a -> f a-vary (Variable a i) = var a i
+ src/Numeric/AD/Jacobian.hs view
@@ -0,0 +1,39 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE DefaultSignatures #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE PatternGuards #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE UndecidableInstances #-}+-----------------------------------------------------------------------------+-- |+-- Copyright   :  (c) Edward Kmett 2010-2014+-- License     :  BSD3+-- Maintainer  :  ekmett@gmail.com+-- Stability   :  experimental+-- Portability :  GHC only+--+-----------------------------------------------------------------------------++module Numeric.AD.Jacobian+  ( Jacobian(..)+  ) where++import Numeric.AD.Mode++-- | 'Jacobian' is useful for defining new AD primitives in a+-- fairly generic way.+class (Mode t, Mode (D t), Num (D t)) => Jacobian t where+  type D t :: *++  unary  :: (Scalar t -> Scalar t) -> D t -> t -> t+  lift1  :: (Scalar t -> Scalar t) -> (D t -> D t) -> t -> t+  lift1_ :: (Scalar t -> Scalar t) -> (D t -> D t -> D t) -> t -> t++  binary :: (Scalar t -> Scalar t -> Scalar t) -> D t -> D t -> t -> t -> t+  lift2  :: (Scalar t -> Scalar t -> Scalar t) -> (D t -> D t -> (D t, D t)) -> t -> t -> t+  lift2_ :: (Scalar t -> Scalar t -> Scalar t) -> (D t -> D t -> D t -> (D t, D t)) -> t -> t -> t
+ src/Numeric/AD/Jet.hs view
@@ -0,0 +1,101 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE FlexibleContexts #-}+#if __GLASGOW_HASKELL__ >= 707+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE StandaloneDeriving #-}+#endif+-----------------------------------------------------------------------------+-- |+-- Copyright   :  (c) Edward Kmett 2010-2014+-- License     :  BSD3+-- Maintainer  :  ekmett@gmail.com+-- Stability   :  experimental+-- Portability :  GHC only+--+-----------------------------------------------------------------------------+module Numeric.AD.Jet+  ( Jet(..)+  , headJet+  , tailJet+  , jet+  ) where++#ifndef MIN_VERSION_base+#define MIN_VERSION_base(x,y,z) 1+#endif++import Control.Applicative+import Data.Foldable+import Data.Traversable+import Data.Monoid+import Data.Typeable+import Control.Comonad.Cofree++infixl 3 :-++-- | A 'Jet' is a tower of all (higher order) partial derivatives of a function+--+-- At each step, a @'Jet' f@ is wrapped in another layer worth of @f@.+--+-- > a :- f a :- f (f a) :- f (f (f a)) :- ...+data Jet f a = a :- Jet f (f a)++-- | Used to sidestep the need for UndecidableInstances.+newtype Showable = Showable (Int -> String -> String)++instance Show Showable where+  showsPrec d (Showable f) = f d++showable :: Show a => a -> Showable+showable a = Showable (`showsPrec` 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++-- | Take the tail of a 'Jet'.+tailJet :: Jet f a -> Jet f (f a)+tailJet (_ :- as) = as+{-# INLINE tailJet #-}++-- | Take the head of a 'Jet'.+headJet :: Jet f a -> a+headJet (a :- _) = a+{-# INLINE headJet #-}++-- | Construct a 'Jet' by unzipping the layers of a 'Cofree' 'Comonad'.+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)++#if __GLASGOW_HASKELL__ >= 707+deriving instance Typeable Jet+#else+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 MIN_VERSION_base(4,4,0)+jetTyCon = mkTyCon3 "ad" "Numeric.AD.Internal.Jet" "Jet"+#else+jetTyCon = mkTyCon "Numeric.AD.Internal.Jet.Jet"+#endif+{-# NOINLINE jetTyCon #-}+#endif
+ src/Numeric/AD/Mode.hs view
@@ -0,0 +1,61 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE PatternGuards #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE UndecidableInstances #-}+-----------------------------------------------------------------------------+-- |+-- Copyright   :  (c) Edward Kmett 2010-2014+-- License     :  BSD3+-- Maintainer  :  ekmett@gmail.com+-- Stability   :  experimental+-- Portability :  GHC only+--+-----------------------------------------------------------------------------++module Numeric.AD.Mode+  (+  -- * AD modes+    Mode(..)+  , Scalar+  ) where++type family Scalar (t :: *) :: *++infixr 7 *^+infixl 7 ^*+infixr 7 ^/++class (Num t, Num (Scalar 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 -> Bool+  isKnownConstant _ = False++  -- | allowed to return False for zero, but we give more NaN's than strictly necessary then+  isKnownZero :: t -> Bool+  isKnownZero _ = False++  -- | Embed a constant+  auto  :: Scalar t -> t++  -- | Scalar-vector multiplication+  (*^) :: Scalar t -> t -> t+  a *^ b = auto a * b++  -- | Vector-scalar multiplication+  (^*) :: t -> Scalar t -> t+  a ^* b = a * auto b++  -- | Scalar division+  (^/) :: Fractional (Scalar t) => t -> Scalar t -> t+  a ^/ b = a ^* recip b++  -- |+  -- @'zero' = 'lift' 0@+  zero :: t+  zero = auto 0
src/Numeric/AD/Mode/Directed.hs view
@@ -1,8 +1,7 @@-{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE RankNTypes #-} ----------------------------------------------------------------------------- -- |--- Module      :  Numeric.AD.Mode.Directed--- Copyright   :  (c) Edward Kmett 2010-12+-- Copyright   :  (c) Edward Kmett 2010-2014 -- License     :  BSD3 -- Maintainer  :  ekmett@gmail.com -- Stability   :  experimental@@ -13,22 +12,22 @@ -----------------------------------------------------------------------------  module Numeric.AD.Mode.Directed-    (-    -- * Gradients-      grad-    , grad'-    -- * Jacobians-    , jacobian-    , jacobian'-    -- * Derivatives-    , diff-    , diff'-    -- * Exposed Types-    , Direction(..)-    ) where+  (+  -- * Gradients+    grad+  , grad'+  -- * Jacobians+  , jacobian+  , jacobian'+  -- * Derivatives+  , diff+  , diff'+  -- * Exposed Types+  , Direction(..)+  ) where  import Prelude hiding (reverse)-import Numeric.AD.Types+import Numeric.AD.Mode import Data.Traversable (Traversable) import qualified Numeric.AD.Mode.Kahn as K import qualified Numeric.AD.Mode.Forward as F@@ -38,57 +37,57 @@ import Data.Ix  data Direction-    = Forward-    | Kahn-    | Reverse-    | Tower-    | Mixed-    deriving (Show, Eq, Ord, Read, Bounded, Enum, Ix)+  = Forward+  | Kahn+  | 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 Kahn    = K.diff-diff Reverse = R.diff-diff Tower   = T.diff-diff Mixed   = F.diff+diff :: Num a => Direction -> (forall t. Mode t => t -> t) -> a -> a+diff Forward f a = F.diff f a+diff Kahn f a    = K.diff f a+diff Reverse f a = R.diff f a+diff Tower f a   = T.diff f a+diff Mixed f a   = F.diff f a {-# INLINE diff #-} -diff' :: Num a => Direction -> (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a)-diff' Forward = F.diff'-diff' Kahn = K.diff'-diff' Reverse = R.diff'-diff' Tower = T.diff'-diff' Mixed = F.diff'+diff' :: Num a => Direction -> (forall t. Mode t => t -> t) -> a -> (a, a)+diff' Forward f a = F.diff' f a+diff' Kahn f a    = K.diff' f a+diff' Reverse f a = R.diff' f a+diff' Tower f a   = T.diff' f a+diff' Mixed f a   = F.diff' f a {-# INLINE diff' #-} -jacobian :: (Traversable f, Traversable g, Num a) => Direction -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f a)-jacobian Forward = F.jacobian-jacobian Kahn    = K.jacobian-jacobian Reverse = R.jacobian-jacobian Tower   = F.jacobian -- error "jacobian Tower: unimplemented"-jacobian Mixed   = M.jacobian+jacobian :: (Traversable f, Traversable g, Num a) => Direction -> (forall t. Mode t => f (t) -> g (t)) -> f a -> g (f a)+jacobian Forward f a = F.jacobian f a+jacobian Kahn f a    = K.jacobian f a+jacobian Reverse f a = R.jacobian f a+jacobian Tower f a   = F.jacobian f a -- error "jacobian Tower: unimplemented"+jacobian Mixed f a   = M.jacobian f a {-# INLINE jacobian #-} -jacobian' :: (Traversable f, Traversable g, Num a) => Direction -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f a)-jacobian' Forward = F.jacobian'-jacobian' Kahn    = K.jacobian'-jacobian' Reverse = R.jacobian'-jacobian' Tower   = F.jacobian' -- error "jacobian' Tower: unimplemented"-jacobian' Mixed   = M.jacobian'+jacobian' :: (Traversable f, Traversable g, Num a) => Direction -> (forall t. Mode t => f (t) -> g (t)) -> f a -> g (a, f a)+jacobian' Forward f a = F.jacobian' f a+jacobian' Kahn f a    = K.jacobian' f a+jacobian' Reverse f a = R.jacobian' f a+jacobian' Tower f a   = F.jacobian' f a -- error "jacobian' Tower: unimplemented"+jacobian' Mixed f a   = M.jacobian' f a {-# INLINE jacobian' #-} -grad :: (Traversable f, Num a) => Direction -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f a-grad Forward = F.grad-grad Kahn    = K.grad-grad Reverse = R.grad-grad Tower   = F.grad -- error "grad Tower: unimplemented"-grad Mixed   = M.grad+grad :: (Traversable f, Num a) => Direction -> (forall t. Mode t => f (t) -> t) -> f a -> f a+grad Forward f a = F.grad f a+grad Kahn f a    = K.grad f a+grad Reverse f a = R.grad f a+grad Tower f a   = F.grad f a -- error "grad Tower: unimplemented"+grad Mixed f a   = M.grad f a {-# INLINE grad #-} -grad' :: (Traversable f, Num a) => Direction -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f a)-grad' Forward = F.grad'-grad' Kahn    = K.grad'-grad' Reverse = R.grad'-grad' Tower   = F.grad' -- error "grad' Tower: unimplemented"-grad' Mixed   = M.grad'+grad' :: (Traversable f, Num a) => Direction -> (forall t. Mode t => f (t) -> t) -> f a -> (a, f a)+grad' Forward f a = F.grad' f a+grad' Kahn f a    = K.grad' f a+grad' Reverse f a = R.grad' f a+grad' Tower f a   = F.grad' f a -- error "grad' Tower: unimplemented"+grad' Mixed f a   = M.grad' f a {-# INLINE grad' #-}
src/Numeric/AD/Mode/Forward.hs view
@@ -1,8 +1,8 @@+{-# LANGUAGE CPP #-} {-# LANGUAGE Rank2Types #-} ----------------------------------------------------------------------------- -- |--- Module      :  Numeric.AD.Mode.Forward--- Copyright   :  (c) Edward Kmett 2010+-- Copyright   :  (c) Edward Kmett 2010-2014 -- License     :  BSD3 -- Maintainer  :  ekmett@gmail.com -- Stability   :  experimental@@ -13,59 +13,57 @@ -----------------------------------------------------------------------------  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+  ( 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+import Numeric.AD.Internal.On  -- | Compute the directional derivative of a function given a zipped up 'Functor' of the input values and their derivatives-du :: (Functor f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f (a, a) -> a+du :: (Functor f, Num a) => (forall s. f (Forward a s) -> Forward a s) -> f (a, a) -> a du f = tangent . f . fmap (uncurry bundle) {-# INLINE du #-}  -- | Compute the answer and directional derivative of a function given a zipped up 'Functor' of the input values and their derivatives-du' :: (Functor f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f (a, a) -> (a, a)+du' :: (Functor f, Num a) => (forall s. f (Forward a s) -> Forward a s) -> f (a, a) -> (a, a) du' f = unbundle . f . fmap (uncurry bundle) {-# INLINE du' #-}  -- | Compute a vector of directional derivatives for a function given a zipped up 'Functor' of the input values and their derivatives.-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 :: (Functor f, Functor g, Num a) => (forall s. f (Forward a s) -> g (Forward a s)) -> f (a, a) -> g a duF f = fmap tangent . f . fmap (uncurry bundle) {-# INLINE duF #-}  -- | Compute a vector of answers and directional derivatives for a function given a zipped up 'Functor' of the input values and their derivatives.-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' :: (Functor f, Functor g, Num a) => (forall s. f (Forward a s) -> g (Forward a s)) -> f (a, a) -> g (a, a) duF' f = fmap unbundle . f . fmap (uncurry bundle) {-# INLINE duF' #-} @@ -73,7 +71,7 @@ -- -- >>> diff sin 0 -- 1.0-diff :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> a+diff :: Num a => (forall s. Forward a s -> Forward a s) -> a -> a diff f a = tangent $ apply f a {-# INLINE diff #-} @@ -90,7 +88,7 @@ -- >>> diff' exp 0 -- (1.0,1.0) -diff' :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a)+diff' :: Num a => (forall s. Forward a s -> Forward a s) -> a -> (a, a) diff' f a = unbundle $ apply f a {-# INLINE diff' #-} @@ -98,7 +96,7 @@ -- -- >>> diffF (\a -> [sin a, cos a]) 0 -- [1.0,-0.0]-diffF :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f a+diffF :: (Functor f, Num a) => (forall s. Forward a s -> f (Forward a s)) -> a -> f a diffF f a = tangent <$> apply f a {-# INLINE diffF #-} @@ -106,78 +104,75 @@ -- -- >>> diffF' (\a -> [sin a, cos a]) 0 -- [(0.0,1.0),(1.0,-0.0)]-diffF' :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f (a, a)+diffF' :: (Functor f, Num a) => (forall s. Forward a s -> f (Forward a s)) -> 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 :: (Traversable f, Functor g, Num a) => (forall s. f (Forward a s) -> g (Forward a s)) -> f a -> f (g a) jacobianT f = bind (fmap tangent . f) {-# INLINE jacobianT #-}  -- | A fast, simple, transposed Jacobian computed with 'Forward' mode 'AD' that combines the output with the input.-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+jacobianWithT :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. f (Forward a s) -> g (Forward a s)) -> f a -> f (g b)+jacobianWithT g f = bindWith g' f where+  g' a ga = g a . tangent <$> ga {-# INLINE jacobianWithT #-}+#ifdef HLINT {-# ANN jacobianWithT "HLint: ignore Eta reduce" #-}+#endif  -- | Compute the Jacobian using 'Forward' mode 'AD'. This must transpose the result, so 'jacobianT' is faster and allows more result types. -- -- -- >>> jacobian (\[x,y] -> [y,x,x+y,x*y,exp x * sin y]) [pi,1] -- [[0.0,1.0],[1.0,0.0],[1.0,1.0],[1.0,3.141592653589793],[19.472221418841606,12.502969588876512]]-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+jacobian :: (Traversable f, Traversable g, Num a) => (forall s . f (Forward a s) -> g (Forward a s)) -> f a -> g (f a)+jacobian f as = transposeWith (const id) t p where+  (p, t) = bind' (fmap tangent . f) as {-# INLINE jacobian #-}  -- | Compute the Jacobian using 'Forward' mode 'AD' and combine the output with the input. This must transpose the result, so 'jacobianWithT' is faster, and allows more result types.-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+jacobianWith :: (Traversable f, Traversable g, Num a) => (a -> a -> b) -> (forall s. f (Forward a s) -> g (Forward a s)) -> 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 #-}  -- | Compute the Jacobian using 'Forward' mode 'AD' along with the actual answer.-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')+jacobian' :: (Traversable f, Traversable g, Num a) => (forall s. f (Forward a s) -> g (Forward a s)) -> 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' #-}  -- | Compute the Jacobian using 'Forward' mode 'AD' combined with the input using a user specified function, along with the actual answer.-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+jacobianWith' :: (Traversable f, Traversable g, Num a) => (a -> a -> b) -> (forall s. f (Forward a s) -> g (Forward a s)) -> 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' #-}  -- | Compute the gradient of a function using forward mode AD. -- -- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Wengert.grad' for @n@ inputs, in exchange for better space utilization.-grad :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f a+grad :: (Traversable f, Num a) => (forall s. f (Forward a s) -> Forward a s) -> f a -> f a grad f = bind (tangent . f) {-# INLINE grad #-}  -- | Compute the gradient and answer to a function using forward mode AD. -- -- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Wengert.grad'' for @n@ inputs, in exchange for better space utilization.-grad' :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f a)-grad' f as = (primal b, tangent <$> bs)-    where-        (b, bs) = bind' f as+grad' :: (Traversable f, Num a) => (forall s. f (Forward a s) -> Forward a s) -> f a -> (a, f a)+grad' f as = (primal b, tangent <$> bs) where+  (b, bs) = bind' f as {-# INLINE grad' #-}  -- | Compute the gradient of a function using forward mode AD and combine the result with the input using a user-specified function. -- -- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Wengert.gradWith' for @n@ inputs, in exchange for better space utilization.-gradWith :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f b+gradWith :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. f (Forward a s) -> Forward a s) -> f a -> f b gradWith g f = bindWith g (tangent . f) {-# INLINE gradWith #-} @@ -188,20 +183,17 @@ -- -- >>> gradWith' (,) sum [0..4] -- (10,[(0,1),(1,1),(2,1),(3,1),(4,1)])-gradWith' :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f b)-gradWith' g f as = (primal $ f (AD . Lift <$> as), bindWith g (tangent . f) as)+gradWith' :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. f (Forward a s) -> Forward a s) -> f a -> (a, f b)+gradWith' g f as = (primal $ f (Lift <$> as), bindWith g (tangent . f) as) {-# 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+hessianProduct :: (Traversable f, Num a) => (forall s s'. f (On (Forward (Forward a s') s)) -> On (Forward (Forward a s') s)) -> f (a, a) -> f a+hessianProduct f = duF $ grad $ off . f . fmap On+{-# INLINE hessianProduct #-}  -- | 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)+hessianProduct' :: (Traversable f, Num a) => (forall s s'. f (On (Forward (Forward a s') s)) -> On (Forward (Forward a s') s)) -> f (a, a) -> f (a, a)+hessianProduct' f = duF' $ grad $ off . f . fmap On+{-# INLINE hessianProduct' #-}
+ src/Numeric/AD/Mode/Forward/Double.hs view
@@ -0,0 +1,170 @@+{-# LANGUAGE RankNTypes #-}+module Numeric.AD.Mode.Forward.Double+  ( ForwardDouble+  -- * Gradient+  , grad+  , grad'+  , gradWith+  , gradWith'+  -- * Jacobian+  , jacobian+  , jacobian'+  , jacobianWith+  , jacobianWith'+  -- * Transposed Jacobian+  , jacobianT+  , jacobianWithT+  -- * Derivatives+  , diff+  , diff'+  , diffF+  , diffF'+  -- * Directional Derivatives+  , du+  , du'+  , duF+  , duF'+  ) where++import Control.Applicative+import Data.Traversable (Traversable)+import Numeric.AD.Mode+import Numeric.AD.Internal.Forward.Double++-- | Compute the directional derivative of a function given a zipped up 'Functor' of the input values and their derivatives+du :: Functor f => (forall s. f (ForwardDouble s) -> ForwardDouble s) -> f (Double, Double) -> Double+du f = tangent . f . fmap (uncurry bundle)+{-# INLINE du #-}++-- | Compute the answer and directional derivative of a function given a zipped up 'Functor' of the input values and their derivatives+du' :: Functor f => (forall s. f (ForwardDouble s) -> ForwardDouble s) -> f (Double, Double) -> (Double, Double)+du' f = unbundle . f . fmap (uncurry bundle)+{-# INLINE du' #-}++-- | Compute a vector of directional derivatives for a function given a zipped up 'Functor' of the input values and their derivatives.+duF :: (Functor f, Functor g) => (forall s. f (ForwardDouble s) -> g (ForwardDouble s)) -> f (Double, Double) -> g Double+duF f = fmap tangent . f . fmap (uncurry bundle)+{-# INLINE duF #-}++-- | Compute a vector of answers and directional derivatives for a function given a zipped up 'Functor' of the input values and their derivatives.+duF' :: (Functor f, Functor g) => (forall s. f (ForwardDouble s) -> g (ForwardDouble s)) -> f (Double, Double) -> g (Double, Double)+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 0+-- 1.0+diff :: (forall s. ForwardDouble s -> ForwardDouble s) -> Double -> Double+diff f a = tangent $ apply f a+{-# INLINE diff #-}++-- | The 'diff'' function calculates the result and first derivative of scalar-to-scalar function by 'Forward' mode 'AD'+--+-- @+-- 'diff'' 'sin' == 'sin' 'Control.Arrow.&&&' 'cos'+-- 'diff'' f = f 'Control.Arrow.&&&' d f+-- @+--+-- >>> diff' sin 0+-- (0.0,1.0)+--+-- >>> diff' exp 0+-- (1.0,1.0)+diff' :: (forall s. ForwardDouble s -> ForwardDouble s) -> Double -> (Double, Double)+diff' f a = unbundle $ apply f a+{-# INLINE diff' #-}++-- | The 'diffF' function calculates the first derivatives of scalar-to-nonscalar function by 'Forward' mode 'AD'+--+-- >>> diffF (\a -> [sin a, cos a]) 0+-- [1.0,-0.0]+diffF :: Functor f => (forall s. ForwardDouble s -> f (ForwardDouble s)) -> Double -> f Double+diffF f a = tangent <$> apply f a+{-# INLINE diffF #-}++-- | The 'diffF'' function calculates the result and first derivatives of a scalar-to-non-scalar function by 'Forward' mode 'AD'+--+-- >>> diffF' (\a -> [sin a, cos a]) 0+-- [(0.0,1.0),(1.0,-0.0)]+diffF' :: Functor f => (forall s. ForwardDouble s -> f (ForwardDouble s)) -> Double -> f (Double, Double)+diffF' f a = unbundle <$> apply f a+{-# INLINE diffF' #-}++-- | A fast, simple, transposed Jacobian computed with forward-mode AD.+jacobianT :: (Traversable f, Functor g) => (forall s. f (ForwardDouble s) -> g (ForwardDouble s)) -> f Double -> f (g Double)+jacobianT f = bind (fmap tangent . f)+{-# INLINE jacobianT #-}++-- | A fast, simple, transposed Jacobian computed with 'Forward' mode 'AD' that combines the output with the input.+jacobianWithT :: (Traversable f, Functor g) => (Double -> Double -> b) -> (forall s. f (ForwardDouble s) -> g (ForwardDouble s)) -> f Double -> f (g b)+jacobianWithT g f = bindWith g' f where+  g' a ga = g a . tangent <$> ga+{-# INLINE jacobianWithT #-}+{-# ANN jacobianWithT "HLint: ignore Eta reduce" #-}++-- | Compute the Jacobian using 'Forward' mode 'AD'. This must transpose the result, so 'jacobianT' is faster and allows more result types.+--+--+-- >>> jacobian (\[x,y] -> [y,x,x+y,x*y,exp x * sin y]) [pi,1]+-- [[0.0,1.0],[1.0,0.0],[1.0,1.0],[1.0,3.141592653589793],[19.472221418841606,12.502969588876512]]+jacobian :: (Traversable f, Traversable g) => (forall s . f (ForwardDouble s) -> g (ForwardDouble s)) -> f Double -> g (f Double)+jacobian f as = transposeWith (const id) t p where+  (p, t) = bind' (fmap tangent . f) as+{-# INLINE jacobian #-}++-- | Compute the Jacobian using 'Forward' mode 'AD' and combine the output with the input. This must transpose the result, so 'jacobianWithT' is faster, and allows more result types.+jacobianWith :: (Traversable f, Traversable g) => (Double -> Double -> b) -> (forall s. f (ForwardDouble s) -> g (ForwardDouble s)) -> f Double -> 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 #-}++-- | Compute the Jacobian using 'Forward' mode 'AD' along with the actual answer.+jacobian' :: (Traversable f, Traversable g) => (forall s. f (ForwardDouble s) -> g (ForwardDouble s)) -> f Double -> g (Double, f Double)+jacobian' f as = transposeWith row t p where+  (p, t) = bind' f as+  row x as' = (primal x, tangent <$> as')+{-# INLINE jacobian' #-}++-- | Compute the Jacobian using 'Forward' mode 'AD' combined with the input using a user specified function, along with the actual answer.+jacobianWith' :: (Traversable f, Traversable g) => (Double -> Double -> b) -> (forall s. f (ForwardDouble s) -> g (ForwardDouble s)) -> f Double -> g (Double, 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' #-}++-- | Compute the gradient of a function using forward mode AD.+--+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Wengert.grad' for @n@ inputs, in exchange for better space utilization.+grad :: (Traversable f) => (forall s. f (ForwardDouble s) -> ForwardDouble s) -> f Double -> f Double+grad f = bind (tangent . f)+{-# INLINE grad #-}++-- | Compute the gradient and answer to a function using forward mode AD.+--+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Wengert.grad'' for @n@ inputs, in exchange for better space utilization.+grad' :: (Traversable f) => (forall s. f (ForwardDouble s) -> ForwardDouble s) -> f Double -> (Double, f Double)+grad' f as = (primal b, tangent <$> bs)+    where+        (b, bs) = bind' f as+{-# INLINE grad' #-}++-- | Compute the gradient of a function using forward mode AD and combine the result with the input using a user-specified function.+--+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Wengert.gradWith' for @n@ inputs, in exchange for better space utilization.+gradWith :: (Traversable f) => (Double -> Double -> b) -> (forall s. f (ForwardDouble s) -> ForwardDouble s) -> f Double -> f b+gradWith g f = bindWith g (tangent . f)+{-# INLINE gradWith #-}++-- | Compute the gradient of a function using forward mode AD and the answer, and combine the result with the input using a+-- user-specified function.+--+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Wengert.gradWith'' for @n@ inputs, in exchange for better space utilization.+--+-- >>> gradWith' (,) sum [0..4]+-- (10.0,[(0.0,1.0),(1.0,1.0),(2.0,1.0),(3.0,1.0),(4.0,1.0)])+gradWith' :: (Traversable f) => (Double -> Double -> b) -> (forall s. f (ForwardDouble s) -> ForwardDouble s) -> f Double -> (Double, f b)+gradWith' g f as = (primal $ f (auto <$> as), bindWith g (tangent . f) as)+{-# INLINE gradWith' #-}
src/Numeric/AD/Mode/Kahn.hs view
@@ -1,8 +1,13 @@-{-# LANGUAGE Rank2Types, TemplateHaskell, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, UndecidableInstances, ScopedTypeVariables #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TemplateHaskell #-}+{-# LANGUAGE UndecidableInstances #-} ----------------------------------------------------------------------------- -- |--- Module      :  Numeric.AD.Mode.Kahn--- Copyright   :  (c) Edward Kmett 2010+-- Copyright   :  (c) Edward Kmett 2010-2014 -- License     :  BSD3 -- Maintainer  :  ekmett@gmail.com -- Stability   :  experimental@@ -19,61 +24,57 @@ -----------------------------------------------------------------------------  module Numeric.AD.Mode.Kahn-    (-    -- * Gradient-      grad-    , grad'-    , gradWith-    , gradWith'--    -- * Jacobian-    , jacobian-    , jacobian'-    , jacobianWith-    , jacobianWith'-    -- * Hessian-    , hessian-    , hessianF-    -- * Derivatives-    , diff-    , diff'-    , diffF-    , diffF'-    -- * Unsafe Variadic Gradient-    , vgrad, vgrad'-    , Grad-    ) where+  ( Kahn+  -- * Gradient+  , grad+  , grad'+  , gradWith+  , gradWith'+  -- * Jacobian+  , jacobian+  , jacobian'+  , jacobianWith+  , jacobianWith'+  -- * Hessian+  , hessian+  , hessianF+  -- * Derivatives+  , diff+  , diff'+  , diffF+  , diffF'+  -- * Unsafe Variadic Gradient+  -- $vgrad+  , vgrad, vgrad'+  , Grad+  ) where  import Control.Applicative ((<$>))+import Data.Functor.Compose import Data.Traversable (Traversable)--import Numeric.AD.Types-import Numeric.AD.Internal.Classes-import Numeric.AD.Internal.Composition+import Numeric.AD.Internal.On import Numeric.AD.Internal.Kahn-import Numeric.AD.Internal.Var --- | The 'grad' function calculates the gradient of a non-scalar-to-scalar function with reverse-mode AD in a single pass.+-- | The 'grad' function calculates the gradient of a non-scalar-to-scalar function with kahn-mode AD in a single pass. -- -- >>> grad (\[x,y,z] -> x*y+z) [1,2,3] -- [2,1,1]--grad :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f a-grad f as = unbind vs (partialArray bds $ f vs)-    where (vs,bds) = bind as+grad :: (Traversable f, Num a) => (forall s. f (Kahn a s) -> Kahn a s) -> f a -> f a+grad f as = unbind vs (partialArray bds $ f vs) where+  (vs,bds) = bind as {-# INLINE grad #-} --- | The 'grad'' function calculates the result and gradient of a non-scalar-to-scalar function with reverse-mode AD in a single pass.+-- | The 'grad'' function calculates the result and gradient of a non-scalar-to-scalar function with kahn-mode AD in a single pass. -- -- >>> grad' (\[x,y,z] -> 4*x*exp y+cos z) [1,2,3] -- (28.566231899122155,[29.5562243957226,29.5562243957226,-0.1411200080598672])-grad' :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f a)-grad' f as = (primal r, unbind vs $ partialArray bds r)-    where (vs, bds) = bind as-          r = f vs+grad' :: (Traversable f, Num a) => (forall s. f (Kahn a s) -> Kahn a s) -> 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.+-- | @'grad' g f@ function calculates the gradient of a non-scalar-to-scalar function @f@ with kahn-mode AD in a single pass. -- The gradient is combined element-wise with the argument using the function @g@. -- -- @@@ -82,46 +83,46 @@ -- @ -- ---gradWith :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f b-gradWith g f as = unbindWith g vs (partialArray bds $ f vs)-    where (vs,bds) = bind as+gradWith :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. f (Kahn a s) -> Kahn a s) -> f a -> f b+gradWith g f as = unbindWith g vs (partialArray bds $ f vs) where+  (vs,bds) = bind as {-# INLINE gradWith #-} --- | @'grad'' g f@ calculates the result and gradient of a non-scalar-to-scalar function @f@ with reverse-mode AD in a single pass+-- | @'grad'' g f@ calculates the result and gradient of a non-scalar-to-scalar function @f@ with kahn-mode AD in a single pass -- the gradient is combined element-wise with the argument using the function @g@. -- -- @'grad'' == 'gradWith'' (\_ dx -> dx)@-gradWith' :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f b)-gradWith' g f as = (primal r, unbindWith g vs $ partialArray bds r)-    where (vs, bds) = bind as-          r = f vs+gradWith' :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. f (Kahn a s) -> Kahn a s) -> 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.+-- | The 'jacobian' function calculates the jacobian of a non-scalar-to-non-scalar function with kahn AD lazily in @m@ passes for @m@ outputs. -- -- >>> jacobian (\[x,y] -> [y,x,x*y]) [2,1] -- [[0,1],[1,0],[1,2]] -- -- >>> jacobian (\[x,y] -> [exp y,cos x,x+y]) [1,2] -- [[0.0,7.38905609893065],[-0.8414709848078965,0.0],[1.0,1.0]]-jacobian :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f a)+jacobian :: (Traversable f, Functor g, Num a) => (forall s. f (Kahn a s) -> g (Kahn a s)) -> f a -> g (f a) jacobian f as = unbind vs . partialArray bds <$> f vs where-    (vs, bds) = bind as+  (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,+-- | The 'jacobian'' function calculates both the result and the Jacobian of a nonscalar-to-nonscalar function, using @m@ invocations of kahn AD, -- where @m@ is the output dimensionality. Applying @fmap snd@ to the result will recover the result of 'jacobian' -- | An alias for 'gradF'' -- -- ghci> jacobian' (\[x,y] -> [y,x,x*y]) [2,1] -- [(1,[0,1]),(2,[1,0]),(2,[1,2])]-jacobian' :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f a)+jacobian' :: (Traversable f, Functor g, Num a) => (forall s. f (Kahn a s) -> g (Kahn a s)) -> 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))+  (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.+-- | 'jacobianWith g f' calculates the Jacobian of a non-scalar-to-non-scalar function @f@ with kahn 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@. --@@ -129,21 +130,21 @@ -- '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 :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. f (Kahn a s) -> g (Kahn a s)) -> f a -> g (f b) jacobianWith g f as = unbindWith g vs . partialArray bds <$> f vs where-    (vs, bds) = bind as+  (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,+-- | 'jacobianWith' g f' calculates both the result and the Jacobian of a nonscalar-to-nonscalar function @f@, using @m@ invocations of kahn 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' :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. f (Kahn a s) -> g (Kahn a s)) -> 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))+  (vs, bds) = bind as+  row a = (primal a, unbindWith g vs (partialArray bds a)) {-# INLINE jacobianWith' #-}  -- | Compute the derivative of a function.@@ -153,7 +154,7 @@ -- -- >>> cos 0 -- 1.0-diff :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> a+diff :: Num a => (forall s. Kahn a s -> Kahn a s) -> a -> a diff f a = derivative $ f (var a 0) {-# INLINE diff #-} @@ -163,7 +164,7 @@ -- -- >>> diff' sin 0 -- (0.0,1.0)-diff' :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a)+diff' :: Num a => (forall s. Kahn a s -> Kahn a s) -> a -> (a, a) diff' f a = derivative' $ f (var a 0) {-# INLINE diff' #-} @@ -171,7 +172,7 @@ -- -- >>> diffF (\a -> [sin a, cos a]) 0 -- [1.0,0.0]-diffF :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f a+diffF :: (Functor f, Num a) => (forall s. Kahn a s -> f (Kahn a s)) -> a -> f a diffF f a = derivative <$> f (var a 0) {-# INLINE diffF #-} @@ -180,25 +181,35 @@ -- -- >>> diffF' (\a -> [sin a, cos a]) 0 -- [(0.0,1.0),(1.0,0.0)]-diffF' :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f (a, a)+diffF' :: (Functor f, Num a) => (forall s. Kahn a s -> f (Kahn a s)) -> 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.++-- | Compute the 'hessian' via the 'jacobian' of the gradient. gradient is computed in 'Kahn' mode and then the 'jacobian' is computed in 'Kahn' mode. -- -- However, since the @'grad' f :: f a -> f a@ is square this is not as fast as using the forward-mode 'jacobian' of a reverse mode gradient provided by 'Numeric.AD.hessian'. -- -- >>> hessian (\[x,y] -> x*y) [1,2] -- [[0,1],[1,0]]-hessian :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f (f a)-hessian f = jacobian (grad (decomposeMode . f . fmap composeMode))+hessian :: (Traversable f, Num a) => (forall s s'. f (On (Kahn (Kahn a s') s)) -> (On (Kahn (Kahn a s') s))) -> f a -> f (f a)+hessian f = jacobian (grad (off . f . fmap On)) --- | 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.+-- | Compute the order 3 Hessian tensor on a non-scalar-to-non-scalar function via the 'Kahn'-mode Jacobian of the 'Kahn'-mode Jacobian of the function. -- -- Less efficient than 'Numeric.AD.Mode.Mixed.hessianF'. -- -- >>> hessianF (\[x,y] -> [x*y,x+y,exp x*cos y]) [1,2] -- [[[0.0,1.0],[1.0,0.0]],[[0.0,0.0],[0.0,0.0]],[[-1.1312043837568135,-2.4717266720048188],[-2.4717266720048188,1.1312043837568135]]]-hessianF :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f (f a))-hessianF f = decomposeFunctor . jacobian (ComposeFunctor . jacobian (fmap decomposeMode . f . fmap composeMode))+hessianF :: (Traversable f, Functor g, Num a) => (forall s s'. f (On (Kahn (Kahn a s') s)) -> g (On (Kahn (Kahn a s') s))) -> f a -> g (f (f a))+hessianF f = getCompose . jacobian (Compose . jacobian (fmap off . f . fmap On)) ++-- $vgrad+--+-- 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 'auto' calls you use when taking the gradient+-- of a function that takes gradients!
src/Numeric/AD/Mode/Reverse.hs view
@@ -1,8 +1,14 @@-{-# LANGUAGE Rank2Types, TemplateHaskell, BangPatterns, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, UndecidableInstances, ScopedTypeVariables #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE TemplateHaskell #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE ScopedTypeVariables #-} ----------------------------------------------------------------------------- -- |--- Module      :  Numeric.AD.Mode.Reverse--- Copyright   :  (c) Edward Kmett 2010+-- Copyright   :  (c) Edward Kmett 2010-2014 -- License     :  BSD3 -- Maintainer  :  ekmett@gmail.com -- Stability   :  experimental@@ -14,56 +20,54 @@ -----------------------------------------------------------------------------  module Numeric.AD.Mode.Reverse-    (-    -- * Gradient-      grad-    , grad'-    , gradWith-    , gradWith'+  ( Reverse+  -- * Gradient+  , grad+  , grad'+  , gradWith+  , gradWith' -    -- * Jacobian-    , jacobian-    , jacobian'-    , jacobianWith-    , jacobianWith'+  -- * Jacobian+  , jacobian+  , jacobian'+  , jacobianWith+  , jacobianWith' -    -- * Hessian-    , hessian-    , hessianF+  -- * Hessian+  , hessian+  , hessianF -    -- * Derivatives-    , diff-    , diff'-    , diffF-    , diffF'-    ) where+  -- * Derivatives+  , diff+  , diff'+  , diffF+  , diffF'+  ) where  import Control.Applicative ((<$>))+import Data.Functor.Compose+import Data.Reflection (Reifies) import Data.Traversable (Traversable)--import Numeric.AD.Types-import Numeric.AD.Internal.Classes-import Numeric.AD.Internal.Composition+import Numeric.AD.Internal.On import Numeric.AD.Internal.Reverse-import Numeric.AD.Internal.Var  -- | The 'grad' function calculates the gradient of a non-scalar-to-scalar function with reverse-mode AD in a single pass. -- -- -- >>> grad (\[x,y,z] -> x*y+z) [1,2,3] -- [2,1,1]-grad :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f a-grad f as = reifyTape (snd bds) $ \p -> unbind vs $! partialArrayOf p bds $! f $ vary <$> vs-  where (vs, bds) = bind as+grad :: (Traversable f, Num a) => (forall s. Reifies s Tape => f (Reverse a s) -> Reverse a s) -> f a -> f a+grad f as = reifyTape (snd bds) $ \p -> unbind vs $! partialArrayOf p bds $! f vs where+  (vs, bds) = bind as {-# INLINE grad #-}  -- | The 'grad'' function calculates the result and gradient of a non-scalar-to-scalar function with reverse-mode AD in a single pass. -- -- >>> grad' (\[x,y,z] -> x*y+z) [1,2,3] -- (5,[2,1,1])-grad' :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f a)-grad' f as = reifyTape (snd bds) $ \p ->-  let r = f (fmap vary vs) in (primal r, unbind vs $! partialArrayOf p bds $! r)+grad' :: (Traversable f, Num a) => (forall s. Reifies s Tape => f (Reverse a s) -> Reverse a s) -> f a -> (a, f a)+grad' f as = reifyTape (snd bds) $ \p -> case f vs of+   r -> (primal r, unbind vs $! partialArrayOf p bds $! r)   where (vs, bds) = bind as {-# INLINE grad' #-} @@ -74,8 +78,8 @@ -- 'grad' == 'gradWith' (\_ dx -> dx) -- 'id' == 'gradWith' 'const' -- @-gradWith :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f b-gradWith g f as = reifyTape (snd bds) $ \p -> unbindWith g vs $! partialArrayOf p bds $! f $ vary <$> vs+gradWith :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Reifies s Tape => f (Reverse a s) -> Reverse a s) -> f a -> f b+gradWith g f as = reifyTape (snd bds) $ \p -> unbindWith g vs $! partialArrayOf p bds $! f vs   where (vs,bds) = bind as {-# INLINE gradWith #-} @@ -85,19 +89,19 @@ -- @ -- 'grad'' == 'gradWith'' (\_ dx -> dx) -- @-gradWith' :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f b)-gradWith' g f as = reifyTape (snd bds) $ \p ->-   let r = f (fmap vary vs) in (primal r, unbindWith g vs $! partialArrayOf p bds $! r)-    where (vs, bds) = bind as+gradWith' :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Reifies s Tape => f (Reverse a s) -> Reverse a s) -> f a -> (a, f b)+gradWith' g f as = reifyTape (snd bds) $ \p -> case f vs of+   r -> (primal r, unbindWith g vs $! partialArrayOf p bds $! r)+  where (vs, bds) = bind as {-# INLINE gradWith' #-}  -- | The 'jacobian' function calculates the jacobian of a non-scalar-to-non-scalar function with reverse AD lazily in @m@ passes for @m@ outputs. -- -- >>> jacobian (\[x,y] -> [y,x,x*y]) [2,1] -- [[0,1],[1,0],[1,2]]-jacobian :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f a)-jacobian f as = reifyTape (snd bds) $ \p -> unbind vs . partialArrayOf p bds <$> f (fmap vary vs)-  where (vs, bds) = bind as+jacobian :: (Traversable f, Functor g, Num a) => (forall s. Reifies s Tape => f (Reverse a s) -> g (Reverse a s)) -> f a -> g (f a)+jacobian f as = reifyTape (snd bds) $ \p -> unbind vs . partialArrayOf p 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,@@ -106,10 +110,10 @@ -- -- >>> jacobian' (\[x,y] -> [y,x,x*y]) [2,1] -- [(1,[0,1]),(2,[1,0]),(2,[1,2])]-jacobian' :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f a)+jacobian' :: (Traversable f, Functor g, Num a) => (forall s. Reifies s Tape => f (Reverse a s) -> g (Reverse a s)) -> f a -> g (a, f a) jacobian' f as = reifyTape (snd bds) $ \p ->   let row a = (primal a, unbind vs $! partialArrayOf p bds $! a)-  in row <$> f (vary <$> vs)+  in row <$> f vs   where (vs, bds) = bind as {-# INLINE jacobian' #-} @@ -121,9 +125,9 @@ -- 'jacobian' == 'jacobianWith' (\_ dx -> dx) -- 'jacobianWith' 'const' == (\f x -> 'const' x '<$>' f x) -- @-jacobianWith :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f b)-jacobianWith g f as = reifyTape (snd bds) $ \p -> unbindWith g vs . partialArrayOf p bds <$> f (fmap vary vs) where-    (vs, bds) = bind as+jacobianWith :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. Reifies s Tape => f (Reverse a s) -> g (Reverse a s)) -> f a -> g (f b)+jacobianWith g f as = reifyTape (snd bds) $ \p -> unbindWith g vs . partialArrayOf p 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,@@ -133,10 +137,10 @@ -- -- @'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' :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. Reifies s Tape => f (Reverse a s) -> g (Reverse a s)) -> f a -> g (a, f b) jacobianWith' g f as = reifyTape (snd bds) $ \p ->   let row a = (primal a, unbindWith g vs $! partialArrayOf p bds $! a)-  in row <$> f (vary <$> vs)+  in row <$> f vs   where (vs, bds) = bind as {-# INLINE jacobianWith' #-} @@ -144,7 +148,7 @@ -- -- >>> diff sin 0 -- 1.0-diff :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> a+diff :: Num a => (forall s. Reifies s Tape => Reverse a s -> Reverse a s) -> a -> a diff f a = reifyTape 1 $ \p -> derivativeOf p $! f (var a 0) {-# INLINE diff #-} @@ -155,7 +159,7 @@ -- -- >>> diff' exp 0 -- (1.0,1.0)-diff' :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a)+diff' :: Num a => (forall s. Reifies s Tape => Reverse a s -> Reverse a s) -> a -> (a, a) diff' f a = reifyTape 1 $ \p -> derivativeOf' p $! f (var a 0) {-# INLINE diff' #-} @@ -164,7 +168,7 @@ -- >>> diffF (\a -> [sin a, cos a]) 0 -- [1.0,0.0] ---diffF :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f a+diffF :: (Functor f, Num a) => (forall s. Reifies s Tape => Reverse a s -> f (Reverse a s)) -> a -> f a diffF f a = reifyTape 1 $ \p -> derivativeOf p <$> f (var a 0) {-# INLINE diffF #-} @@ -172,7 +176,7 @@ -- -- >>> diffF' (\a -> [sin a, cos a]) 0 -- [(0.0,1.0),(1.0,0.0)]-diffF' :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f (a, a)+diffF' :: (Functor f, Num a) => (forall s. Reifies s Tape => Reverse a s -> f (Reverse a s)) -> a -> f (a, a) diffF' f a = reifyTape 1 $ \p -> derivativeOf' p <$> f (var a 0) {-# INLINE diffF' #-} @@ -182,8 +186,9 @@ -- -- >>> hessian (\[x,y] -> x*y) [1,2] -- [[0,1],[1,0]]-hessian :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f (f a)-hessian f = jacobian (grad (decomposeMode . f . fmap composeMode))+hessian :: (Traversable f, Num a) => (forall s s'. (Reifies s Tape, Reifies s' Tape) => f (On (Reverse (Reverse a s') s)) -> (On (Reverse (Reverse a s') s))) -> f a -> f (f a)+hessian f = jacobian (grad (off . f . fmap On))+{-# INLINE hessian #-}  -- | 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. --@@ -191,5 +196,6 @@ -- -- >>> hessianF (\[x,y] -> [x*y,x+y,exp x*cos y]) [1,2] -- [[[0.0,1.0],[1.0,0.0]],[[0.0,0.0],[0.0,0.0]],[[-1.1312043837568135,-2.4717266720048188],[-2.4717266720048188,1.1312043837568135]]]-hessianF :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f (f a))-hessianF f = decomposeFunctor . jacobian (ComposeFunctor . jacobian (fmap decomposeMode . f . fmap composeMode))+hessianF :: (Traversable f, Functor g, Num a) => (forall s s'. (Reifies s Tape, Reifies s' Tape) => f (On (Reverse (Reverse a s') s)) -> g (On (Reverse (Reverse a s') s))) -> f a -> g (f (f a))+hessianF f = getCompose . jacobian (Compose . jacobian (fmap off . f . fmap On))+{-# INLINE hessianF #-}
src/Numeric/AD/Mode/Sparse.hs view
@@ -1,8 +1,7 @@ {-# LANGUAGE Rank2Types #-} ----------------------------------------------------------------------------- -- |--- Module      : Numeric.AD.Mode.Sparse--- Copyright   : (c) Edward Kmett 2010+-- Copyright   : (c) Edward Kmett 2010-2014 -- License     : BSD3 -- Maintainer  : ekmett@gmail.com -- Stability   : experimental@@ -13,41 +12,42 @@ -----------------------------------------------------------------------------  module Numeric.AD.Mode.Sparse-    (-    -- * Sparse Gradients-      grad-    , grad'-    , gradWith-    , gradWith'-    , grads--    -- * Sparse Jacobians (synonyms)-    , jacobian-    , jacobian'-    , jacobianWith-    , jacobianWith'-    , jacobians+  ( Sparse+  -- * Sparse Gradients+  , grad+  , grad'+  , gradWith+  , gradWith'+  -- * Variadic Gradients+  -- $vgrad+  , Grad+  , vgrad+  -- * Higher-Order Gradients+  , grads+  -- * Variadic Higher-Order Gradients+  , Grads+  , vgrads -    -- * Sparse Hessians-    , hessian-    , hessian'+  -- * Sparse Jacobians (synonyms)+  , jacobian+  , jacobian'+  , jacobianWith+  , jacobianWith'+  , jacobians -    , hessianF-    , hessianF'+  -- * Sparse Hessians+  , hessian+  , hessian' -    -- * Unsafe gradients-    , vgrad-    , vgrads+  , hessianF+  , hessianF' -    -- * Exposed Types-    , Grad-    , Grads-    ) where+  ) where  import Control.Comonad import Data.Traversable import Control.Comonad.Cofree-import Numeric.AD.Types+import Numeric.AD.Jet import Numeric.AD.Internal.Sparse import Numeric.AD.Internal.Combinators @@ -55,43 +55,43 @@ 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 :: (Traversable f, Num a) => (forall s. f (Sparse a s) -> Sparse a s) -> 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' :: (Traversable f, Num a) => (forall s. f (Sparse a s) -> Sparse a s) -> 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 :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. f (Sparse a s) -> Sparse a s) -> 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' :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. f (Sparse a s) -> Sparse a s) -> 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 :: (Traversable f, Functor g, Num a) => (forall s. f (Sparse a s) -> g (Sparse a s)) -> 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' :: (Traversable f, Functor g, Num a) => (forall s. f (Sparse a s) -> g (Sparse a s)) -> 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 :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. f (Sparse a s) -> g (Sparse a s)) -> 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' :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. f (Sparse a s) -> g (Sparse a s)) -> 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 :: (Traversable f, Num a) => (forall s. f (Sparse a s) -> Sparse a s) -> 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 :: (Traversable f, Functor g, Num a) => (forall s. f (Sparse a s) -> g (Sparse a s)) -> f a -> g (Cofree f a) jacobians f as = ds as <$> apply f as {-# INLINE jacobians #-} @@ -103,18 +103,27 @@ 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 :: (Traversable f, Num a) => (forall s. f (Sparse a s) -> Sparse a s) -> 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' :: (Traversable f, Num a) => (forall s. f (Sparse a s) -> Sparse a s) -> 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 :: (Traversable f, Functor g, Num a) => (forall s. f (Sparse a s) -> g (Sparse a s)) -> 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' :: (Traversable f, Functor g, Num a) => (forall s. f (Sparse a s) -> g (Sparse a s)) -> f a -> g (a, f (a, f a)) hessianF' f as = d2' <$> jacobians f as {-# INLINE hessianF' #-}++-- $vgrad+--+-- 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 'auto' calls you use when taking the gradient+-- of a function that takes gradients!
src/Numeric/AD/Mode/Tower.hs view
@@ -1,8 +1,8 @@-{-# LANGUAGE Rank2Types, BangPatterns #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE BangPatterns #-} ----------------------------------------------------------------------------- -- |--- Module      : Numeric.AD.Mode.Tower--- Copyright   : (c) Edward Kmett 2010+-- Copyright   : (c) Edward Kmett 2010-2014 -- License     : BSD3 -- Maintainer  : ekmett@gmail.com -- Stability   : experimental@@ -13,111 +13,103 @@ -----------------------------------------------------------------------------  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+  ( 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 :: Num a => (forall s. Tower a s -> Tower a s) -> 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 :: Num a => (forall s. Tower a s -> Tower a s) -> 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 :: (Functor f, Num a) => (forall s. Tower a s -> f (Tower a s)) -> 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 :: (Functor f, Num a) => (forall s. Tower a s -> f (Tower a s)) -> 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 _ _ [] = []+taylor :: Fractional a => (forall s. Tower a s -> Tower a s) -> 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 :: Fractional a => (forall s. Tower a s -> Tower a s) -> 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 :: Fractional a => (forall s. Tower a s -> Tower a s) -> 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 :: Fractional a => (forall s. Tower a s -> Tower a s) -> 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 :: Num a => (forall s. Tower a s -> Tower a s) -> 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' :: Num a => (forall s. Tower a s -> Tower a s) -> 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 :: (Functor f, Num a) => (forall s. f (Tower a s) -> Tower a s) -> 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' :: (Functor f, Num a) => (forall s. f (Tower a s) -> Tower a s) -> 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 :: (Functor f, Functor g, Num a) => (forall s. f (Tower a s) -> g (Tower a s)) -> 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' :: (Functor f, Functor g, Num a) => (forall s. f (Tower a s) -> g (Tower a s)) -> 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 :: (Functor f, Num a) => (forall s. f (Tower a s) -> Tower a s) -> 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 :: (Functor f, Num a) => (forall s. f (Tower a s) -> Tower a s) -> 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 :: (Functor f, Functor g, Num a) => (forall s. f (Tower a s) -> g (Tower a s)) -> 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 :: (Functor f, Functor g, Num a) => (forall s. f (Tower a s) -> g (Tower a s)) -> 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)
src/Numeric/AD/Newton.hs view
@@ -1,7 +1,10 @@-{-# LANGUAGE Rank2Types, BangPatterns, ScopedTypeVariables #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeFamilies #-} ----------------------------------------------------------------------------- -- |--- Module      :  Numeric.AD.Newton -- Copyright   :  (c) Edward Kmett 2010 -- License     :  BSD3 -- Maintainer  :  ekmett@gmail.com@@ -11,29 +14,34 @@ -----------------------------------------------------------------------------  module Numeric.AD.Newton-    (-    -- * Newton's Method (Forward AD)-      findZero-    , inverse-    , fixedPoint-    , extremum-    -- * Gradient Ascent/Descent (Reverse AD)-    , gradientDescent-    , gradientAscent-    , conjugateGradientDescent-    , conjugateGradientAscent-    ) where+  (+  -- * Newton's Method (Forward AD)+    findZero+  , inverse+  , fixedPoint+  , extremum+  -- * Gradient Ascent/Descent (Reverse AD)+  , gradientDescent+  , gradientAscent+  , conjugateGradientDescent+  , conjugateGradientAscent+  ) where  import Prelude hiding (all, mapM, sum)-import Data.Functor import Data.Foldable (all, sum)+import Data.Reflection (Reifies) import Data.Traversable-import Numeric.AD.Types+import Numeric.AD.Mode import Numeric.AD.Mode.Forward (diff, diff') import Numeric.AD.Mode.Reverse (grad, gradWith') import Numeric.AD.Internal.Combinators-import Numeric.AD.Internal.Composition+import Numeric.AD.Internal.Forward (Forward)+import Numeric.AD.Internal.On+import Numeric.AD.Internal.Reverse (Reverse, Tape) +-- $setup+-- >>> import Data.Complex+ -- | 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.) If the stream becomes constant@@ -44,10 +52,9 @@ -- >>> take 10 $ findZero (\x->x^2-4) 1 -- [1.0,2.5,2.05,2.000609756097561,2.0000000929222947,2.000000000000002,2.0] ----- >>> import Data.Complex -- >>> last $ take 10 $ findZero ((+1).(^2)) (1 :+ 1) -- 0.0 :+ 1.0-findZero :: (Fractional a, Eq a) => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]+findZero :: (Fractional a, Eq a) => (forall s. Forward a s -> Forward a s) -> a -> [a] findZero f = go where   go x = x : if x == xn then [] else go xn where     (y,y') = diff' f x@@ -63,7 +70,7 @@ -- -- >>> last $ take 10 $ inverse sqrt 1 (sqrt 10) -- 10.0-inverse :: (Fractional a, Eq a) => (forall s. Mode s => AD s a -> AD s a) -> a -> a -> [a]+inverse :: (Fractional a, Eq a) => (forall s. Forward a s -> Forward a s) -> a -> a -> [a] inverse f x0 y = findZero (\x -> f x - auto y) x0 {-# INLINE inverse  #-} @@ -76,7 +83,7 @@ -- -- >>> last $ take 10 $ fixedPoint cos 1 -- 0.7390851332151607-fixedPoint :: (Fractional a, Eq a) => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]+fixedPoint :: (Fractional a, Eq a) => (forall s. Forward a s -> Forward a s) -> a -> [a] fixedPoint f = findZero (\x -> f x - x) {-# INLINE fixedPoint #-} @@ -87,8 +94,8 @@ -- -- >>> last $ take 10 $ extremum cos 1 -- 0.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))+extremum :: (Fractional a, Eq a) => (forall s s'. On (Forward (Forward a s') s) -> On (Forward (Forward a s') s)) -> a -> [a]+extremum f = findZero (diff (off . f . On)) {-# INLINE extremum #-}  -- | The 'gradientDescent' function performs a multivariate@@ -98,44 +105,53 @@ -- 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 :: (Traversable f, Fractional a, Ord a) => (forall s. Reifies s Tape => f (Reverse a s) -> Reverse a s) -> 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+  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 #-}  -- | Perform a gradient descent using reverse mode automatic differentiation to compute the gradient.-gradientAscent :: (Traversable f, Fractional a, Ord a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> [f a]+gradientAscent :: (Traversable f, Fractional a, Ord a) => (forall s. Reifies s Tape => f (Reverse a s) -> Reverse a s) -> f a -> [f a] gradientAscent f = gradientDescent (negate . f) {-# INLINE gradientAscent #-} --- | Perform a conjugate gradient descent using reverse mode automatic differentiation to compute the gradient.-conjugateGradientDescent :: (Traversable f, Fractional a, Ord a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> [f a]-conjugateGradientDescent f x0 = takeWhile (all (\a -> a == a)) (go x0 d0 d0)+-- | Perform a conjugate gradient descent using reverse mode automatic differentiation to compute the gradient, and using forward-on-forward mode for computing extrema.+--+-- >>> let sq x = x * x+-- >>> let rosenbrock [x,y] = sq (1 - x) + 100 * sq (y - sq x)+-- >>> rosenbrock [0,0]+-- 1+-- >>> rosenbrock (conjugateGradientDescent rosenbrock [0, 0] !! 5) < 0.1+-- True+conjugateGradientDescent :: (Traversable f, Ord a, Fractional a) => (forall t. (Mode t, a ~ Scalar t, Num t) => f t -> t) -> f a -> [f a]+conjugateGradientDescent f = conjugateGradientAscent (negate . f)+{-# INLINE conjugateGradientDescent #-}++-- | Perform a conjugate gradient ascent using reverse mode automatic differentiation to compute the gradient.+conjugateGradientAscent :: (Traversable f, Ord a, Fractional a) => (forall t. (Mode t, a ~ Scalar t, Num t) => f t -> t) -> f a -> [f a]+conjugateGradientAscent f x0 = takeWhile (all (\a -> a == a)) (go x0 d0 d0 delta0)   where     dot x y = sum $ zipWithT (*) x y-    d0 = negate <$> grad f x0-    go xi ri di = xi : go xi1 ri1 di1+    d0 = grad f x0+    delta0 = dot d0 d0+    go xi _ri di deltai = xi : go xi1 ri1 di1 deltai1       where-        ai  = last $ take 20 $ extremum (\a -> f $ zipWithT (\x d -> auto x + a * auto d) xi di) 0+        ai = last $ take 20 $ extremum (\a -> f $ zipWithT (\x d -> auto x + a * auto d) xi di) 0         xi1 = zipWithT (\x d -> x + ai*d) xi di-        ri1 = negate <$> grad f xi1-        bi1 = max 0 $ dot ri1 (zipWithT (-) ri1 ri) / dot ri1 ri1-        di1 = zipWithT (\r d -> r * bi1*d) ri1 di-{-# INLINE conjugateGradientDescent #-}---- | Perform a conjugate gradient ascent using reverse mode automatic differentiation to compute the gradient.-conjugateGradientAscent :: (Traversable f, Fractional a, Ord a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> [f a]-conjugateGradientAscent f = conjugateGradientDescent (negate . f)+        ri1 = grad f xi1+        deltai1 = dot ri1 ri1+        bi1 = deltai1 / deltai+        di1 = zipWithT (\r d -> r + bi1 * d) ri1 di {-# INLINE conjugateGradientAscent #-}
− src/Numeric/AD/Types.hs 
@@ -1,50 +0,0 @@-{-# LANGUAGE Rank2Types #-}--------------------------------------------------------------------------------- |--- Module      :  Numeric.AD.Types--- Copyright   :  (c) Edward Kmett 2010-12--- 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---- | Evaluate a scalar-to-scalar function in the trivial identity AD mode.-lowerUU :: (forall s. Mode s => AD s a -> AD s a) -> a -> a-lowerUU f = unprobe . f . probe-{-# INLINE lowerUU #-}---- | Evaluate a scalar-to-nonscalar function in the trivial identity AD mode.-lowerUF :: (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f a-lowerUF f = unprobed . f . probe-{-# INLINE lowerUF #-}---- | Evaluate a nonscalar-to-scalar function in the trivial identity AD mode.-lowerFU :: (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> a-lowerFU f = unprobe . f . probed-{-# INLINE lowerFU #-}---- | Evaluate a nonscalar-to-nonscalar function in the trivial identity AD mode.-lowerFF :: (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g a-lowerFF f = unprobed . f . probed-{-# INLINE lowerFF #-}
− src/Numeric/AD/Variadic.hs 
@@ -1,28 +0,0 @@--------------------------------------------------------------------------------- |--- 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.Kahn-import Numeric.AD.Variadic.Sparse (Grads, vgrads)
− src/Numeric/AD/Variadic/Kahn.hs 
@@ -1,26 +0,0 @@--------------------------------------------------------------------------------- |--- Module      :  Numeric.AD.Variadic.Kahn--- Copyright   :  (c) Edward Kmett 2010-2012--- License     :  BSD3--- Maintainer  :  ekmett@gmail.com--- Stability   :  experimental--- Portability :  non-portable------ Variadic combinators for reverse-mode automatic differentiation.------ Unfortunately, variadicity comes at the expense of being able to use--- quantification to avoid sensitivity confusion, so be careful when--- counting the number of @lift@ you use when taking the gradient of a--- function that takes gradients!-----------------------------------------------------------------------------------module Numeric.AD.Variadic.Kahn-    (-    -- * Unsafe Variadic Gradient-      vgrad, vgrad'-    , Grad-    ) where--import Numeric.AD.Internal.Kahn
− src/Numeric/AD/Variadic/Sparse.hs 
@@ -1,26 +0,0 @@--------------------------------------------------------------------------------- |--- 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
tests/doctests.hs view
@@ -14,6 +14,8 @@   : "-idist/build/autogen"   : "-optP-include"   : "-optPdist/build/autogen/cabal_macros.h"+  : "-optP-I"+  : "-optPinclude"   : "-hide-all-packages"   : map ("-package="++) deps ++ sources