diff --git a/.ghci b/.ghci --- a/.ghci +++ b/.ghci @@ -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 diff --git a/.travis.yml b/.travis.yml --- a/.travis.yml +++ b/.travis.yml @@ -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" diff --git a/CHANGELOG.markdown b/CHANGELOG.markdown --- a/CHANGELOG.markdown +++ b/CHANGELOG.markdown @@ -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 ----- diff --git a/README.markdown b/README.markdown --- a/README.markdown +++ b/README.markdown @@ -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 -------- diff --git a/ad.cabal b/ad.cabal --- a/ad.cabal +++ b/ad.cabal @@ -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 diff --git a/bench/BlackScholes.hs b/bench/BlackScholes.hs new file mode 100644 --- /dev/null +++ b/bench/BlackScholes.hs @@ -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 +-- ] + ] diff --git a/src/Numeric/AD.hs b/src/Numeric/AD.hs --- a/src/Numeric/AD.hs +++ b/src/Numeric/AD.hs @@ -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 , 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! diff --git a/src/Numeric/AD/Halley.hs b/src/Numeric/AD/Halley.hs --- a/src/Numeric/AD/Halley.hs +++ b/src/Numeric/AD/Halley.hs @@ -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 #-} diff --git a/src/Numeric/AD/Internal/Classes.hs b/src/Numeric/AD/Internal/Classes.hs deleted file mode 100644 --- a/src/Numeric/AD/Internal/Classes.hs +++ /dev/null @@ -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' diff --git a/src/Numeric/AD/Internal/Combinators.hs b/src/Numeric/AD/Internal/Combinators.hs --- a/src/Numeric/AD/Internal/Combinators.hs +++ b/src/Numeric/AD/Internal/Combinators.hs @@ -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 diff --git a/src/Numeric/AD/Internal/Composition.hs b/src/Numeric/AD/Internal/Composition.hs deleted file mode 100644 --- a/src/Numeric/AD/Internal/Composition.hs +++ /dev/null @@ -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 - diff --git a/src/Numeric/AD/Internal/Dense.hs b/src/Numeric/AD/Internal/Dense.hs --- a/src/Numeric/AD/Internal/Dense.hs +++ b/src/Numeric/AD/Internal/Dense.hs @@ -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" diff --git a/src/Numeric/AD/Internal/Forward.hs b/src/Numeric/AD/Internal/Forward.hs --- a/src/Numeric/AD/Internal/Forward.hs +++ b/src/Numeric/AD/Internal/Forward.hs @@ -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 diff --git a/src/Numeric/AD/Internal/Forward/Double.hs b/src/Numeric/AD/Internal/Forward/Double.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AD/Internal/Forward/Double.hs @@ -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 diff --git a/src/Numeric/AD/Internal/Identity.hs b/src/Numeric/AD/Internal/Identity.hs --- a/src/Numeric/AD/Internal/Identity.hs +++ b/src/Numeric/AD/Internal/Identity.hs @@ -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) diff --git a/src/Numeric/AD/Internal/Jet.hs b/src/Numeric/AD/Internal/Jet.hs deleted file mode 100644 --- a/src/Numeric/AD/Internal/Jet.hs +++ /dev/null @@ -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 #-} diff --git a/src/Numeric/AD/Internal/Kahn.hs b/src/Numeric/AD/Internal/Kahn.hs --- a/src/Numeric/AD/Internal/Kahn.hs +++ b/src/Numeric/AD/Internal/Kahn.hs @@ -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 + +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 diff --git a/src/Numeric/AD/Internal/On.hs b/src/Numeric/AD/Internal/On.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AD/Internal/On.hs @@ -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) diff --git a/src/Numeric/AD/Internal/Reverse.hs b/src/Numeric/AD/Internal/Reverse.hs --- a/src/Numeric/AD/Internal/Reverse.hs +++ b/src/Numeric/AD/Internal/Reverse.hs @@ -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 diff --git a/src/Numeric/AD/Internal/Sparse.hs b/src/Numeric/AD/Internal/Sparse.hs --- a/src/Numeric/AD/Internal/Sparse.hs +++ b/src/Numeric/AD/Internal/Sparse.hs @@ -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 #-} - diff --git a/src/Numeric/AD/Internal/Tower.hs b/src/Numeric/AD/Internal/Tower.hs --- a/src/Numeric/AD/Internal/Tower.hs +++ b/src/Numeric/AD/Internal/Tower.hs @@ -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 diff --git a/src/Numeric/AD/Internal/Types.hs b/src/Numeric/AD/Internal/Types.hs deleted file mode 100644 --- a/src/Numeric/AD/Internal/Types.hs +++ /dev/null @@ -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 diff --git a/src/Numeric/AD/Internal/Var.hs b/src/Numeric/AD/Internal/Var.hs deleted file mode 100644 --- a/src/Numeric/AD/Internal/Var.hs +++ /dev/null @@ -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 diff --git a/src/Numeric/AD/Jacobian.hs b/src/Numeric/AD/Jacobian.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AD/Jacobian.hs @@ -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 diff --git a/src/Numeric/AD/Jet.hs b/src/Numeric/AD/Jet.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AD/Jet.hs @@ -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 diff --git a/src/Numeric/AD/Mode.hs b/src/Numeric/AD/Mode.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AD/Mode.hs @@ -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 diff --git a/src/Numeric/AD/Mode/Directed.hs b/src/Numeric/AD/Mode/Directed.hs --- a/src/Numeric/AD/Mode/Directed.hs +++ b/src/Numeric/AD/Mode/Directed.hs @@ -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' #-} diff --git a/src/Numeric/AD/Mode/Forward.hs b/src/Numeric/AD/Mode/Forward.hs --- a/src/Numeric/AD/Mode/Forward.hs +++ b/src/Numeric/AD/Mode/Forward.hs @@ -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' #-} diff --git a/src/Numeric/AD/Mode/Forward/Double.hs b/src/Numeric/AD/Mode/Forward/Double.hs new file mode 100644 --- /dev/null +++ b/src/Numeric/AD/Mode/Forward/Double.hs @@ -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' #-} diff --git a/src/Numeric/AD/Mode/Kahn.hs b/src/Numeric/AD/Mode/Kahn.hs --- a/src/Numeric/AD/Mode/Kahn.hs +++ b/src/Numeric/AD/Mode/Kahn.hs @@ -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! diff --git a/src/Numeric/AD/Mode/Reverse.hs b/src/Numeric/AD/Mode/Reverse.hs --- a/src/Numeric/AD/Mode/Reverse.hs +++ b/src/Numeric/AD/Mode/Reverse.hs @@ -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 #-} diff --git a/src/Numeric/AD/Mode/Sparse.hs b/src/Numeric/AD/Mode/Sparse.hs --- a/src/Numeric/AD/Mode/Sparse.hs +++ b/src/Numeric/AD/Mode/Sparse.hs @@ -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! diff --git a/src/Numeric/AD/Mode/Tower.hs b/src/Numeric/AD/Mode/Tower.hs --- a/src/Numeric/AD/Mode/Tower.hs +++ b/src/Numeric/AD/Mode/Tower.hs @@ -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) diff --git a/src/Numeric/AD/Newton.hs b/src/Numeric/AD/Newton.hs --- a/src/Numeric/AD/Newton.hs +++ b/src/Numeric/AD/Newton.hs @@ -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 #-} diff --git a/src/Numeric/AD/Types.hs b/src/Numeric/AD/Types.hs deleted file mode 100644 --- a/src/Numeric/AD/Types.hs +++ /dev/null @@ -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 #-} diff --git a/src/Numeric/AD/Variadic.hs b/src/Numeric/AD/Variadic.hs deleted file mode 100644 --- a/src/Numeric/AD/Variadic.hs +++ /dev/null @@ -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) diff --git a/src/Numeric/AD/Variadic/Kahn.hs b/src/Numeric/AD/Variadic/Kahn.hs deleted file mode 100644 --- a/src/Numeric/AD/Variadic/Kahn.hs +++ /dev/null @@ -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 diff --git a/src/Numeric/AD/Variadic/Sparse.hs b/src/Numeric/AD/Variadic/Sparse.hs deleted file mode 100644 --- a/src/Numeric/AD/Variadic/Sparse.hs +++ /dev/null @@ -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 diff --git a/tests/doctests.hs b/tests/doctests.hs --- a/tests/doctests.hs +++ b/tests/doctests.hs @@ -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