ad 4.4 → 4.4.1

raw patch · 13 files changed

+672/−85 lines, 13 filesdep ~basedep ~doctestPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependency ranges changed: base, doctest

API changes (from Hackage documentation)

Files

• .travis.yml +49/−57
• CHANGELOG.markdown +11/−0
• ad.cabal +8/−3
• src/Numeric/AD.hs +9/−2
• src/Numeric/AD/Internal/Forward/Double.hs +4/−4
• src/Numeric/AD/Internal/Kahn.hs +2/−2
• src/Numeric/AD/Internal/Reverse/Double.hs +356/−0
• src/Numeric/AD/Jet.hs +1/−1
• src/Numeric/AD/Mode/Forward.hs +4/−4
• src/Numeric/AD/Mode/Forward/Double.hs +4/−4
• src/Numeric/AD/Mode/Reverse/Double.hs +216/−0
• src/Numeric/AD/Rank1/Forward.hs +4/−4
• src/Numeric/AD/Rank1/Forward/Double.hs +4/−4

@@ -2,11 +2,17 @@ # #   haskell-ci '--output=.travis.yml' '--config=cabal.haskell-ci' 'cabal.project' #+# To regenerate the script (for example after adjusting tested-with) run+#+#   haskell-ci regenerate+# # For more information, see https://github.com/haskell-CI/haskell-ci #-# version: 0.5.20180830+# version: 0.10 #+version: ~> 1.0 language: c+os: linux dist: xenial git:   # whether to recursively clone submodules@@ -17,11 +23,12 @@       - irc.freenode.org#haskell-lens     skip_join: true     template:-      - "\"\\x0313ad\\x03/\\x0306%{branch}\\x03 \\x0314%{commit}\\x03 %{build_url} %{message}\""+      - "\x0313ad\x03/\x0306%{branch}\x03 \x0314%{commit}\x03 %{build_url} %{message}" cache:   directories:     - $HOME/.cabal/packages     - $HOME/.cabal/store+    - $HOME/.hlint before_cache:   - rm -fv $CABALHOME/packages/hackage.haskell.org/build-reports.log   # remove files that are regenerated by 'cabal update'@@ -31,30 +38,38 @@   - rm -fv $CABALHOME/packages/hackage.haskell.org/01-index.tar   - rm -fv $CABALHOME/packages/hackage.haskell.org/01-index.tar.idx   - rm -rfv $CABALHOME/packages/head.hackage-matrix:+jobs:   include:-    - compiler: ghc-8.8.1-      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-8.8.1","cabal-install-3.0"]}}+    - compiler: ghc-8.10.1+      addons: {"apt":{"sources":[{"sourceline":"deb http://ppa.launchpad.net/hvr/ghc/ubuntu xenial main","key_url":"https://keyserver.ubuntu.com/pks/lookup?op=get&search=0x063dab2bdc0b3f9fcebc378bff3aeacef6f88286"}],"packages":["ghc-8.10.1","cabal-install-3.2"]}}+      os: linux+    - compiler: ghc-8.8.3+      addons: {"apt":{"sources":[{"sourceline":"deb http://ppa.launchpad.net/hvr/ghc/ubuntu xenial main","key_url":"https://keyserver.ubuntu.com/pks/lookup?op=get&search=0x063dab2bdc0b3f9fcebc378bff3aeacef6f88286"}],"packages":["ghc-8.8.3","cabal-install-3.2"]}}+      os: linux     - compiler: ghc-8.6.5-      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-8.6.5","cabal-install-3.0"]}}+      addons: {"apt":{"sources":[{"sourceline":"deb http://ppa.launchpad.net/hvr/ghc/ubuntu xenial main","key_url":"https://keyserver.ubuntu.com/pks/lookup?op=get&search=0x063dab2bdc0b3f9fcebc378bff3aeacef6f88286"}],"packages":["ghc-8.6.5","cabal-install-3.2"]}}+      os: linux     - compiler: ghc-8.4.4-      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-8.4.4","cabal-install-3.0"]}}+      addons: {"apt":{"sources":[{"sourceline":"deb http://ppa.launchpad.net/hvr/ghc/ubuntu xenial main","key_url":"https://keyserver.ubuntu.com/pks/lookup?op=get&search=0x063dab2bdc0b3f9fcebc378bff3aeacef6f88286"}],"packages":["ghc-8.4.4","cabal-install-3.2"]}}+      os: linux     - compiler: ghc-8.2.2-      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-8.2.2","cabal-install-3.0"]}}+      addons: {"apt":{"sources":[{"sourceline":"deb http://ppa.launchpad.net/hvr/ghc/ubuntu xenial main","key_url":"https://keyserver.ubuntu.com/pks/lookup?op=get&search=0x063dab2bdc0b3f9fcebc378bff3aeacef6f88286"}],"packages":["ghc-8.2.2","cabal-install-3.2"]}}+      os: linux     - compiler: ghc-8.0.2-      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-8.0.2","cabal-install-3.0"]}}+      addons: {"apt":{"sources":[{"sourceline":"deb http://ppa.launchpad.net/hvr/ghc/ubuntu xenial main","key_url":"https://keyserver.ubuntu.com/pks/lookup?op=get&search=0x063dab2bdc0b3f9fcebc378bff3aeacef6f88286"}],"packages":["ghc-8.0.2","cabal-install-3.2"]}}+      os: linux     - compiler: ghc-7.10.3-      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-7.10.3","cabal-install-3.0"]}}+      addons: {"apt":{"sources":[{"sourceline":"deb http://ppa.launchpad.net/hvr/ghc/ubuntu xenial main","key_url":"https://keyserver.ubuntu.com/pks/lookup?op=get&search=0x063dab2bdc0b3f9fcebc378bff3aeacef6f88286"}],"packages":["ghc-7.10.3","cabal-install-3.2"]}}+      os: linux     - compiler: ghc-7.8.4-      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-7.8.4","cabal-install-3.0"]}}+      addons: {"apt":{"sources":[{"sourceline":"deb http://ppa.launchpad.net/hvr/ghc/ubuntu xenial main","key_url":"https://keyserver.ubuntu.com/pks/lookup?op=get&search=0x063dab2bdc0b3f9fcebc378bff3aeacef6f88286"}],"packages":["ghc-7.8.4","cabal-install-3.2"]}}+      os: linux     - compiler: ghc-7.6.3-      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-7.6.3","cabal-install-3.0"]}}+      addons: {"apt":{"sources":[{"sourceline":"deb http://ppa.launchpad.net/hvr/ghc/ubuntu xenial main","key_url":"https://keyserver.ubuntu.com/pks/lookup?op=get&search=0x063dab2bdc0b3f9fcebc378bff3aeacef6f88286"}],"packages":["ghc-7.6.3","cabal-install-3.2"]}}+      os: linux     - compiler: ghc-7.4.2-      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-7.4.2","cabal-install-3.0"]}}-    - compiler: ghc-head-      addons: {"apt":{"sources":["hvr-ghc"],"packages":["ghc-head","cabal-install-head"]}}-  allow_failures:-    - compiler: ghc-head+      addons: {"apt":{"sources":[{"sourceline":"deb http://ppa.launchpad.net/hvr/ghc/ubuntu xenial main","key_url":"https://keyserver.ubuntu.com/pks/lookup?op=get&search=0x063dab2bdc0b3f9fcebc378bff3aeacef6f88286"}],"packages":["ghc-7.4.2","cabal-install-3.2"]}}+      os: linux before_install:   - HC=$(echo "/opt/$CC/bin/ghc" | sed 's/-/\//')   - WITHCOMPILER="-w $HC"@@ -67,33 +82,11 @@   - TOP=$(pwd)   - "HCNUMVER=$(${HC} --numeric-version|perl -ne '/^(\\d+)\\.(\\d+)\\.(\\d+)(\\.(\\d+))?$/; print(10000 * $1 + 100 * $2 + ($3 == 0 ? $5 != 1 : $3))')"   - echo $HCNUMVER-  - CABAL="$CABAL -vnormal+nowrap+markoutput"+  - CABAL="$CABAL -vnormal+nowrap"   - set -o pipefail-  - |-    echo 'function blue(s) { printf "\033[0;34m" s "\033[0m " }'           >> .colorful.awk-    echo 'BEGIN { state = "output"; }'                                     >> .colorful.awk-    echo '/^-----BEGIN CABAL OUTPUT-----$/ { state = "cabal" }'            >> .colorful.awk-    echo '/^-----END CABAL OUTPUT-----$/ { state = "output" }'             >> .colorful.awk-    echo '!/^(-----BEGIN CABAL OUTPUT-----|-----END CABAL OUTPUT-----)/ {' >> .colorful.awk-    echo '  if (state == "cabal") {'                                       >> .colorful.awk-    echo '    print blue($0)'                                              >> .colorful.awk-    echo '  } else {'                                                      >> .colorful.awk-    echo '    print $0'                                                    >> .colorful.awk-    echo '  }'                                                             >> .colorful.awk-    echo '}'                                                               >> .colorful.awk-  - cat .colorful.awk-  - |-    color_cabal_output () {-      awk -f $TOP/.colorful.awk-    }-  - echo text | color_cabal_output-install:-  - ${CABAL} --version-  - echo "$(${HC} --version) [$(${HC} --print-project-git-commit-id 2> /dev/null || echo '?')]"   - TEST=--enable-tests   - BENCH=--enable-benchmarks   - HEADHACKAGE=false-  - if [ $HCNUMVER -gt 80801 ] ; then HEADHACKAGE=true ; fi   - rm -f $CABALHOME/config   - |     echo "verbose: normal +nowrap +markoutput"          >> $CABALHOME/config@@ -111,17 +104,12 @@     echo "  prefix: $CABALHOME"                         >> $CABALHOME/config     echo "repository hackage.haskell.org"               >> $CABALHOME/config     echo "  url: http://hackage.haskell.org/"           >> $CABALHOME/config+install:+  - ${CABAL} --version+  - echo "$(${HC} --version) [$(${HC} --print-project-git-commit-id 2> /dev/null || echo '?')]"   - |-    if $HEADHACKAGE; then-    echo "allow-newer: $($HCPKG list --simple-output | sed -E 's/([a-zA-Z-]+)-[0-9.]+/*:\1/g')" >> $CABALHOME/config-    echo "repository head.hackage.ghc.haskell.org"                                        >> $CABALHOME/config-    echo "   url: https://ghc.gitlab.haskell.org/head.hackage/"                           >> $CABALHOME/config-    echo "   secure: True"                                                                >> $CABALHOME/config-    echo "   root-keys: 7541f32a4ccca4f97aea3b22f5e593ba2c0267546016b992dfadcd2fe944e55d" >> $CABALHOME/config-    echo "              26021a13b401500c8eb2761ca95c61f2d625bfef951b939a8124ed12ecf07329" >> $CABALHOME/config-    echo "              f76d08be13e9a61a377a85e2fb63f4c5435d40f8feb3e12eb05905edb8cdea89" >> $CABALHOME/config-    echo "   key-threshold: 3"                                                            >> $CABALHOME/config-    fi+    echo "program-default-options"                >> $CABALHOME/config+    echo "  ghc-options: $GHCJOBS +RTS -M6G -RTS" >> $CABALHOME/config   - cat $CABALHOME/config   - rm -fv cabal.project cabal.project.local cabal.project.freeze   - travis_retry ${CABAL} v2-update -v@@ -130,19 +118,21 @@   - touch cabal.project   - |     echo "packages: ." >> cabal.project+  - if [ $HCNUMVER -ge 80200 ] ; then echo 'package ad' >> cabal.project ; fi+  - "if [ $HCNUMVER -ge 80200 ] ; then echo '  ghc-options: -Werror=missing-methods' >> cabal.project ; fi"   - |   - "for pkg in $($HCPKG list --simple-output); do echo $pkg | sed 's/-[^-]*$//' | (grep -vE -- '^(ad)$' || true) | sed 's/^/constraints: /' | sed 's/$/ installed/' >> cabal.project.local; done"   - cat cabal.project || true   - cat cabal.project.local || true   - if [ -f "./configure.ac" ]; then (cd "." && autoreconf -i); fi-  - ${CABAL} v2-freeze $WITHCOMPILER ${TEST} ${BENCH} | color_cabal_output+  - ${CABAL} v2-freeze $WITHCOMPILER ${TEST} ${BENCH}   - "cat cabal.project.freeze | sed -E 's/^(constraints: *| *)//' | sed 's/any.//'"   - rm  cabal.project.freeze-  - ${CABAL} v2-build $WITHCOMPILER ${TEST} ${BENCH} --dep -j2 all | color_cabal_output+  - travis_wait 40 ${CABAL} v2-build $WITHCOMPILER ${TEST} ${BENCH} --dep -j2 all script:   - DISTDIR=$(mktemp -d /tmp/dist-test.XXXX)   # Packaging...-  - ${CABAL} v2-sdist all | color_cabal_output+  - ${CABAL} v2-sdist all   # Unpacking...   - mv dist-newstyle/sdist/*.tar.gz ${DISTDIR}/   - cd ${DISTDIR} || false@@ -154,19 +144,21 @@   - touch cabal.project   - |     echo "packages: ${PKGDIR_ad}" >> cabal.project+  - if [ $HCNUMVER -ge 80200 ] ; then echo 'package ad' >> cabal.project ; fi+  - "if [ $HCNUMVER -ge 80200 ] ; then echo '  ghc-options: -Werror=missing-methods' >> cabal.project ; fi"   - |   - "for pkg in $($HCPKG list --simple-output); do echo $pkg | sed 's/-[^-]*$//' | (grep -vE -- '^(ad)$' || true) | sed 's/^/constraints: /' | sed 's/$/ installed/' >> cabal.project.local; done"   - cat cabal.project || true   - cat cabal.project.local || true   # Building with tests and benchmarks...   # build & run tests, build benchmarks-  - ${CABAL} v2-build $WITHCOMPILER ${TEST} ${BENCH} all | color_cabal_output+  - ${CABAL} v2-build $WITHCOMPILER ${TEST} ${BENCH} all   # Testing...-  - ${CABAL} v2-test $WITHCOMPILER ${TEST} ${BENCH} all | color_cabal_output+  - ${CABAL} v2-test $WITHCOMPILER ${TEST} ${BENCH} all   # cabal check...   - (cd ${PKGDIR_ad} && ${CABAL} -vnormal check)   # haddock...-  - ${CABAL} v2-haddock $WITHCOMPILER --with-haddock $HADDOCK ${TEST} ${BENCH} all | color_cabal_output+  - ${CABAL} v2-haddock $WITHCOMPILER --with-haddock $HADDOCK ${TEST} ${BENCH} all -# REGENDATA ["--output=.travis.yml","--config=cabal.haskell-ci","cabal.project"]+# REGENDATA ("0.10",["--output=.travis.yml","--config=cabal.haskell-ci","cabal.project"]) # EOF

@@ -1,3 +1,14 @@+4.4.1 [2020.10.13]+------------------+* Change the fixity of `:-` in `Numeric.AD.Jet` to be right-associative.+  Previously, it was `infixl`, which made things like `x :- y :- z` nearly+  unusable.+* Fix backpropagation error in Kahn mode.+* Fix bugs in the `Erf` instance for `ForwardDouble`.+* Add `Numeric.AD.Mode.Reverse.Double`, a variant of `Numeric.AD.Mode.Reverse`+  that is specialized to `Double`.+* Re-export `Jet(..)`, `headJet`, `tailJet` and `jet` from `Numeric.AD`.+ 4.4 [2020.02.03] ---------------- * Generalize the type of `stochasticGradientDescent`:

@@ -1,5 +1,5 @@ name:          ad-version:       4.4+version:       4.4.1 license:       BSD3 license-File:  LICENSE copyright:     (c) Edward Kmett 2010-2015,@@ -20,7 +20,8 @@              , GHC == 8.2.2              , GHC == 8.4.4              , GHC == 8.6.5-             , GHC == 8.8.1+             , GHC == 8.8.3+             , GHC == 8.10.1 synopsis:      Automatic Differentiation extra-source-files:   .gitignore@@ -148,6 +149,7 @@     Numeric.AD.Internal.On     Numeric.AD.Internal.Or     Numeric.AD.Internal.Reverse+    Numeric.AD.Internal.Reverse.Double     Numeric.AD.Internal.Sparse     Numeric.AD.Internal.Tower     Numeric.AD.Internal.Type@@ -158,6 +160,7 @@     Numeric.AD.Mode.Forward.Double     Numeric.AD.Mode.Kahn     Numeric.AD.Mode.Reverse+    Numeric.AD.Mode.Reverse.Double     Numeric.AD.Mode.Sparse     Numeric.AD.Mode.Tower     Numeric.AD.Newton@@ -176,6 +179,8 @@     Numeric.AD.Internal.Combinators    ghc-options: -Wall+  if impl(ghc >= 8.6)+    ghc-options: -Wno-star-is-type    ghc-options: -fspec-constr -fdicts-cheap -O2 @@ -188,7 +193,7 @@     ad,     base,     directory,-    doctest >= 0.9.0.1 && < 0.17,+    doctest >= 0.9.0.1 && < 0.18,     filepath   ghc-options: -Wall -threaded   hs-source-dirs: tests

@@ -129,6 +129,12 @@   , conjugateGradientDescent   , conjugateGradientAscent   , stochasticGradientDescent++  -- * Working with towers+  , Jet(..)+  , headJet+  , tailJet+  , jet   ) where  import Data.Functor.Compose@@ -143,6 +149,7 @@ import Numeric.AD.Internal.Sparse (Sparse, Grads, vgrads)  import Numeric.AD.Internal.Type+import Numeric.AD.Jet import Numeric.AD.Mode  import qualified Numeric.AD.Rank1.Forward as Forward1@@ -176,7 +183,7 @@ -- -- >>> import Numeric.AD.Internal.Doctest --- | @'hessianProduct' f wv@ computes the product of the hessian @H@ of a non-scalar-to-scalar function @f@ at @w = 'fst' <$> wv@ with a vector @v = snd <$> wv@ using \"Pearlmutter\'s method\" from <http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.29.6143>, which states:+-- | @'hessianProduct' f wv@ computes the product of the hessian @H@ of a non-scalar-to-scalar function @f@ at @w = 'fst' '<$>' wv@ with a vector @v = snd '<$>' wv@ using \"Pearlmutter\'s method\" from <http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.29.6143>, which states: -- -- > H v = (d/dr) grad_w (w + r v) | r = 0 --@@ -185,7 +192,7 @@ hessianProduct :: (Traversable f, Num a) => (forall s. Reifies s Tape => f (On (Reverse s (Forward a))) -> On (Reverse s (Forward a))) -> f (a, a) -> f a hessianProduct f = Forward1.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.+-- | @'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 --

@@ -194,12 +194,12 @@ instance Erf ForwardDouble 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)+  normcdf = lift1 normcdf $ \x -> (recip $ sqrt (2 * pi)) * exp (- x * x / 2)  instance InvErf ForwardDouble 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)+  inverf = lift1_ inverf $ \x _ -> sqrt pi / 2 * exp (x * x)+  inverfc = lift1_ inverfc $ \x _ -> negate (sqrt pi / 2) * exp (x * x)+  invnormcdf = lift1_ invnormcdf $ \x _ -> sqrt (2 * pi) * exp (x * x / 2)  bind :: Traversable f => (f ForwardDouble -> b) -> f Double -> f b bind f as = snd $ mapAccumL outer (0 :: Int) as where

@@ -231,9 +231,9 @@   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 :: Tape a Int -> [Int]   successors (Unary _ _ b) = [b]-  successors (Binary _ _ _ b c) = [b,c]+  successors (Binary _ _ _ b c) = if b == c then [b] else [b,c]   successors _ = []  -- | Return an 'Array' of 'partials' given bounds for the variable IDs.

@@ -0,0 +1,356 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE UndecidableInstances #-}+{-# OPTIONS_GHC -fno-full-laziness #-}+{-# OPTIONS_HADDOCK not-home #-}++module Numeric.AD.Internal.Reverse.Double+  ( ReverseDouble(..)+  , 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, findWithDefault)+import Data.Number.Erf+import Data.Proxy+import Data.Reflection+#if __GLASGOW_HASKELL__ < 710+import Data.Traversable (Traversable, mapM)+#else+import Data.Traversable (mapM)+#endif+import Data.Typeable+import Numeric.AD.Internal.Combinators+import Numeric.AD.Internal.Identity+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+data Cells where+  Nil    :: Cells+  Unary  :: {-# UNPACK #-} !Int -> Double -> Cells -> Cells+  Binary :: {-# UNPACK #-} !Int -> {-# UNPACK #-} !Int -> Double -> Double -> Cells -> Cells+#endif++dropCells :: Int -> Cells -> Cells+dropCells 0 xs = xs+dropCells _ Nil = Nil+dropCells n (Unary _ _ xs)      = (dropCells $! n - 1) xs+dropCells n (Binary _ _ _ _ xs) = (dropCells $! n - 1) xs++data Head = Head {-# UNPACK #-} !Int Cells++newtype Tape = Tape { getTape :: IORef Head }++un :: Int -> Double -> Head -> (Head, Int)+un i di (Head r t) = h `seq` r' `seq` (h, r') where+  r' = r + 1+  h = Head r' (Unary i di t)+{-# INLINE un #-}++bin :: Int -> Int -> Double -> Double -> Head -> (Head, Int)+bin i j di dj (Head r t) = h `seq` r' `seq` (h, r') where+  r' = r + 1+  h = Head r' (Binary i j di dj t)+{-# INLINE bin #-}++modifyTape :: Reifies s Tape => p s -> (Head -> (Head, r)) -> IO r+modifyTape p = atomicModifyIORef (getTape (reflect p))+{-# INLINE modifyTape #-}++-- | 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. Reifies s Tape => (Double -> Double) -> Double -> Int -> Double -> ReverseDouble s+unarily f di i b = ReverseDouble (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. Reifies s Tape => (Double -> Double -> Double) -> Double -> Double -> Int -> Double -> Int -> Double -> ReverseDouble s+binarily f di dj i b j c = ReverseDouble (unsafePerformIO (modifyTape (Proxy :: Proxy s) (bin i j di dj))) $! f b c+{-# INLINE binarily #-}++#ifndef HLINT+data ReverseDouble s where+  Zero :: ReverseDouble s+  Lift :: Double -> ReverseDouble s+  ReverseDouble :: {-# UNPACK #-} !Int -> Double -> ReverseDouble s+  deriving (Show, Typeable)+#endif++instance (Reifies s Tape) => Mode (ReverseDouble s) where+  type Scalar (ReverseDouble s) = Double++  isKnownZero Zero = True+  isKnownZero _    = False++  isKnownConstant ReverseDouble{} = False+  isKnownConstant _ = True++  auto = Lift+  zero = Zero+  a *^ b = lift1 (a *) (\_ -> auto a) b+  a ^* b = lift1 (* b) (\_ -> auto b) a+  a ^/ b = lift1 (/ b) (\_ -> auto (recip b)) a++(<+>) :: (Reifies s Tape) => ReverseDouble s -> ReverseDouble s -> ReverseDouble s+(<+>)  = binary (+) 1 1++(<**>) :: (Reifies s Tape) => ReverseDouble s -> ReverseDouble s -> ReverseDouble 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 * xi ** (yi - 1), z * log xi)) x y++primal :: ReverseDouble s -> Double+primal Zero = 0+primal (Lift a) = a+primal (ReverseDouble _ a) = a++instance (Reifies s Tape) => Jacobian (ReverseDouble s) where+  type D (ReverseDouble s) = Id Double++  unary f _         (Zero)   = Lift (f 0)+  unary f _         (Lift a) = Lift (f a)+  unary f (Id dadi) (ReverseDouble i b) = unarily f dadi i b++  lift1 f df b = unary f (df (Id pb)) b where+    pb = primal b++  lift1_ f df b = unary (const a) (df (Id a) (Id pb)) b where+    pb = primal b+    a = f pb++  binary f _         _         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        (ReverseDouble i c) = unarily (f 0) dadc i c+  binary f _         (Id dadc) (Lift b)    (ReverseDouble i c) = unarily (f b) dadc i c+  binary f (Id dadb) _         (ReverseDouble i b) Zero        = unarily (`f` 0) dadb i b+  binary f (Id dadb) _         (ReverseDouble i b) (Lift c)    = unarily (`f` c) dadb i b+  binary f (Id dadb) (Id dadc) (ReverseDouble i b) (ReverseDouble j c) = binarily f dadb dadc i b j c++  lift2 f df b c = binary f dadb dadc b c where+    (dadb, dadc) = df (Id (primal b)) (Id (primal c))++  lift2_ f df b c = binary (\_ _ -> a) dadb dadc b c where+    pb = primal b+    pc = primal c+    a = f pb pc+    (dadb, dadc) = df (Id a) (Id pb) (Id pc)+++instance (Reifies s Tape) => Eq (ReverseDouble s) where+  a == b = primal a == primal b++instance (Reifies s Tape) => Ord (ReverseDouble s) where+  compare a b = compare (primal a) (primal b)++instance (Reifies s Tape) => Num (ReverseDouble 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 (Reifies s Tape) => Fractional (ReverseDouble s) where+  fromRational 0 = zero+  fromRational r = auto (fromRational r)+  x / y        = x * recip y+  recip        = lift1_ recip (const . negate . join (*))++instance (Reifies s Tape) => Floating (ReverseDouble 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))+  (**)     = (<**>)+  --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+  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 (Reifies s Tape) => Enum (ReverseDouble s) where+  succ             = lift1 succ (const 1)+  pred             = lift1 pred (const 1)+  toEnum           = auto . toEnum+  fromEnum a       = fromEnum (primal a)+  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 (Reifies s Tape) => Real (ReverseDouble s) where+  toRational      = toRational . primal++instance (Reifies s Tape) => RealFloat (ReverseDouble 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 (Reifies s Tape) => RealFrac (ReverseDouble 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 (Reifies s Tape) => Erf (ReverseDouble 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 -> (recip $ sqrt (2 * pi)) * exp (- x * x / 2)++instance (Reifies s Tape) => InvErf (ReverseDouble s) where+  inverf = lift1_ inverf $ \x _ -> sqrt pi / 2 * exp (x * x)+  inverfc = lift1_ inverfc $ \x _ -> negate (sqrt pi / 2) * exp (x * x)+  invnormcdf = lift1_ invnormcdf $ \x _ -> sqrt (2 * pi) * exp (x * x / 2)++-- | Helper that extracts the derivative of a chain when the chain was constructed with 1 variable.+derivativeOf :: (Reifies s Tape) => Proxy s -> ReverseDouble s -> Double+derivativeOf _ = sum . partials+{-# INLINE derivativeOf #-}++-- | Helper that extracts both the primal and derivative of a chain when the chain was constructed with 1 variable.+derivativeOf' :: (Reifies s Tape) => Proxy s -> ReverseDouble s -> (Double, Double)+derivativeOf' p r = (primal r, derivativeOf p r)+{-# INLINE derivativeOf' #-}++-- | Used internally to push sensitivities down the chain.+backPropagate :: Int -> Cells -> STArray s Int Double -> ST s Int+backPropagate k Nil _ = return k+backPropagate k (Unary i g xs) ss = do+  da <- readArray ss k+  db <- readArray ss i+  writeArray ss i $! db + unsafeCoerce g*da+  (backPropagate $! k - 1) xs ss+backPropagate k (Binary i j g h xs) ss = do+  da <- readArray ss k+  db <- readArray ss i+  writeArray ss i $! db + unsafeCoerce g*da+  dc <- readArray ss j+  writeArray ss j $! dc + unsafeCoerce h*da+  (backPropagate $! k - 1) xs ss++-- | Extract the partials from the current chain for a given AD variable.+partials :: forall s. (Reifies s Tape) => ReverseDouble s -> [Double]+partials Zero        = []+partials (Lift _)    = []+partials (ReverseDouble 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) => Proxy s -> (Int, Int) -> ReverseDouble s -> Array Int Double+partialArrayOf _ vbounds = accumArray (+) 0 vbounds . zip [0..] . partials+{-# INLINE partialArrayOf #-}++-- | Return an 'IntMap' of sparse partials+partialMapOf :: (Reifies s Tape) => Proxy s -> ReverseDouble s-> IntMap Double+partialMapOf _ = fromDistinctAscList . zip [0..] . partials+{-# INLINE partialMapOf #-}++-- | Construct a tape that starts with @n@ variables.+reifyTape :: Int -> (forall s. Reifies s Tape => Proxy s -> r) -> r+reifyTape vs k = unsafePerformIO $ do+  h <- newIORef (Head vs Nil)+  return (reify (Tape h) k)+{-# NOINLINE reifyTape #-}++var :: Double -> Int -> ReverseDouble s+var a v = ReverseDouble v a++varId :: ReverseDouble s -> Int+varId (ReverseDouble v _) = v+varId _ = error "varId: not a Var"++bind :: Traversable f => f Double -> (f (ReverseDouble 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 (ReverseDouble s) -> Array Int Double -> f Double+unbind xs ys = fmap (\v -> ys ! varId v) xs++unbindWith :: (Functor f) => (Double -> b -> c) -> f (ReverseDouble s) -> Array Int b -> f c+unbindWith f xs ys = fmap (\v -> f (primal v) (ys ! varId v)) xs++unbindMap :: (Functor f) => f (ReverseDouble s) -> IntMap Double -> f Double+unbindMap xs ys = fmap (\v -> findWithDefault 0 (varId v) ys) xs++unbindMapWithDefault :: (Functor f) => b -> (Double -> b -> c) -> f (ReverseDouble s) -> IntMap b -> f c+unbindMapWithDefault z f xs ys = fmap (\v -> f (primal v) $ findWithDefault z (varId v) ys) xs

@@ -36,7 +36,7 @@ import Data.Typeable import Control.Comonad.Cofree -infixl 3 :-+infixr 3 :-  -- | A 'Jet' is a tower of all (higher order) partial derivatives of a function --

@@ -150,21 +150,21 @@  -- | 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.+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Reverse.grad' for @n@ inputs, in exchange for better space utilization. grad :: (Traversable f, Num a) => (forall s. f (AD s (Forward a)) -> AD s (Forward a)) -> f a -> f a grad f = Rank1.grad (runAD.f.fmap AD) {-# 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.+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Reverse.grad'' for @n@ inputs, in exchange for better space utilization. grad' :: (Traversable f, Num a) => (forall s. f (AD s (Forward a)) -> AD s (Forward a)) -> f a -> (a, f a) grad' f = Rank1.grad' (runAD.f.fmap AD) {-# 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.+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Reverse.gradWith' for @n@ inputs, in exchange for better space utilization. gradWith :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. f (AD s (Forward a)) -> AD s (Forward a)) -> f a -> f b gradWith g f = Rank1.gradWith g (runAD.f.fmap AD) {-# INLINE gradWith #-}@@ -172,7 +172,7 @@ -- | 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.+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Reverse.gradWith'' for @n@ inputs, in exchange for better space utilization. -- -- >>> gradWith' (,) sum [0..4] -- (10,[(0,1),(1,1),(2,1),(3,1),(4,1)])

@@ -144,21 +144,21 @@  -- | 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.+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Reverse.grad' for @n@ inputs, in exchange for better space utilization. grad :: Traversable f => (forall s. f (AD s ForwardDouble) -> AD s ForwardDouble) -> f Double -> f Double grad f = Rank1.grad (runAD.f.fmap AD) {-# 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.+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Reverse.grad'' for @n@ inputs, in exchange for better space utilization. grad' :: Traversable f => (forall s. f (AD s ForwardDouble) -> AD s ForwardDouble) -> f Double -> (Double, f Double) grad' f = Rank1.grad' (runAD.f.fmap AD) {-# 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.+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Reverse.gradWith' for @n@ inputs, in exchange for better space utilization. gradWith :: Traversable f => (Double -> Double -> b) -> (forall s. f (AD s ForwardDouble) -> AD s ForwardDouble) -> f Double -> f b gradWith g f = Rank1.gradWith g (runAD.f.fmap AD) {-# INLINE gradWith #-}@@ -166,7 +166,7 @@ -- | 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.+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Reverse.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)])

@@ -0,0 +1,216 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE ScopedTypeVariables #-}+-----------------------------------------------------------------------------+-- |+-- Copyright   :  (c) Edward Kmett 2010-2015+-- License     :  BSD3+-- Maintainer  :  ekmett@gmail.com+-- Stability   :  experimental+-- Portability :  GHC only+--+-- Reverse-mode automatic differentiation using Wengert lists and+-- Data.Reflection+--+-- This version is specialized to `Double` enabling the entire+-- structure+--+-----------------------------------------------------------------------------++module Numeric.AD.Mode.Reverse.Double+  ( ReverseDouble, auto+  -- * Gradient+  , grad+  , grad'+  , gradWith+  , gradWith'++  -- * Jacobian+  , jacobian+  , jacobian'+  , jacobianWith+  , jacobianWith'++  -- * Hessian+  , hessian+  , hessianF++  -- * Derivatives+  , diff+  , diff'+  , diffF+  , diffF'+  ) where++#if __GLASGOW_HASKELL__ < 710+import Data.Functor ((<$>))+import Data.Traversable (Traversable)+#endif+import Data.Functor.Compose+import Data.Reflection (Reifies)+import Numeric.AD.Internal.On+import qualified Numeric.AD.Internal.Reverse as R+import qualified Numeric.AD.Mode.Reverse as M+import Numeric.AD.Internal.Reverse.Double+import Numeric.AD.Mode++-- $setup+--+-- >>> import Numeric.AD.Internal.Doctest++-- | 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.0,1.0,1.0]+--+-- >>> grad (\[x,y] -> x**y) [0,2]+-- [0.0,NaN]+grad :: (Traversable f) => (forall s. Reifies s Tape => f (ReverseDouble s) -> ReverseDouble s) -> f Double -> f Double+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.0,[2.0,1.0,1.0])+grad' :: (Traversable f) => (forall s. Reifies s Tape => f (ReverseDouble s) -> ReverseDouble s) -> f Double -> (Double, f Double)+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' #-}++-- | @'grad' g f@ function calculates the gradient of a non-scalar-to-scalar function @f@ with reverse-mode AD in a single pass.+-- The gradient is combined element-wise with the argument using the function @g@.+--+-- @+-- 'grad' == 'gradWith' (\_ dx -> dx)+-- 'id' == 'gradWith' 'const'+-- @+gradWith :: (Traversable f) => (Double -> Double -> b) -> (forall s. Reifies s Tape => f (ReverseDouble s) -> ReverseDouble s) -> f Double -> 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 #-}++-- | @'grad'' g f@ calculates the result and gradient of a non-scalar-to-scalar function @f@ with reverse-mode AD in a single pass+-- the gradient is combined element-wise with the argument using the function @g@.+--+-- @+-- 'grad'' == 'gradWith'' (\_ dx -> dx)+-- @+gradWith' :: (Traversable f) => (Double -> Double -> b) -> (forall s. Reifies s Tape => f (ReverseDouble s) -> ReverseDouble s) -> f Double -> (Double, 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.0,1.0],[1.0,0.0],[1.0,2.0]]+jacobian :: (Traversable f, Functor g) => (forall s. Reifies s Tape => f (ReverseDouble s) -> g (ReverseDouble s)) -> f Double -> g (f Double)+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,+-- where @m@ is the output dimensionality. Applying @fmap snd@ to the result will recover the result of 'jacobian'+-- | An alias for 'gradF''+--+-- >>> jacobian' (\[x,y] -> [y,x,x*y]) [2,1]+-- [(1.0,[0.0,1.0]),(2.0,[1.0,0.0]),(2.0,[1.0,2.0])]+jacobian' :: (Traversable f, Functor g) => (forall s. Reifies s Tape => f (ReverseDouble s) -> g (ReverseDouble s)) -> f Double -> g (Double, f Double)+jacobian' f as = reifyTape (snd bds) $ \p ->+  let row a = (primal a, unbind vs $! partialArrayOf p bds $! a)+  in row <$> f vs+  where (vs, bds) = bind as+{-# INLINE jacobian' #-}++-- | 'jacobianWith g f' calculates the Jacobian of a non-scalar-to-non-scalar function @f@ with reverse AD lazily in @m@ passes for @m@ outputs.+--+-- Instead of returning the Jacobian matrix, the elements of the matrix are combined with the input using the @g@.+--+-- @+-- 'jacobian' == 'jacobianWith' (\_ dx -> dx)+-- 'jacobianWith' 'const' == (\f x -> 'const' x '<$>' f x)+-- @+jacobianWith :: (Traversable f, Functor g) => (Double -> Double -> b) -> (forall s. Reifies s Tape => f (ReverseDouble s) -> g (ReverseDouble s)) -> f Double -> 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,+-- 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) => (Double -> Double -> b) -> (forall s. Reifies s Tape => f (ReverseDouble s) -> g (ReverseDouble s)) -> f Double -> g (Double, 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 vs+  where (vs, bds) = bind as+{-# INLINE jacobianWith' #-}++-- | Compute the derivative of a function.+--+-- >>> diff sin 0+-- 1.0+diff :: (forall s. Reifies s Tape => ReverseDouble s -> ReverseDouble s) -> Double -> Double+diff f a = reifyTape 1 $ \p -> derivativeOf p $! f (var a 0)+{-# INLINE diff #-}++-- | The 'diff'' function calculates the result and derivative, as a pair, of a scalar-to-scalar function.+--+-- >>> diff' sin 0+-- (0.0,1.0)+--+-- >>> diff' exp 0+-- (1.0,1.0)+diff' :: (forall s. Reifies s Tape => ReverseDouble s -> ReverseDouble s) -> Double -> (Double, Double)+diff' f a = reifyTape 1 $ \p -> derivativeOf' p $! f (var a 0)+{-# INLINE diff' #-}++-- | Compute the derivatives of each result of a scalar-to-vector function with regards to its input.+--+-- >>> diffF (\a -> [sin a, cos a]) 0+-- [1.0,0.0]+--+diffF :: (Functor f) => (forall s. Reifies s Tape => ReverseDouble s -> f (ReverseDouble s)) -> Double -> f Double+diffF f a = reifyTape 1 $ \p -> derivativeOf p <$> f (var a 0)+{-# INLINE diffF #-}++-- | Compute the derivatives of each result of a scalar-to-vector function with regards to its input along with the answer.+--+-- >>> diffF' (\a -> [sin a, cos a]) 0+-- [(0.0,1.0),(1.0,0.0)]+diffF' :: (Functor f) => (forall s. Reifies s Tape => ReverseDouble s -> f (ReverseDouble s)) -> Double -> f (Double, Double)+diffF' f a = reifyTape 1 $ \p -> derivativeOf' p <$> f (var a 0)+{-# INLINE diffF' #-}++-- | Compute the hessian via the jacobian of the gradient. gradient is computed in reverse mode and then the jacobian is computed in reverse mode.+--+-- However, since the @'grad' f :: f a -> f a@ is square this is not as fast as using the forward-mode Jacobian of a reverse mode gradient provided by 'Numeric.AD.hessian'.+--+-- >>> hessian (\[x,y] -> x*y) [1,2]+-- [[0.0,1.0],[1.0,0.0]]+hessian :: (Traversable f) => (forall s s'. (Reifies s R.Tape, Reifies s' Tape) => f (On (R.Reverse s (ReverseDouble s'))) -> (On (R.Reverse s (ReverseDouble s')))) -> f Double -> f (f Double)+hessian f = jacobian (M.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.+--+-- Less efficient than 'Numeric.AD.Mode.Mixed.hessianF'.+--+-- >>> hessianF (\[x,y] -> [x*y,x+y,exp x*cos y]) [1,2 :: Double]+-- [[[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) => (forall s s'. (Reifies s R.Tape, Reifies s' Tape) => f (On (R.Reverse s (ReverseDouble s'))) -> g (On (R.Reverse s (ReverseDouble s')))) -> f Double -> g (f (f Double))+hessianF f = getCompose . jacobian (Compose . M.jacobian (fmap off . f . fmap On))+{-# INLINE hessianF #-}

@@ -160,14 +160,14 @@  -- | 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.+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Reverse.grad' for @n@ inputs, in exchange for better space utilization. grad :: (Traversable f, Num a) => (f (Forward a) -> Forward a) -> 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.+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Reverse.grad'' for @n@ inputs, in exchange for better space utilization. grad' :: (Traversable f, Num a) => (f (Forward a) -> Forward a) -> f a -> (a, f a) grad' f as = (primal b, tangent <$> bs) where   (b, bs) = bind' f as@@ -175,7 +175,7 @@  -- | 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.+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Reverse.gradWith' for @n@ inputs, in exchange for better space utilization. gradWith :: (Traversable f, Num a) => (a -> a -> b) -> (f (Forward a) -> Forward a) -> f a -> f b gradWith g f = bindWith g (tangent . f) {-# INLINE gradWith #-}@@ -183,7 +183,7 @@ -- | 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.+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Reverse.gradWith'' for @n@ inputs, in exchange for better space utilization. -- -- >>> gradWith' (,) sum [0..4] -- (10,[(0,1),(1,1),(2,1),(3,1),(4,1)])

@@ -139,14 +139,14 @@  -- | 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.+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Reverse.grad' for @n@ inputs, in exchange for better space utilization. grad :: Traversable f => (f ForwardDouble -> ForwardDouble) -> 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.+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Reverse.grad'' for @n@ inputs, in exchange for better space utilization. grad' :: Traversable f => (f ForwardDouble -> ForwardDouble) -> f Double -> (Double, f Double) grad' f as = (primal b, tangent <$> bs)     where@@ -155,7 +155,7 @@  -- | 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.+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Reverse.gradWith' for @n@ inputs, in exchange for better space utilization. gradWith :: Traversable f => (Double -> Double -> b) -> (f ForwardDouble -> ForwardDouble) -> f Double -> f b gradWith g f = bindWith g (tangent . f) {-# INLINE gradWith #-}@@ -163,7 +163,7 @@ -- | 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.+-- Note, this performs /O(n)/ worse than 'Numeric.AD.Mode.Reverse.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)])