ad 4.3.6 → 4.5.6
raw patch · 72 files changed
Files
- .gitignore +1/−0
- .hlint.yaml +15/−0
- .travis.yml +0/−151
- CHANGELOG.markdown +90/−0
- README.markdown +1/−1
- Setup.lhs +5/−32
- Warning.hs +0/−5
- ad.cabal +66/−44
- cbits/tape.c +128/−0
- include/instances.h +28/−20
- include/internal_kahn.h +327/−0
- include/rank1_kahn.h +270/−0
- src/Numeric/AD.hs +43/−11
- src/Numeric/AD/Double.hs +246/−0
- src/Numeric/AD/Halley.hs +0/−118
- src/Numeric/AD/Halley/Double.hs +83/−0
- src/Numeric/AD/Internal/Combinators.hs +5/−8
- src/Numeric/AD/Internal/Dense.hs +13/−19
- src/Numeric/AD/Internal/Dense/Representable.hs +178/−0
- src/Numeric/AD/Internal/Doctest.hs +1/−1
- src/Numeric/AD/Internal/Forward.hs +7/−19
- src/Numeric/AD/Internal/Forward/Double.hs +13/−101
- src/Numeric/AD/Internal/Identity.hs +6/−5
- src/Numeric/AD/Internal/Kahn.hs +12/−17
- src/Numeric/AD/Internal/Kahn/Double.hs +9/−0
- src/Numeric/AD/Internal/Kahn/Float.hs +9/−0
- src/Numeric/AD/Internal/On.hs +5/−6
- src/Numeric/AD/Internal/Or.hs +4/−33
- src/Numeric/AD/Internal/Reverse.hs +18/−24
- src/Numeric/AD/Internal/Reverse/Double.hs +335/−0
- src/Numeric/AD/Internal/Sparse.hs +15/−59
- src/Numeric/AD/Internal/Sparse/Common.hs +54/−0
- src/Numeric/AD/Internal/Sparse/Double.hs +258/−0
- src/Numeric/AD/Internal/Tower.hs +11/−12
- src/Numeric/AD/Internal/Tower/Double.hs +271/−0
- src/Numeric/AD/Internal/Type.hs +2/−1
- src/Numeric/AD/Jacobian.hs +1/−2
- src/Numeric/AD/Jet.hs +9/−32
- src/Numeric/AD/Mode.hs +38/−34
- src/Numeric/AD/Mode/Dense.hs +68/−0
- src/Numeric/AD/Mode/Dense/Representable.hs +79/−0
- src/Numeric/AD/Mode/Forward.hs +5/−9
- src/Numeric/AD/Mode/Forward/Double.hs +5/−9
- src/Numeric/AD/Mode/Kahn.hs +1/−6
- src/Numeric/AD/Mode/Kahn/Double.hs +240/−0
- src/Numeric/AD/Mode/Reverse.hs +95/−34
- src/Numeric/AD/Mode/Reverse/Double.hs +274/−0
- src/Numeric/AD/Mode/Sparse.hs +1/−5
- src/Numeric/AD/Mode/Sparse/Double.hs +162/−0
- src/Numeric/AD/Mode/Tower.hs +1/−2
- src/Numeric/AD/Mode/Tower/Double.hs +115/−0
- src/Numeric/AD/Newton.hs +3/−8
- src/Numeric/AD/Newton/Double.hs +4/−32
- src/Numeric/AD/Rank1/Dense.hs +102/−0
- src/Numeric/AD/Rank1/Dense/Representable.hs +103/−0
- src/Numeric/AD/Rank1/Forward.hs +109/−34
- src/Numeric/AD/Rank1/Forward/Double.hs +107/−28
- src/Numeric/AD/Rank1/Halley.hs +4/−2
- src/Numeric/AD/Rank1/Halley/Double.hs +122/−0
- src/Numeric/AD/Rank1/Kahn.hs +13/−216
- src/Numeric/AD/Rank1/Kahn/Double.hs +22/−0
- src/Numeric/AD/Rank1/Kahn/Float.hs +22/−0
- src/Numeric/AD/Rank1/Newton.hs +1/−5
- src/Numeric/AD/Rank1/Newton/Double.hs +2/−8
- src/Numeric/AD/Rank1/Sparse.hs +1/−5
- src/Numeric/AD/Rank1/Sparse/Double.hs +192/−0
- src/Numeric/AD/Rank1/Tower.hs +93/−24
- src/Numeric/AD/Rank1/Tower/Double.hs +199/−0
- tests/Regression.hs +158/−0
- tests/doctests.hs +0/−25
- travis/cabal-apt-install +0/−27
- travis/config +0/−16
.gitignore view
@@ -31,3 +31,4 @@ cabal.project.local~ .HTF/ .ghc.environment.*+Makefile
+ .hlint.yaml view
@@ -0,0 +1,15 @@+- arguments: [-XCPP, --cpp-define=HLINT, --cpp-ansi, --cpp-include=include]++- fixity: "infixr 8 **!, <**>"+- fixity: "infixl 7 *!, /!, ^*, *^, ^/"+- fixity: "infixl 6 +!, -!, <+>"+- fixity: "infixr 5 :!"+- fixity: "infix 4 ==!"+- fixity: "infixl 3 :-"++# this doesn't work well with Rank2Types+- ignore: {name: Eta reduce}+# Numeric.AD.Rank1.Kahn's use of CPP makes it difficult to avoid redundant parentheses+- ignore: {name: Redundant bracket, within: [Numeric.AD.Rank1.Kahn]}+- ignore: {name: Unused LANGUAGE pragma}+- ignore: {name: Reduce duplication}
− .travis.yml
@@ -1,151 +0,0 @@-# This Travis job script has been generated by a script via-#-# runghc make_travis_yml_2.hs '-o' '.travis.yml' '--ghc-head' '--irc-channel=irc.freenode.org#haskell-lens' '--no-no-tests-no-bench' '--no-unconstrained' 'cabal.project'-#-# For more information, see https://github.com/haskell-CI/haskell-ci-#-language: c-sudo: false--git:- submodules: false # whether to recursively clone submodules--notifications:- irc:- channels:- - "irc.freenode.org#haskell-lens"- skip_join: true- template:- - "\x0313ad\x03/\x0306%{branch}\x03 \x0314%{commit}\x03 %{build_url} %{message}"--cache:- directories:- - $HOME/.cabal/packages- - $HOME/.cabal/store--before_cache:- - rm -fv $HOME/.cabal/packages/hackage.haskell.org/build-reports.log- # remove files that are regenerated by 'cabal update'- - rm -fv $HOME/.cabal/packages/hackage.haskell.org/00-index.*- - rm -fv $HOME/.cabal/packages/hackage.haskell.org/*.json- - rm -fv $HOME/.cabal/packages/hackage.haskell.org/01-index.cache- - rm -fv $HOME/.cabal/packages/hackage.haskell.org/01-index.tar- - rm -fv $HOME/.cabal/packages/hackage.haskell.org/01-index.tar.idx-- - rm -rfv $HOME/.cabal/packages/head.hackage--matrix:- include:- - compiler: "ghc-8.6.3"- # env: TEST=--disable-tests BENCH=--disable-benchmarks- addons: {apt: {packages: [ghc-ppa-tools,cabal-install-2.4,ghc-8.6.3], sources: [hvr-ghc]}}- - compiler: "ghc-8.4.4"- # env: TEST=--disable-tests BENCH=--disable-benchmarks- addons: {apt: {packages: [ghc-ppa-tools,cabal-install-2.4,ghc-8.4.4], sources: [hvr-ghc]}}- - compiler: "ghc-8.2.2"- # env: TEST=--disable-tests BENCH=--disable-benchmarks- addons: {apt: {packages: [ghc-ppa-tools,cabal-install-2.4,ghc-8.2.2], sources: [hvr-ghc]}}- - compiler: "ghc-8.0.2"- # env: TEST=--disable-tests BENCH=--disable-benchmarks- addons: {apt: {packages: [ghc-ppa-tools,cabal-install-2.4,ghc-8.0.2], sources: [hvr-ghc]}}- - compiler: "ghc-7.10.3"- # env: TEST=--disable-tests BENCH=--disable-benchmarks- addons: {apt: {packages: [ghc-ppa-tools,cabal-install-2.4,ghc-7.10.3], sources: [hvr-ghc]}}- - compiler: "ghc-7.8.4"- # env: TEST=--disable-tests BENCH=--disable-benchmarks- addons: {apt: {packages: [ghc-ppa-tools,cabal-install-2.4,ghc-7.8.4], sources: [hvr-ghc]}}- - compiler: "ghc-7.6.3"- # env: TEST=--disable-tests BENCH=--disable-benchmarks- addons: {apt: {packages: [ghc-ppa-tools,cabal-install-2.4,ghc-7.6.3], sources: [hvr-ghc]}}- - compiler: "ghc-7.4.2"- # env: TEST=--disable-tests BENCH=--disable-benchmarks- addons: {apt: {packages: [ghc-ppa-tools,cabal-install-2.4,ghc-7.4.2], sources: [hvr-ghc]}}- - compiler: "ghc-head"- env: GHCHEAD=true- addons: {apt: {packages: [ghc-ppa-tools,cabal-install-head,ghc-head], sources: [hvr-ghc]}}-- allow_failures:- - compiler: "ghc-head"--before_install:- - HC=${CC}- - HCPKG=${HC/ghc/ghc-pkg}- - unset CC- - ROOTDIR=$(pwd)- - mkdir -p $HOME/.local/bin- - "PATH=/opt/ghc/bin:/opt/ghc-ppa-tools/bin:$HOME/local/bin:$PATH"- - HCNUMVER=$(( $(${HC} --numeric-version|sed -E 's/([0-9]+)\.([0-9]+)\.([0-9]+).*/\1 * 10000 + \2 * 100 + \3/') ))- - echo $HCNUMVER--install:- - cabal --version- - echo "$(${HC} --version) [$(${HC} --print-project-git-commit-id 2> /dev/null || echo '?')]"- - BENCH=${BENCH---enable-benchmarks}- - TEST=${TEST---enable-tests}- - HADDOCK=${HADDOCK-true}- - UNCONSTRAINED=${UNCONSTRAINED-true}- - NOINSTALLEDCONSTRAINTS=${NOINSTALLEDCONSTRAINTS-false}- - GHCHEAD=${GHCHEAD-false}- - travis_retry cabal update -v- - "sed -i.bak 's/^jobs:/-- jobs:/' ${HOME}/.cabal/config"- - rm -fv cabal.project cabal.project.local- # Overlay Hackage Package Index for GHC HEAD: https://github.com/hvr/head.hackage- - |- if $GHCHEAD; then- sed -i 's/-- allow-newer: .*/allow-newer: *:base/' ${HOME}/.cabal/config- for pkg in $($HCPKG list --simple-output); do pkg=$(echo $pkg | sed 's/-[^-]*$//'); sed -i "s/allow-newer: /allow-newer: *:$pkg, /" ${HOME}/.cabal/config; done-- echo 'repository head.hackage' >> ${HOME}/.cabal/config- echo ' url: http://head.hackage.haskell.org/' >> ${HOME}/.cabal/config- echo ' secure: True' >> ${HOME}/.cabal/config- echo ' root-keys: 07c59cb65787dedfaef5bd5f987ceb5f7e5ebf88b904bbd4c5cbdeb2ff71b740' >> ${HOME}/.cabal/config- echo ' 2e8555dde16ebd8df076f1a8ef13b8f14c66bad8eafefd7d9e37d0ed711821fb' >> ${HOME}/.cabal/config- echo ' 8f79fd2389ab2967354407ec852cbe73f2e8635793ac446d09461ffb99527f6e' >> ${HOME}/.cabal/config- echo ' key-threshold: 3' >> ${HOME}/.cabal.config-- grep -Ev -- '^\s*--' ${HOME}/.cabal/config | grep -Ev '^\s*$'-- cabal new-update head.hackage -v- fi- - grep -Ev -- '^\s*--' ${HOME}/.cabal/config | grep -Ev '^\s*$'- - "printf 'packages: \".\"\\n' > cabal.project"- - "printf 'write-ghc-environment-files: always\\n' >> cabal.project"- - touch cabal.project.local- - "if ! $NOINSTALLEDCONSTRAINTS; then for pkg in $($HCPKG list --simple-output); do echo $pkg | grep -vw -- ad | sed 's/^/constraints: /' | sed 's/-[^-]*$/ installed/' >> cabal.project.local; done; fi"- - cat cabal.project || true- - cat cabal.project.local || true- - if [ -f "./configure.ac" ]; then- (cd "." && autoreconf -i);- fi- - rm -f cabal.project.freeze- - cabal new-build -w ${HC} ${TEST} ${BENCH} --project-file="cabal.project" --dep -j2 all- - rm -rf .ghc.environment.* "."/dist- - DISTDIR=$(mktemp -d /tmp/dist-test.XXXX)--# Here starts the actual work to be performed for the package under test;-# any command which exits with a non-zero exit code causes the build to fail.-script:- # test that source-distributions can be generated- - cabal new-sdist all- - mv dist-newstyle/sdist/*.tar.gz ${DISTDIR}/- - cd ${DISTDIR} || false- - find . -maxdepth 1 -name '*.tar.gz' -exec tar -xvf '{}' \;- - "printf 'packages: ad-*/*.cabal\\n' > cabal.project"- - "printf 'write-ghc-environment-files: always\\n' >> cabal.project"- - touch cabal.project.local- - "if ! $NOINSTALLEDCONSTRAINTS; then for pkg in $($HCPKG list --simple-output); do echo $pkg | grep -vw -- ad | sed 's/^/constraints: /' | sed 's/-[^-]*$/ installed/' >> cabal.project.local; done; fi"- - cat cabal.project || true- - cat cabal.project.local || true-- # build & run tests, build benchmarks- - cabal new-build -w ${HC} ${TEST} ${BENCH} all- - if [ "x$TEST" = "x--enable-tests" ]; then cabal new-test -w ${HC} ${TEST} ${BENCH} all; fi-- # cabal check- - (cd ad-* && cabal check)-- # haddock- - if $HADDOCK; then cabal new-haddock -w ${HC} ${TEST} ${BENCH} all; else echo "Skipping haddock generation";fi--# REGENDATA ["-o",".travis.yml","--ghc-head","--irc-channel=irc.freenode.org#haskell-lens","--no-no-tests-no-bench","--no-unconstrained","cabal.project"]-# EOF
CHANGELOG.markdown view
@@ -1,3 +1,93 @@+4.5.6 [2024.05.01]+------------------+* Add specialized implementations of `log1p`, `expm1`, `log1pexp`, and+ `log1mexp` in `Floating` instances.++4.5.5 [2024.01.28]+------------------+* `Numeric.AD.Mode.Reverse.Double` now handles IEEE floating-point special+ values (e.g., `NaN` and `Inf`) correctly when `ad` is compiled with `+ffi`.+ Note that this increase in floating-point accuracy may come at a slight+ performance penalty in certain applications. If this negatively impacts your+ application, please mention this at https://github.com/ekmett/ad/issues/106.++4.5.4 [2023.02.19]+------------------+* Add a `Num (Scalar (Scalar t))` constraint to `On`'s `Mode` instance, which is+ required to make it typecheck with GHC 9.6.++ (Note that this constraint was already present implicitly due to superclass+ expansion, so this is not a breaking change. The only reason that it must be+ added explicitly with GHC 9.6 or later is due to 9.6 being more conservative+ with superclass expansion.)++4.5.3 [2023.01.21]+------------------+* Support building with GHC 9.6.++4.5.2 [2022.06.17]+------------------+* Fix a bug that would cause `Numeric.AD.Mode.Reverse.diff` and+ `Numeric.AD.Mode.Reverse.Double.diff` to compute different answers under+ certain circumstances when `ad` was compiled with the `+ffi` flag.++4.5.1 [2022.05.18]+------------------+* Allow building with `transformers-0.6.*`.++4.5 [2021.11.07]+----------------+* The build-type has been changed from `Custom` to `Simple`.+ To achieve this, the `doctests` test suite has been removed in favor of using+ [`cabal-docspec`](https://github.com/phadej/cabal-extras/tree/master/cabal-docspec)+ to run the doctests.+* Expose `Dense` mode AD again.+* Add a `Dense.Representable` mode, which is a variant of `Dense` that exploits+ `Representable` functors rather than `Traversable` functors.+* `Representable` can now also be useful as it can allow us to `unjet` to convert+ a value of type `Jet f a` safely back into `Cofree f a`.+* Improve `Reverse.Double` mode performance by increasing strictness and using an FFI-based tape.+* Reverse mode AD uses `reifyTypeable` internally. This means the region parameter/infinitesimals+ that mark each tape are `Typeable`, allowing you to do things like define instances of `Exception`+ that name the region parameter and perform similar shenanigans.+* Drastically reduce code duplication in `Double`-based modes, enabling more of them.+* Fixed a number of modes that were handling `(**)` improperly due to the aforementioned code+ duplication problem.+* Add a `Tower.Double` mode (internally) that uses lazy lists of strict doubles.+* Add a `Kahn.Double` mode (internally) that holds strict doubles in the graph.+* Switch to using pattern synonyms internally for detecting "known" zeros.+* Drop support for versions of GHC before 8.0+* The `.Double` modes have been modified to exploit the fact that we can definitely check a Double for equality with 0.+ In future releases we may require a typeclass that offers the ability to check for known zeroes for all types you+ process. This will allow us to improve the quality of the results, but may require you to either write an small instance+ declaration if you are processing some esoteric data type of your own, or put on/off a newtype that indicates to skip+ known zero optimizations or to use Eq. If there are particularly common types with tricky cases, a future `ad-instances`+ package might be the right way forward for them to find a home.+* Add `Numeric.AD.Double`, which tries to mix and match between all the different AD modes to produce optimal results+ but uses the various `.Double` specializations to reduce the amount of boxing and indirection on the heap.+* Add `Numeric.AD.Halley.Double`.+* Removed the `fooNoEq` variants from `Newton.Double`, `Double`s always have an `Eq` instance.++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`:++ ```diff+ -stochasticGradientDescent :: (Traversable f, Fractional a, Ord a) => (forall s. Reifies s Tape => f (Scalar a) -> f (Reverse s a) -> Reverse s a) -> [f (Scalar a)] -> f a -> [f a]+ +stochasticGradientDescent :: (Traversable f, Fractional a, Ord a) => (forall s. Reifies s Tape => e -> f (Reverse s a) -> Reverse s a) -> [e] -> f a -> [f a]+ ```+ 4.3.6 [2019.02.28] ------------------ * Make the test suite pass when built against `musl` `libc`.
README.markdown view
@@ -1,7 +1,7 @@ ad == -[](https://hackage.haskell.org/package/ad) [](http://travis-ci.org/ekmett/ad)+[](https://hackage.haskell.org/package/ad) [](https://github.com/ekmett/ad/actions?query=workflow%3AHaskell-CI) A package that provides an intuitive API for [Automatic Differentiation](http://en.wikipedia.org/wiki/Automatic_differentiation) (AD) in Haskell. Automatic differentiation provides a means to calculate the derivatives of a function while evaluating it. Unlike numerical methods based on running the program with multiple inputs or symbolic approaches, automatic differentiation typically only decreases performance by a small multiplier.
Setup.lhs view
@@ -1,34 +1,7 @@-\begin{code}-{-# LANGUAGE CPP #-}-{-# OPTIONS_GHC -Wall #-}-module Main (main) where--#ifndef MIN_VERSION_cabal_doctest-#define MIN_VERSION_cabal_doctest(x,y,z) 0-#endif--#if MIN_VERSION_cabal_doctest(1,0,0)--import Distribution.Extra.Doctest ( defaultMainWithDoctests )-main :: IO ()-main = defaultMainWithDoctests "doctests"--#else--#ifdef MIN_VERSION_Cabal--- If the macro is defined, we have new cabal-install,--- but for some reason we don't have cabal-doctest in package-db------ Probably we are running cabal sdist, when otherwise using new-build--- workflow-import Warning ()-#endif--import Distribution.Simple--main :: IO ()-main = defaultMain+#!/usr/bin/runhaskell+> module Main (main) where -#endif+> import Distribution.Simple -\end{code}+> main :: IO ()+> main = defaultMain
− Warning.hs
@@ -1,5 +0,0 @@-module Warning- {-# WARNING ["You are configuring this package without cabal-doctest installed.",- "The doctests test-suite will not work as a result.",- "To fix this, install cabal-doctest before configuring."] #-}- () where
ad.cabal view
@@ -1,8 +1,8 @@ name: ad-version: 4.3.6+version: 4.5.6 license: BSD3 license-File: LICENSE-copyright: (c) Edward Kmett 2010-2015,+copyright: (c) Edward Kmett 2010-2021, (c) Barak Pearlmutter and Jeffrey Mark Siskind 2008-2009 author: Edward Kmett maintainer: ekmett@gmail.com@@ -10,27 +10,28 @@ category: Math homepage: http://github.com/ekmett/ad bug-reports: http://github.com/ekmett/ad/issues-build-type: Custom+build-type: Simple cabal-version: >= 1.10-tested-with: GHC == 7.4.2- , GHC == 7.6.3- , GHC == 7.8.4- , GHC == 7.10.3- , GHC == 8.0.2+tested-with: GHC == 8.0.2 , GHC == 8.2.2 , GHC == 8.4.4- , GHC == 8.6.3+ , GHC == 8.6.5+ , GHC == 8.8.4+ , GHC == 8.10.7+ , GHC == 9.0.2+ , GHC == 9.2.7+ , GHC == 9.4.5+ , GHC == 9.6.2 synopsis: Automatic Differentiation extra-source-files: .gitignore- .travis.yml+ .hlint.yaml .vim.custom CHANGELOG.markdown README.markdown- Warning.hs- travis/cabal-apt-install- travis/config include/instances.h+ include/rank1_kahn.h+ include/internal_kahn.h description: Forward-, reverse- and mixed- mode automatic differentiation combinators with a common API. .@@ -88,17 +89,14 @@ default: False manual: True -custom-setup- setup-depends:- base >= 4.3 && <5,- Cabal >= 1.10,- cabal-doctest >= 1 && <1.1+flag ffi+ default: False+ manual: True library hs-source-dirs: src include-dirs: include default-language: Haskell2010- other-extensions: BangPatterns DeriveDataTypeable@@ -108,25 +106,26 @@ GeneralizedNewtypeDeriving MultiParamTypeClasses PatternGuards+ PatternSynonyms Rank2Types ScopedTypeVariables- TemplateHaskell TypeFamilies TypeOperators UndecidableInstances build-depends:- array >= 0.2 && < 0.6,- base >= 4.3 && < 5,- comonad >= 4 && < 6,- containers >= 0.2 && < 0.7,- data-reify >= 0.6 && < 0.7,- erf >= 2.0 && < 2.1,- free >= 4.6.1 && < 6,- nats >= 0.1.2 && < 2,- reflection >= 1.4 && < 3,- semigroups >= 0.16 && < 1,- transformers >= 0.3 && < 0.6+ adjunctions >= 4.4 && < 5,+ array >= 0.4 && < 0.6,+ base >= 4.9 && < 5,+ comonad >= 4 && < 6,+ containers >= 0.5 && < 0.8,+ data-reify >= 0.6 && < 0.7,+ erf >= 2.0 && < 2.1,+ free >= 4.6.1 && < 6,+ nats >= 0.1.2 && < 2,+ reflection >= 1.4 && < 3,+ semigroups >= 0.16 && < 1,+ transformers >= 0.5.2.0 && < 0.7 if impl(ghc < 7.8) build-depends: tagged >= 0.7 && < 1@@ -138,59 +137,82 @@ exposed-modules: Numeric.AD- Numeric.AD.Halley+ Numeric.AD.Double+ Numeric.AD.Halley.Double Numeric.AD.Internal.Dense+ Numeric.AD.Internal.Dense.Representable+ Numeric.AD.Internal.Doctest Numeric.AD.Internal.Forward Numeric.AD.Internal.Forward.Double Numeric.AD.Internal.Identity Numeric.AD.Internal.Kahn+ Numeric.AD.Internal.Kahn.Double+ Numeric.AD.Internal.Kahn.Float Numeric.AD.Internal.On Numeric.AD.Internal.Or Numeric.AD.Internal.Reverse+ Numeric.AD.Internal.Reverse.Double Numeric.AD.Internal.Sparse+ Numeric.AD.Internal.Sparse.Common+ Numeric.AD.Internal.Sparse.Double Numeric.AD.Internal.Tower+ Numeric.AD.Internal.Tower.Double Numeric.AD.Internal.Type Numeric.AD.Jacobian Numeric.AD.Jet Numeric.AD.Mode+ Numeric.AD.Mode.Dense+ Numeric.AD.Mode.Dense.Representable Numeric.AD.Mode.Forward Numeric.AD.Mode.Forward.Double Numeric.AD.Mode.Kahn+ Numeric.AD.Mode.Kahn.Double Numeric.AD.Mode.Reverse+ Numeric.AD.Mode.Reverse.Double Numeric.AD.Mode.Sparse+ Numeric.AD.Mode.Sparse.Double Numeric.AD.Mode.Tower+ Numeric.AD.Mode.Tower.Double Numeric.AD.Newton Numeric.AD.Newton.Double+ Numeric.AD.Rank1.Dense+ Numeric.AD.Rank1.Dense.Representable Numeric.AD.Rank1.Forward Numeric.AD.Rank1.Forward.Double Numeric.AD.Rank1.Halley+ Numeric.AD.Rank1.Halley.Double Numeric.AD.Rank1.Kahn+ Numeric.AD.Rank1.Kahn.Double+ Numeric.AD.Rank1.Kahn.Float Numeric.AD.Rank1.Newton Numeric.AD.Rank1.Newton.Double Numeric.AD.Rank1.Sparse+ Numeric.AD.Rank1.Sparse.Double Numeric.AD.Rank1.Tower+ Numeric.AD.Rank1.Tower.Double + if flag(ffi)+ other-extensions: ForeignFunctionInterface+ c-sources: cbits/tape.c+ cpp-options: -DAD_FFI+ other-modules:- Numeric.AD.Internal.Doctest 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+ x-docspec-extra-packages: distributive --- Verify the results of the examples-test-suite doctests+test-suite regression default-language: Haskell2010- type: exitcode-stdio-1.0- main-is: doctests.hs- build-depends:- ad,- base,- directory,- doctest >= 0.9.0.1 && < 0.17,- filepath- ghc-options: -Wall -threaded+ type: exitcode-stdio-1.0+ main-is: Regression.hs hs-source-dirs: tests+ build-depends: ad, base, tasty, tasty-hunit+ ghc-options: -fspec-constr -fdicts-cheap -O2 benchmark blackscholes default-language: Haskell2010
+ cbits/tape.c view
@@ -0,0 +1,128 @@+#include <stdlib.h>+#include <math.h>+#include <string.h>+#include <stdio.h>++typedef struct Tape+{+ int idx;+ int offset;+ int size;+ int variables;+ double *val;+ int *lnk;+ struct Tape* prev; +} tape_t;++int tape_variables(void *p)+{+ tape_t *pTape = (tape_t*)p;+ return pTape->variables;+}++void* tape_alloc(int variables, int size)+{+ void* p = malloc(sizeof(tape_t) + size*2*(sizeof(double) + sizeof(int)) );+ tape_t *pTape = (tape_t*)p;++ pTape->size = size;+ pTape->idx = 0;+ pTape->offset = variables;+ pTape->variables = variables;++ pTape->val = p + sizeof(tape_t);+ pTape->lnk = p + sizeof(tape_t) + size*2*sizeof(double);+ pTape->prev = 0;++ return pTape;+}++int tape_push(void* p, int i_l, int i_r, double d_l, double d_r)+{+ tape_t *pTape = (tape_t*)p;++ // time to allocate new block?+ if (pTape->idx >= pTape->size)+ {+ int newSize = pTape->size * 2;++ p = malloc( sizeof(tape_t) + newSize*2*(sizeof(double) + sizeof(int)) );+ + tape_t *pNew = (tape_t*)p;+ *pNew = *pTape;++ pTape->idx = 0;+ pTape->val = p + sizeof(tape_t);+ pTape->lnk = p + sizeof(tape_t) + newSize*2*sizeof(double);+ pTape->offset = pNew->offset + pNew->size;+ pTape->size = newSize;+ pTape->prev = pNew;+ }++ int i = pTape->idx++;++ pTape->val[i*2] = d_l;+ pTape->val[i*2 + 1] = d_r;++ pTape->lnk[i*2] = i_l;+ pTape->lnk[i*2 + 1] = i_r;++ return (i + pTape->offset);+}++void tape_backPropagate(void* p, int start, double* out)+{+ tape_t *pTape = (tape_t*)p;++ int variables = pTape->variables;++ double* buffer = calloc( pTape->offset + pTape->idx, sizeof(double) );+ buffer[start] = 1.0;+ + int idx = 1 + start;++ while (pTape)+ {+ idx -= pTape->offset;++ while (--idx >= 0)+ {+ double v = buffer[idx + pTape->offset];++ // TODO: if we do not care about handling IEEE floating point special values (NaN, Inf) correctly+ // then we can skip the rest of the loop body in case v == 0+ // see also https://github.com/ekmett/ad/issues/106++ int i = pTape->lnk[idx*2];+ if (i >= 0)+ {+ double x = v * pTape->val[idx*2];+ if (x != 0) buffer[i] += x;+ }++ int j = pTape->lnk[idx*2 + 1]; + if (j >= 0)+ {+ double y = v * pTape->val[idx*2 + 1];+ if (y != 0) buffer[j] += y;+ }+ }+ idx += 1 + pTape->offset;+ pTape = pTape->prev;+ }+ + memcpy(out, buffer, variables * sizeof(double) );+ free(buffer);+}++void tape_free(void* p)+{+ tape_t *pTape = (tape_t*)p;++ while (pTape)+ {+ p = pTape;+ pTape = pTape->prev;+ free(p);+ }+}
include/instances.h view
@@ -1,22 +1,24 @@ #ifndef BODY1-#define BODY1(x) x+#define BODY1(x) x => #endif #ifndef BODY2-#define BODY2(x,y) (x,y)+#define BODY2(x,y) (x,y) => #endif -instance BODY2(Num a, Eq a) => Eq (HEAD) where+instance BODY2(Num a, Eq a) Eq HEAD where a == b = primal a == primal b -instance BODY2(Num a, Ord a) => Ord (HEAD) where+instance BODY2(Num a, Ord a) Ord HEAD where compare a b = compare (primal a) (primal b) -instance BODY2(Num a, Bounded a) => Bounded (HEAD) where+#ifndef NO_Bounded+instance BODY2(Num a, Bounded a) Bounded HEAD where maxBound = auto maxBound minBound = auto minBound+#endif -instance BODY1(Num a) => Num (HEAD) where+instance BODY1(Num a) Num HEAD where fromInteger 0 = zero fromInteger n = auto (fromInteger n) (+) = (<+>) -- binary (+) 1 1@@ -26,23 +28,24 @@ abs = lift1 abs signum signum a = lift1 signum (const zero) a -instance BODY1(Fractional a) => Fractional (HEAD) where+instance BODY1(Fractional a) Fractional HEAD where fromRational 0 = zero fromRational r = auto (fromRational r) x / y = x * recip y recip = lift1_ recip (const . negate . join (*)) -instance BODY1(Floating a) => Floating (HEAD) where+instance BODY1(Floating a) Floating HEAD 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++ KnownZero ** y = auto (0 ** primal y)+ _ ** KnownZero = 1+ x ** Auto y = lift1 (**y) (\z -> y *^ z ** auto (y-1)) x+ x ** y = lift2_ (**) (\z xi yi -> (yi * xi ** (yi - 1), z * log xi)) x y+ sin = lift1 sin cos cos = lift1 cos $ negate . sin tan = lift1 tan $ recip . join (*) . cos@@ -56,7 +59,12 @@ acosh = lift1 acosh $ \x -> recip (sqrt (join (*) x - 1)) atanh = lift1 atanh $ \x -> recip (1 - join (*) x) -instance BODY2(Num a, Enum a) => Enum (HEAD) where+ log1p = lift1 log1p $ recip . (+) 1+ expm1 = lift1 expm1 exp+ log1pexp = lift1 log1pexp $ recip . (+) 1 . exp . negate+ log1mexp = lift1 log1mexp $ recip . negate . expm1 . negate++instance BODY2(Num a, Enum a) Enum HEAD where succ = lift1 succ (const 1) pred = lift1 pred (const 1) toEnum = auto . toEnum@@ -66,10 +74,10 @@ 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 BODY1(Real a) => Real (HEAD) where+instance BODY1(Real a) Real HEAD where toRational = toRational . primal -instance BODY1(RealFloat a) => RealFloat (HEAD) where+instance BODY1(RealFloat a) RealFloat HEAD where floatRadix = floatRadix . primal floatDigits = floatDigits . primal floatRange = floatRange . primal@@ -85,7 +93,7 @@ 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 BODY1(RealFrac a) => RealFrac (HEAD) where+instance BODY1(RealFrac a) RealFrac HEAD where properFraction a = (w, a `withPrimal` pb) where pa = primal a (w, pb) = properFraction pa@@ -94,12 +102,12 @@ ceiling = ceiling . primal floor = floor . primal -instance BODY1(Erf a) => Erf (HEAD) where+instance BODY1(Erf a) Erf HEAD 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)+ normcdf = lift1 normcdf $ \x -> recip (sqrt (2 * pi)) * exp (- x * x / 2) -instance BODY1(InvErf a) => InvErf (HEAD) where+instance BODY1(InvErf a) InvErf HEAD 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)
+ include/internal_kahn.h view
@@ -0,0 +1,327 @@++{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE UndecidableInstances #-}+{-# OPTIONS_GHC -fno-full-laziness #-}+{-# OPTIONS_HADDOCK not-home #-}++-----------------------------------------------------------------------------+-- |+-- Copyright : (c) Edward Kmett 2010-2021+-- License : BSD3+-- Maintainer : ekmett@gmail.com+-- Stability : experimental+-- Portability : GHC only+--+-- This module provides reverse-mode Automatic Differentiation implementation using+-- linear time topological sorting after the fact.+--+-- For this form of reverse-mode AD we use 'System.Mem.StableName.StableName' to recover+-- sharing information from the tape to avoid combinatorial explosion, and thus+-- run asymptotically faster than it could without such sharing information, but the use+-- of side-effects contained herein is benign.+--+-----------------------------------------------------------------------------++MODULE+ ( AD_EXPORT+ , Tape(..)+ , partials+ , partialArray+ , partialMap+ , derivative+ , derivative'+ , vgrad, vgrad'+ , Grad(..)+ , bind+ , unbind+ , unbindMap+ , unbindWithUArray+ , unbindWithArray+ , unbindMapWithDefault+ , primal+ , var+ , varId+ ) where++import Control.Monad.ST+import Control.Monad hiding (mapM)+import Control.Monad.Trans.State+import qualified Data.List as List (foldl')+import Data.Array.ST+import Data.Array.IArray+import qualified Data.Array as A+import Data.Array.Unboxed (UArray)+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 Data.Data (Data)+import Data.Typeable (Typeable)+import qualified GHC.Exts as Exts+import Numeric+import Numeric.AD.Internal.Combinators+import Numeric.AD.Internal.Identity+import Numeric.AD.Jacobian+import Numeric.AD.Mode+IMPORTS++-- | A @Tape@ records the information needed back propagate from the output to each input during reverse 'Mode' AD.+data Tape t+ = Zero+ | Lift {-# UNPACK #-} !SCALAR_TYPE+ | Var {-# UNPACK #-} !SCALAR_TYPE {-# UNPACK #-} !Int+ | Binary {-# UNPACK #-} !SCALAR_TYPE {-# UNPACK #-} !SCALAR_TYPE {-# UNPACK #-} !SCALAR_TYPE t t+ | Unary {-# UNPACK #-} !SCALAR_TYPE {-# UNPACK #-} !SCALAR_TYPE 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 AD_TYPE = Kahn (Tape AD_TYPE) deriving (Show, Typeable)++instance MuRef AD_TYPE where+ type DeRef AD_TYPE = Tape++ 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 Mode AD_TYPE where+ type Scalar AD_TYPE = SCALAR_TYPE++ isKnownZero (Kahn Zero) = True+ isKnownZero (Kahn (Lift 0)) = True+ isKnownZero _ = False++ asKnownConstant (Kahn Zero) = Just 0+ asKnownConstant (Kahn (Lift n)) = Just n+ asKnownConstant _ = Nothing++ isKnownConstant (Kahn Zero) = True+ isKnownConstant (Kahn (Lift _)) = True+ isKnownConstant _ = False++ 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++(<+>) :: AD_TYPE -> AD_TYPE -> AD_TYPE+(<+>) = binary (+) 1 1++primal :: AD_TYPE -> SCALAR_TYPE+primal (Kahn Zero) = 0+primal (Kahn (Lift a)) = a+primal (Kahn (Var a _)) = a+primal (Kahn (Binary a _ _ _ _)) = a+primal (Kahn (Unary a _ _)) = a++instance Jacobian AD_TYPE where+ type D AD_TYPE = Id SCALAR_TYPE++ 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 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 _ _ (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 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)+++mul :: AD_TYPE -> AD_TYPE -> AD_TYPE+mul = lift2 (*) (\x y -> (y, x))++#define HEAD AD_TYPE+#define BODY1(x)+#define BODY2(x,y)+#define NO_Bounded+#include <instances.h>++derivative+ :: AD_TYPE -> SCALAR_TYPE+derivative = sum . map snd . partials+{-# INLINE derivative #-}++derivative'+ :: AD_TYPE -> (SCALAR_TYPE, SCALAR_TYPE)+derivative' r = (primal r, derivative r)+{-# INLINE derivative' #-}++-- | back propagate sensitivities along a tape.+backPropagate :: (Vertex -> (Tape Int, Int, [Int])) -> STUArray s Int SCALAR_TYPE -> 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)+ Binary _ gb gc b c -> do+ da <- readArray ss i+ db <- readArray ss b+ writeArray ss b (db + gb*da)+ dc <- readArray ss c+ writeArray ss c (dc + gc*da)+ _ -> return ()+ where+ (node, i, _) = vmap v+ -- this isn't _quite_ right, as it should allow negative zeros to multiply through++topSortAcyclic :: Graph -> [Vertex]+topSortAcyclic g = 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 []++-- | This returns a list of contributions to the partials.+-- The variable ids returned in the list are likely /not/ unique!+partials :: AD_TYPE -> [(Int, SCALAR_TYPE)]++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 = A.array xsBounds xs+ vmap i = (vertexMap ! i, i, [])+ xsBounds = sbounds xs++ sensitivities = runSTUArray $ do+ ss <- newArray xsBounds 0+ writeArray ss start 1+ forM_ (topSortAcyclic g) $+ backPropagate vmap ss+ return ss++ sbounds ((a,_):as) = List.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 Int -> [Int]+ successors (Unary _ _ b) = [b]+ successors (Binary _ _ _ b c) = if b == c then [b] else [b,c]+ successors _ = []++-- | Return an 'Array' of 'partials' given bounds for the variable IDs.+partialArray :: (Int, Int) -> AD_TYPE -> UArray Int SCALAR_TYPE+partialArray vbounds tape = accumArray (+) 0 vbounds (partials tape)+{-# INLINE partialArray #-}++-- | Return an 'IntMap' of sparse partials+partialMap :: AD_TYPE -> IntMap SCALAR_TYPE+partialMap = fromListWith (+) . partials+{-# INLINE partialMap #-}++-- strict list of scalars+data List = Nil | Cons !SCALAR_TYPE !List++instance Exts.IsList List where+ type Item List = SCALAR_TYPE+ fromList (x:xs) = Cons x (Exts.fromList xs)+ fromList [] = Nil+ toList Nil = []+ toList (Cons x xs) = x : Exts.toList xs++class Grad i o o' | i -> o o', o -> i o', o' -> i o where+ pack :: i -> [AD_TYPE] -> AD_TYPE+ unpack :: (List -> List) -> o+ unpack'+ :: (List -> (SCALAR_TYPE, List))+ -> o'++instance Grad AD_TYPE List (SCALAR_TYPE, List) where+ pack i _ = i+ unpack f = f Nil+ unpack' f = f Nil++instance Grad i o o'+ => Grad (AD_TYPE -> i) (SCALAR_TYPE -> o) (SCALAR_TYPE -> o') where+ pack f (a:as) = pack (f a) as+ pack _ [] = error "Grad.pack: logic error"+ unpack f a = unpack (f . Cons a)+ unpack' f a = unpack' (f . Cons a)++vgrad :: Grad i o o' => i -> o+vgrad i = unpack (unsafeGrad (pack i)) where+ unsafeGrad f as = unbinds vs (partialArray bds $ f vs) where+ (vs,bds) = binds as++vgrad' :: Grad i o o' => i -> o'+vgrad' i = unpack' (unsafeGrad' (pack i)) where+ unsafeGrad' f as = (primal r, unbinds vs (partialArray bds r)) where+ r = f vs+ (vs,bds) = binds as++var :: SCALAR_TYPE -> Int -> AD_TYPE+var a v = Kahn (Var a v)++varId :: AD_TYPE -> Int+varId (Kahn (Var _ v)) = v+varId _ = error "varId: not a Var"++bind :: Traversable f => f SCALAR_TYPE -> (f AD_TYPE, (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')++binds :: List -> ([AD_TYPE], (Int,Int))+binds = bind . Exts.toList++unbind :: Functor f => f AD_TYPE -> UArray Int SCALAR_TYPE -> f SCALAR_TYPE+unbind xs ys = fmap (\v -> ys ! varId v) xs++unbinds :: Foldable f => f AD_TYPE -> UArray Int SCALAR_TYPE -> List+unbinds xs ys = foldr (\v r -> Cons (ys ! varId v) r) Nil xs++unbindWithUArray :: (Functor f, IArray UArray b) => (SCALAR_TYPE -> b -> c) -> f AD_TYPE -> UArray Int b -> f c+unbindWithUArray f xs ys = fmap (\v -> f (primal v) (ys ! varId v)) xs++unbindWithArray :: Functor f => (SCALAR_TYPE -> b -> c) -> f AD_TYPE -> Array Int b -> f c+unbindWithArray f xs ys = fmap (\v -> f (primal v) (ys ! varId v)) xs++unbindMap :: Functor f => f AD_TYPE -> IntMap SCALAR_TYPE -> f SCALAR_TYPE+unbindMap xs ys = fmap (\v -> findWithDefault 0 (varId v) ys) xs++unbindMapWithDefault :: Functor f => b -> (SCALAR_TYPE -> b -> c) -> f AD_TYPE -> IntMap b -> f c+unbindMapWithDefault z f xs ys = fmap (\v -> f (primal v) $ findWithDefault z (varId v) ys) xs
+ include/rank1_kahn.h view
@@ -0,0 +1,270 @@+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE UndecidableInstances #-}+-----------------------------------------------------------------------------+-- |+-- Copyright : (c) Edward Kmett 2010-2021+-- License : BSD3+-- Maintainer : ekmett@gmail.com+-- Stability : experimental+-- Portability : GHC only+--+-- This module provides reverse-mode Automatic Differentiation using post-hoc linear time+-- topological sorting.+--+-- For reverse mode AD we use 'System.Mem.StableName.StableName' to recover sharing information from+-- the tape to avoid combinatorial explosion, and thus run asymptotically faster+-- than it could without such sharing information, but the use of side-effects+-- contained herein is benign.+--+-----------------------------------------------------------------------------++MODULE+ ( AD_EXPORT+ , auto+ -- * 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 Data.Functor.Compose+import Numeric.AD.Internal.On+import Numeric.AD.Mode+IMPORTS++-- $setup+--+-- >>> import Numeric.AD.Internal.Doctest++-- | 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.0,2.0,3.0]+-- [2.0,1.0,1.0]+grad+ :: BASE1_1(Traversable f, Num a)+ => (f AD_TYPE -> AD_TYPE)+ -> f SCALAR_TYPE+ -> f SCALAR_TYPE+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 kahn-mode AD in a single pass.+--+-- >>> grad' (\[x,y,z] -> x*y+z) [1.0,2.0,3.0]+-- (5.0,[2.0,1.0,1.0])+grad'+ :: BASE1_1(Traversable f, Num a)+ => (f AD_TYPE -> AD_TYPE)+ -> f SCALAR_TYPE+ -> (SCALAR_TYPE, f SCALAR_TYPE)+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 kahn-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+ :: BASE1_1(Traversable f, Num a)+ => (SCALAR_TYPE -> SCALAR_TYPE -> b)+ -> (f AD_TYPE -> AD_TYPE)+ -> f SCALAR_TYPE+ -> 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 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'+ :: BASE1_1(Traversable f, Num a)+ => (SCALAR_TYPE -> SCALAR_TYPE -> b)+ -> (f AD_TYPE -> AD_TYPE)+ -> f SCALAR_TYPE+ -> (SCALAR_TYPE, 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 kahn AD lazily in @m@ passes for @m@ outputs.+--+-- >>> jacobian (\[x,y] -> [y,x,x*y]) [2.0,1.0]+-- [[0.0,1.0],[1.0,0.0],[1.0,2.0]]+jacobian+ :: BASE2_1(Traversable f, Functor g, Num a)+ => (f AD_TYPE -> g AD_TYPE)+ -> f SCALAR_TYPE+ -> g (f SCALAR_TYPE)+jacobian f as = unbind vs . partialArray bds <$> f vs where+ (vs, bds) = bind as+{-# INLINE jacobian #-}++-- | The 'jacobian'' function calculates both the result and the Jacobian of a nonscalar-to-nonscalar function, using @m@ invocations of 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.0,1.0]+-- [(1.0,[0.0,1.0]),(2.0,[1.0,0.0]),(2.0,[1.0,2.0])]+jacobian'+ :: BASE2_1(Traversable f, Functor g, Num a)+ => (f AD_TYPE -> g AD_TYPE)+ -> f SCALAR_TYPE+ -> g (SCALAR_TYPE, f SCALAR_TYPE)+jacobian' f as = row <$> f vs where+ (vs, bds) = bind as+ row a = (primal a, unbind vs (partialArray bds a))+{-# INLINE jacobian' #-}++-- | 'jacobianWith g f' calculates the Jacobian of a non-scalar-to-non-scalar function @f@ with 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@.+--+-- @+-- 'jacobian' = 'jacobianWith' (\_ dx -> dx)+-- 'jacobianWith' 'const' = (\f x -> 'const' x '<$>' f x)+-- @+jacobianWith+ :: BASE2_1(Traversable f, Functor g, Num a)+ => (SCALAR_TYPE -> SCALAR_TYPE -> b)+ -> (f AD_TYPE -> g AD_TYPE)+ -> f SCALAR_TYPE+ -> g (f b)+jacobianWith g f as = UNBINDWITH g vs . partialArray bds <$> f vs where+ (vs, bds) = bind as+{-# INLINE jacobianWith #-}++-- | 'jacobianWith' g f' calculates both the result and the Jacobian of a nonscalar-to-nonscalar function @f@, using @m@ invocations of 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'+ :: BASE2_1(Traversable f, Functor g, Num a)+ => (SCALAR_TYPE -> SCALAR_TYPE -> b)+ -> (f AD_TYPE -> g AD_TYPE)+ -> f SCALAR_TYPE+ -> g (SCALAR_TYPE, f b)+jacobianWith' g f as = row <$> f vs where+ (vs, bds) = bind as+ row a = (primal a, UNBINDWITH g vs (partialArray bds a))+{-# INLINE jacobianWith' #-}++-- | Compute the derivative of a function.+--+-- >>> diff sin 0+-- 1.0+--+-- >>> cos 0+-- 1.0+diff :: BASE0_1(Num a)+ (AD_TYPE -> AD_TYPE)+ -> SCALAR_TYPE+ -> SCALAR_TYPE+diff f a = derivative $ f (var a 0)+{-# INLINE diff #-}++-- | The 'diff'' function calculates the value and derivative, as a+-- pair, of a scalar-to-scalar function.+--+--+-- >>> diff' sin 0+-- (0.0,1.0)+diff'+ :: BASE0_1(Num a)+ (AD_TYPE -> AD_TYPE)+ -> SCALAR_TYPE+ -> (SCALAR_TYPE, SCALAR_TYPE)+diff' f a = derivative' $ f (var a 0)+{-# INLINE diff' #-}++-- | Compute the derivatives of a function that returns a vector with regards to its single input.+--+-- >>> diffF (\a -> [sin a, cos a]) 0+-- [1.0,0.0]+diffF+ :: BASE1_1(Functor f, Num a)+ => (AD_TYPE -> f AD_TYPE)+ -> SCALAR_TYPE+ -> f SCALAR_TYPE+diffF f a = derivative <$> f (var a 0)+{-# INLINE diffF #-}++-- | Compute the derivatives of a function that returns a vector with regards to its single input+-- as well as the primal answer.+--+-- >>> diffF' (\a -> [sin a, cos a]) 0+-- [(0.0,1.0),(1.0,0.0)]+diffF'+ :: BASE1_1(Functor f, Num a)+ => (AD_TYPE -> f AD_TYPE)+ -> SCALAR_TYPE+ -> f (SCALAR_TYPE, SCALAR_TYPE)+diffF' f a = derivative' <$> f (var a 0)+{-# INLINE diffF' #-}++-- | 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.0,2.0]+-- [[0.0,1.0],[1.0,0.0]]+hessian+ :: BASE1_1(Traversable f, Num a)+ => (f (On (Kahn AD_TYPE))+ -> On (Kahn AD_TYPE))+ -> f SCALAR_TYPE+ -> f (f SCALAR_TYPE)+hessian f = jacobian (GRAD (off . f . fmap On))++-- | 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+ :: BASE2_1(Traversable f, Functor g, Num a)+ => (f (On (Kahn AD_TYPE)) -> g (On (Kahn AD_TYPE)))+ -> f SCALAR_TYPE+ -> g (f (f SCALAR_TYPE))+hessianF f = getCompose . jacobian (Compose . JACOBIAN (fmap off . f . fmap On))++-- $vgrad+--+-- Variadic combinators for variadic mixed-mode automatic differentiation.+--+-- Unfortunately, variadicity comes at the expense of being able to use+-- quantification to avoid sensitivity confusion, so be careful when+-- counting the number of 'auto' calls you use when taking the gradient+-- of a function that takes gradients!
src/Numeric/AD.hs view
@@ -1,11 +1,10 @@-{-# LANGUAGE CPP #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE Rank2Types #-} {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE PatternGuards #-} ----------------------------------------------------------------------------- -- |--- Copyright : (c) Edward Kmett 2010-2015+-- Copyright : (c) Edward Kmett 2010-2021 -- License : BSD3 -- Maintainer : ekmett@gmail.com -- Stability : experimental@@ -129,13 +128,17 @@ , conjugateGradientDescent , conjugateGradientAscent , stochasticGradientDescent++ -- * Working with towers+ , Jet(..)+ , headJet+ , tailJet+ , jet ) where import Data.Functor.Compose-#if __GLASGOW_HASKELL__ < 710-import Data.Traversable (Traversable)-#endif import Data.Reflection (Reifies)+import Data.Typeable import Numeric.AD.Internal.Forward (Forward) import Numeric.AD.Internal.Kahn (Grad, vgrad, vgrad') import Numeric.AD.Internal.On@@ -143,6 +146,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,35 +180,63 @@ -- -- >>> 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 -- -- 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. Reifies s Tape => f (On (Reverse s (Forward a))) -> On (Reverse s (Forward a))) -> f (a, a) -> f a+hessianProduct+ :: (Traversable f, Num a)+ => ( forall s. (Reifies s Tape, Typeable s)+ => 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 -- -- 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. Reifies s Tape => f (On (Reverse s (Forward a))) -> On (Reverse s (Forward a))) -> f (a, a) -> f (a, a)+hessianProduct'+ :: (Traversable f, Num a)+ => ( forall s. (Reifies s Tape, Typeable s)+ => f (On (Reverse s (Forward a)))+ -> On (Reverse s (Forward a))+ )+ -> f (a, a)+ -> f (a, a) hessianProduct' f = Forward1.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 (\[x,y] -> x*y) [1,2] -- [[0,1],[1,0]]-hessian :: (Traversable f, Num a) => (forall s. Reifies s Tape => f (On (Reverse s (Sparse a))) -> On (Reverse s (Sparse a))) -> f a -> f (f a)+hessian+ :: (Traversable f, Num a)+ => ( forall s. (Reifies s Tape, Typeable s)+ => f (On (Reverse s (Sparse a)))+ -> On (Reverse s (Sparse a))+ )+ -> f a+ -> f (f a) hessian f = Sparse1.jacobian (grad (off . f . fmap On)) -- | 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 :: RDouble] -- [[[0.0,1.0],[1.0,0.0]],[[0.0,0.0],[0.0,0.0]],[[-1.131204383757,-2.471726672005],[-2.471726672005,1.131204383757]]]-hessianF :: (Traversable f, Functor g, Num a) => (forall s. Reifies s Tape => f (On (Reverse s (Sparse a))) -> g (On (Reverse s (Sparse a)))) -> f a -> g (f (f a))+hessianF+ :: (Traversable f, Functor g, Num a)+ => ( forall s. (Reifies s Tape, Typeable s)+ => f (On (Reverse s (Sparse a)))+ -> g (On (Reverse s (Sparse a)))+ )+ -> f a+ -> g (f (f a)) hessianF f as = getCompose $ Sparse1.jacobian (Compose . Reverse.jacobian (fmap off . f . fmap On)) as -- $vgrad
+ src/Numeric/AD/Double.hs view
@@ -0,0 +1,246 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE PatternGuards #-}+-----------------------------------------------------------------------------+-- |+-- Copyright : (c) Edward Kmett 2010-2021+-- License : BSD3+-- Maintainer : ekmett@gmail.com+-- Stability : experimental+-- Portability : GHC only+--+-- Mixed-Mode Automatic Differentiation, specialized to doubles.+--+-- Each combinator exported from this module chooses an appropriate AD mode.+-- The following basic operations are supported, modified as appropriate by the suffixes below:+--+-- * 'grad' computes the gradient (partial derivatives) of a function at a point+--+-- * 'jacobian' computes the Jacobian matrix of a function at a point+--+-- * 'diff' computes the derivative of a function at a point+--+-- * 'du' computes a directional derivative of a function at a point+--+-- * 'hessian' compute the Hessian matrix (matrix of second partial derivatives) of a function at a point+--+-- The suffixes have the following meanings:+--+-- * @\'@ -- also return the answer+--+-- * @With@ lets the user supply a function to blend the input with the output+--+-- * @F@ is a version of the base function lifted to return a 'Traversable' (or 'Functor') result+--+-- * @s@ means the function returns all higher derivatives in a list or f-branching 'Stream'+--+-- * @T@ means the result is transposed with respect to the traditional formulation.+--+-- * @0@ means that the resulting derivative list is padded with 0s at the end.+-----------------------------------------------------------------------------++module Numeric.AD.Double+ ( AD++ -- * AD modes+ , Mode(auto)+ , Scalar++ -- * Gradients (Reverse Mode)+ , grad+ , grad'+ , gradWith+ , gradWith'++ -- * Higher Order Gradients (Sparse-on-Reverse)+ , grads++ -- * Variadic Gradients (Sparse or Kahn)+ -- $vgrad+ , Grad , vgrad, vgrad'+ , Grads, vgrads++ -- * Jacobians (Sparse or Reverse)+ , jacobian+ , jacobian'+ , jacobianWith+ , jacobianWith'++ -- * Higher Order Jacobian (Sparse-on-Reverse)+ , jacobians++ -- * Transposed Jacobians (Forward Mode)+ , jacobianT+ , jacobianWithT++ -- * Hessian (Sparse-On-Reverse)+ , hessian+ , hessian'++ -- * Hessian Tensors (Sparse or Sparse-On-Reverse)+ , hessianF++ -- * Hessian Tensors (Sparse)+ , hessianF'++ -- * Hessian Vector Products (Forward-On-Reverse)+ , hessianProduct+ , hessianProduct'++ -- * Derivatives (Forward Mode)+ , diff+ , diffF++ , diff'+ , diffF'++ -- * Derivatives (Tower)+ , diffs+ , diffsF++ , diffs0+ , diffs0F++ -- * Directional Derivatives (Forward Mode)+ , du+ , du'+ , duF+ , duF'++ -- * Directional Derivatives (Tower)+ , dus+ , dus0+ , dusF+ , dus0F++ -- * Taylor Series (Tower)+ , taylor+ , taylor0++ -- * Maclaurin Series (Tower)+ , maclaurin+ , maclaurin0++ -- * Gradient Descent+ -- , gradientDescent+ -- , gradientAscent+ , conjugateGradientDescent+ , conjugateGradientAscent+ -- , stochasticGradientDescent++ -- * Working with towers+ , Jet(..)+ , headJet+ , tailJet+ , jet+ ) where++import Data.Functor.Compose+import Data.Reflection (Reifies)+import Data.Typeable+import Numeric.AD.Internal.Forward.Double (ForwardDouble)+import Numeric.AD.Internal.Kahn.Double (Grad, vgrad, vgrad')+import Numeric.AD.Internal.On+import Numeric.AD.Internal.Reverse (Reverse, Tape)+import Numeric.AD.Internal.Sparse.Double (SparseDouble, Grads, vgrads)++import Numeric.AD.Internal.Type+import Numeric.AD.Jet+import Numeric.AD.Mode++import qualified Numeric.AD.Rank1.Forward.Double as ForwardDouble1+import Numeric.AD.Mode.Forward.Double+ ( diff, diff', diffF, diffF'+ , du, du', duF, duF'+ , jacobianT, jacobianWithT+ )++import Numeric.AD.Mode.Tower.Double+ ( 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.Double+ ( grad, grad', gradWith, gradWith'+ , jacobian, jacobian', jacobianWith, jacobianWith'+ )++-- temporary until we make a full sparse mode+import qualified Numeric.AD.Rank1.Sparse.Double as SparseDouble1+import Numeric.AD.Mode.Sparse.Double+ ( grads, jacobians, hessian', hessianF'+ )++import Numeric.AD.Newton.Double++-- $setup+--+-- >>> 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:+--+-- > H v = (d/dr) grad_w (w + r v) | r = 0+--+-- Or in other words, we take the directional derivative of the gradient. The gradient is calculated in reverse mode, then the directional derivative is calculated in forward mode.+--+hessianProduct+ :: Traversable f+ => ( forall s. (Reifies s Tape, Typeable s)+ => f (On (Reverse s ForwardDouble))+ -> On (Reverse s ForwardDouble)+ )+ -> f (Double, Double)+ -> f Double+hessianProduct f = ForwardDouble1.duF (Reverse.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+ => ( forall s. (Reifies s Tape, Typeable s)+ => f (On (Reverse s ForwardDouble))+ -> On (Reverse s ForwardDouble)+ )+ -> f (Double, Double)+ -> f (Double, Double)+hessianProduct' f = ForwardDouble1.duF' (Reverse.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 (\[x,y] -> x*y) [1,2]+-- [[0.0,1.0],[1.0,0.0]]+hessian+ :: Traversable f+ => ( forall s. (Reifies s Tape, Typeable s)+ => f (On (Reverse s SparseDouble))+ -> On (Reverse s SparseDouble)+ )+ -> f Double+ -> f (f Double)+hessian f = SparseDouble1.jacobian (Reverse.grad (off . f . fmap On))++-- | Compute the order 3 Hessian tensor on a non-scalar-to-non-scalar function using 'Sparse'-on-'Reverse'+hessianF+ :: (Traversable f, Functor g)+ => ( forall s. (Reifies s Tape, Typeable s)+ => f (On (Reverse s SparseDouble))+ -> g (On (Reverse s SparseDouble))+ )+ -> f Double+ -> g (f (f Double))+hessianF f as = getCompose $ SparseDouble1.jacobian (Compose . Reverse.jacobian (fmap off . f . fmap On)) as++-- $vgrad+--+-- Variadic combinators for variadic mixed-mode automatic differentiation.+--+-- Unfortunately, variadicity comes at the expense of being able to use+-- quantification to avoid sensitivity confusion, so be careful when+-- counting the number of 'auto' calls you use when taking the gradient+-- of a function that takes gradients!
− src/Numeric/AD/Halley.hs
@@ -1,118 +0,0 @@-{-# LANGUAGE Rank2Types #-}-{-# LANGUAGE ScopedTypeVariables #-}--------------------------------------------------------------------------------- |--- Copyright : (c) Edward Kmett 2010-2015--- License : BSD3--- Maintainer : ekmett@gmail.com--- Stability : experimental--- Portability : GHC only------ Root finding using Halley's rational method (the second in--- the class of Householder methods). Assumes the function is three--- times continuously differentiable and converges cubically when--- progress can be made.-----------------------------------------------------------------------------------module Numeric.AD.Halley- (- -- * Halley's Method (Tower AD)- findZero- , findZeroNoEq- , inverse- , inverseNoEq- , fixedPoint- , fixedPointNoEq- , extremum- , extremumNoEq- ) where--import Prelude-import Numeric.AD.Internal.Forward (Forward)-import Numeric.AD.Internal.On-import Numeric.AD.Internal.Tower (Tower)-import Numeric.AD.Internal.Type (AD(..))-import qualified Numeric.AD.Rank1.Halley as Rank1---- $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--- ("it converges"), no further elements are returned.------ Examples:------ >>> take 10 $ findZero (\x->x^2-4) 1--- [1.0,1.8571428571428572,1.9997967892704736,1.9999999999994755,2.0]------ >>> last $ take 10 $ findZero ((+1).(^2)) (1 :+ 1)--- 0.0 :+ 1.0-findZero :: (Fractional a, Eq a) => (forall s. AD s (Tower a) -> AD s (Tower a)) -> a -> [a]-findZero f = Rank1.findZero (runAD.f.AD)-{-# INLINE findZero #-}---- | The 'findZeroNoEq' function behaves the same as 'findZero' except that it--- doesn't truncate the list once the results become constant. This means it--- can be used with types without an 'Eq' instance.-findZeroNoEq :: Fractional a => (forall s. AD s (Tower a) -> AD s (Tower a)) -> a -> [a]-findZeroNoEq f = Rank1.findZeroNoEq (runAD.f.AD)-{-# INLINE findZeroNoEq #-}---- | The 'inverse' function inverts 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--- ("it converges"), no further elements are returned.------ 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. AD s (Tower a) -> AD s (Tower a)) -> a -> a -> [a]-inverse f = Rank1.inverse (runAD.f.AD)-{-# INLINE inverse #-}---- | The 'inverseNoEq' function behaves the same as 'inverse' except that it--- doesn't truncate the list once the results become constant. This means it--- can be used with types without an 'Eq' instance.-inverseNoEq :: Fractional a => (forall s. AD s (Tower a) -> AD s (Tower a)) -> a -> a -> [a]-inverseNoEq f = Rank1.inverseNoEq (runAD.f.AD)-{-# INLINE inverseNoEq #-}---- | The 'fixedPoint' function find a fixedpoint of a scalar--- function using Halley's method; its output is a stream of--- increasingly accurate results. (Modulo the usual caveats.)------ If the stream becomes constant ("it converges"), no further--- elements are returned.------ >>> last $ take 10 $ fixedPoint cos 1--- 0.7390851332151607-fixedPoint :: (Fractional a, Eq a) => (forall s. AD s (Tower a) -> AD s (Tower a)) -> a -> [a]-fixedPoint f = Rank1.fixedPoint (runAD.f.AD)-{-# INLINE fixedPoint #-}---- | The 'fixedPointNoEq' function behaves the same as 'fixedPoint' except that--- it doesn't truncate the list once the results become constant. This means it--- can be used with types without an 'Eq' instance.-fixedPointNoEq :: Fractional a => (forall s. AD s (Tower a) -> AD s (Tower a)) -> a -> [a]-fixedPointNoEq f = Rank1.fixedPointNoEq (runAD.f.AD)-{-# INLINE fixedPointNoEq #-}---- | 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--- constant ("it converges"), no further elements are returned.------ >>> take 10 $ extremum cos 1--- [1.0,0.29616942658570555,4.59979519460002e-3,1.6220740159042513e-8,0.0]-extremum :: (Fractional a, Eq a) => (forall s. AD s (On (Forward (Tower a))) -> AD s (On (Forward (Tower a)))) -> a -> [a]-extremum f = Rank1.extremum (runAD.f.AD)-{-# INLINE extremum #-}---- | The 'extremumNoEq' function behaves the same as 'extremum' except that it--- doesn't truncate the list once the results become constant. This means it--- can be used with types without an 'Eq' instance.-extremumNoEq :: Fractional a => (forall s. AD s (On (Forward (Tower a))) -> AD s (On (Forward (Tower a)))) -> a -> [a]-extremumNoEq f = Rank1.extremumNoEq (runAD.f.AD)-{-# INLINE extremumNoEq #-}
+ src/Numeric/AD/Halley/Double.hs view
@@ -0,0 +1,83 @@+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE ScopedTypeVariables #-}+-----------------------------------------------------------------------------+-- |+-- Copyright : (c) Edward Kmett 2010-2021+-- License : BSD3+-- Maintainer : ekmett@gmail.com+-- Stability : experimental+-- Portability : GHC only+--+-- Root finding using Halley's rational method (the second in+-- the class of Householder methods). Assumes the function is three+-- times continuously differentiable and converges cubically when+-- progress can be made.+--+-----------------------------------------------------------------------------++module Numeric.AD.Halley.Double+ (+ -- * Halley's Method (Tower AD)+ findZero+ , inverse+ , fixedPoint+ , extremum+ ) where++import Prelude+import Numeric.AD.Internal.Forward (Forward)+import Numeric.AD.Internal.On+import Numeric.AD.Internal.Tower.Double (TowerDouble)+import Numeric.AD.Internal.Type (AD(..))+import qualified Numeric.AD.Rank1.Halley.Double as Rank1++-- $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+-- ("it converges"), no further elements are returned.+--+-- Examples:+--+-- >>> take 10 $ findZero (\x->x^2-4) 1+-- [1.0,1.8571428571428572,1.9997967892704736,1.9999999999994755,2.0]+findZero :: (forall s. AD s TowerDouble -> AD s TowerDouble) -> Double -> [Double]+findZero f = Rank1.findZero (runAD.f.AD)+{-# INLINE findZero #-}++-- | The 'inverse' function inverts a scalar function using+-- Halley's method; its output is a stream of increasingly accurate+-- results. (Modulo the usual caveats.) If the stream becomes constant+-- ("it converges"), no further elements are returned.+--+-- 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 :: (forall s. AD s TowerDouble -> AD s TowerDouble) -> Double -> Double -> [Double]+inverse f = Rank1.inverse (runAD.f.AD)+{-# INLINE inverse #-}++-- | The 'fixedPoint' function find a fixedpoint of a scalar+-- function using Halley's method; its output is a stream of+-- increasingly accurate results. (Modulo the usual caveats.)+--+-- If the stream becomes constant ("it converges"), no further+-- elements are returned.+--+-- >>> last $ take 10 $ fixedPoint cos 1+-- 0.7390851332151607+fixedPoint :: (forall s. AD s TowerDouble -> AD s TowerDouble) -> Double -> [Double]+fixedPoint f = Rank1.fixedPoint (runAD.f.AD)+{-# 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+-- constant ("it converges"), no further elements are returned.+--+-- >>> take 10 $ extremum cos 1+-- [1.0,0.29616942658570555,4.59979519460002e-3,1.6220740159042513e-8,0.0]+extremum :: (forall s. AD s (On (Forward TowerDouble)) -> AD s (On (Forward TowerDouble))) -> Double -> [Double]+extremum f = Rank1.extremum (runAD.f.AD)+{-# INLINE extremum #-}
src/Numeric/AD/Internal/Combinators.hs view
@@ -1,9 +1,8 @@-{-# LANGUAGE CPP #-} {-# OPTIONS_GHC -fno-warn-incomplete-patterns #-} {-# OPTIONS_HADDOCK not-home #-} ----------------------------------------------------------------------------- -- |--- Copyright : (c) Edward Kmett 2010-2015+-- Copyright : (c) Edward Kmett 2010-2021 -- License : BSD3 -- Maintainer : ekmett@gmail.com -- Stability : experimental@@ -19,19 +18,17 @@ , takeWhileDifferent ) where -#if __GLASGOW_HASKELL__ < 710-import Data.Traversable (Traversable, mapAccumL)-import Data.Foldable (Foldable, toList)-#else import Data.Traversable (mapAccumL) import Data.Foldable (toList)-#endif 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-zipWithT f as = snd . mapAccumL (\(a:as') b -> (as', f a b)) (toList as)+zipWithT f as = snd . mapAccumL f' (toList as)+ where+ f' (a:as') b = (as', f a b)+ f' [] _ = error "zipWithT: second argument contains less entries than third argument" -- | 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
src/Numeric/AD/Internal/Dense.hs view
@@ -11,7 +11,7 @@ ----------------------------------------------------------------------------- -- |--- Copyright : (c) Edward Kmett 2010-2015+-- Copyright : (c) Edward Kmett 2010-2021 -- License : BSD3 -- Maintainer : ekmett@gmail.com -- Stability : experimental@@ -41,18 +41,11 @@ ) where import Control.Monad (join)-#if __GLASGOW_HASKELL__ < 710-import Data.Functor-#endif import Data.Typeable ()-import Data.Traversable- ( mapAccumL-#if __GLASGOW_HASKELL__ < 710- , Traversable-#endif- )+import Data.Traversable (mapAccumL) import Data.Data () import Data.Number.Erf+import Numeric import Numeric.AD.Internal.Combinators import Numeric.AD.Internal.Identity import Numeric.AD.Jacobian@@ -97,6 +90,13 @@ instance (Num a, Traversable f) => Mode (Dense f a) where type Scalar (Dense f a) = a+ asKnownConstant (Lift a) = Just a+ asKnownConstant Zero = Just 0+ asKnownConstant _ = Nothing+ isKnownConstant Dense{} = False+ isKnownConstant _ = True+ isKnownZero Zero = True+ isKnownZero _ = False auto = Lift zero = Zero _ *^ Zero = Zero@@ -117,12 +117,6 @@ Dense a da <+> Lift b = Dense (a + b) da Dense a da <+> Dense b db = Dense (a + b) $ zipWithT (+) da db -(<**>) :: (Traversable f, Floating a) => Dense f a -> Dense f a -> Dense f 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 * log xi)) x y- instance (Traversable f, Num a) => Jacobian (Dense f a) where type D (Dense f a) = Id a unary f _ Zero = Lift (f 0)@@ -188,7 +182,7 @@ mul :: (Traversable f, Num a) => Dense f a -> Dense f a -> Dense f a mul = lift2 (*) (\x y -> (y, x)) -#define BODY1(x) (Traversable f, x)-#define BODY2(x,y) (Traversable f, x, y)-#define HEAD Dense f a+#define BODY1(x) (Traversable f, x) =>+#define BODY2(x,y) (Traversable f, x, y) =>+#define HEAD (Dense f a) #include "instances.h"
+ src/Numeric/AD/Internal/Dense/Representable.hs view
@@ -0,0 +1,178 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE UndecidableInstances #-}+{-# OPTIONS_HADDOCK not-home #-}++-----------------------------------------------------------------------------+-- |+-- Copyright : (c) Edward Kmett 2010-2021+-- License : BSD3+-- Maintainer : ekmett@gmail.com+-- Stability : experimental+-- Portability : GHC only+--+-- A dense forward AD based on representable functors. This allows for much larger+-- forward mode data types than 'Numeric.AD.Internal.Dense, as we only need+-- the ability to compare the representation of a functor for equality, rather+-- than put the representation on in a straight line like you have to with+-- 'Traversable'.+-----------------------------------------------------------------------------++module Numeric.AD.Internal.Dense.Representable+ ( Repr(..)+ , ds+ , ds'+ , vars+ , apply+ ) where++import Control.Monad (join)+import Data.Functor.Rep+import Data.Typeable ()+import Data.Data ()+import Data.Number.Erf+import Numeric+import Numeric.AD.Internal.Combinators+import Numeric.AD.Internal.Identity+import Numeric.AD.Jacobian+import Numeric.AD.Mode++data Repr f a+ = Lift !a+ | Repr !a (f a)+ | Zero++instance Show a => Show (Repr f a) where+ showsPrec d (Lift a) = showsPrec d a+ showsPrec d (Repr a _) = showsPrec d a+ showsPrec _ Zero = showString "0"++ds :: f a -> Repr f a -> f a+ds _ (Repr _ da) = da+ds z _ = z+{-# INLINE ds #-}++ds' :: Num a => f a -> Repr f a -> (a, f a)+ds' _ (Repr a da) = (a, da)+ds' z (Lift a) = (a, z)+ds' z Zero = (0, z)+{-# INLINE ds' #-}++-- Bind variables and count inputs+vars :: (Representable f, Eq (Rep f), Num a) => f a -> f (Repr f a)+vars = imapRep $ \i a -> Repr a $ tabulate $ \j -> if i == j then 1 else 0+{-# INLINE vars #-}++apply :: (Representable f, Eq (Rep f), Num a) => (f (Repr f a) -> b) -> f a -> b+apply f as = f (vars as)+{-# INLINE apply #-}++primal :: Num a => Repr f a -> a+primal Zero = 0+primal (Lift a) = a+primal (Repr a _) = a++instance (Representable f, Num a) => Mode (Repr f a) where+ type Scalar (Repr f a) = a+ asKnownConstant (Lift a) = Just a+ asKnownConstant Zero = Just 0+ asKnownConstant _ = Nothing+ isKnownConstant Repr{} = False+ isKnownConstant _ = True+ isKnownZero Zero = True+ isKnownZero _ = False+ auto = Lift+ zero = Zero+ _ *^ Zero = Zero+ a *^ Lift b = Lift (a * b)+ a *^ Repr b db = Repr (a * b) $ fmap (a*) db+ Zero ^* _ = Zero+ Lift a ^* b = Lift (a * b)+ Repr a da ^* b = Repr (a * b) $ fmap (*b) da+ Zero ^/ _ = Zero+ Lift a ^/ b = Lift (a / b)+ Repr a da ^/ b = Repr (a / b) $ fmap (/b) da++(<+>) :: (Representable f, Num a) => Repr f a -> Repr f a -> Repr f a+Zero <+> a = a+a <+> Zero = a+Lift a <+> Lift b = Lift (a + b)+Lift a <+> Repr b db = Repr (a + b) db+Repr a da <+> Lift b = Repr (a + b) da+Repr a da <+> Repr b db = Repr (a + b) $ liftR2 (+) da db++instance (Representable f, Num a) => Jacobian (Repr f a) where+ type D (Repr f a) = Id a+ unary f _ Zero = Lift (f 0)+ unary f _ (Lift b) = Lift (f b)+ unary f (Id dadb) (Repr b db) = Repr (f b) (fmap (dadb *) db)++ lift1 f _ Zero = Lift (f 0)+ lift1 f _ (Lift b) = Lift (f b)+ lift1 f df (Repr b db) = Repr (f b) (fmap (dadb *) db) where+ Id dadb = df (Id b)++ lift1_ f _ Zero = Lift (f 0)+ lift1_ f _ (Lift b) = Lift (f b)+ lift1_ f df (Repr b db) = Repr a (fmap (dadb *) db) where+ a = f b+ Id dadb = df (Id a) (Id b)++ binary f _ _ Zero Zero = Lift (f 0 0)+ binary f _ _ Zero (Lift c) = Lift (f 0 c)+ binary f _ _ (Lift b) Zero = Lift (f b 0)+ binary f _ _ (Lift b) (Lift c) = Lift (f b c)+ binary f _ (Id dadc) Zero (Repr c dc) = Repr (f 0 c) $ fmap (* dadc) dc+ binary f _ (Id dadc) (Lift b) (Repr c dc) = Repr (f b c) $ fmap (* dadc) dc+ binary f (Id dadb) _ (Repr b db) Zero = Repr (f b 0) $ fmap (dadb *) db+ binary f (Id dadb) _ (Repr b db) (Lift c) = Repr (f b c) $ fmap (dadb *) db+ binary f (Id dadb) (Id dadc) (Repr b db) (Repr c dc) = Repr (f b c) $ liftR2 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 (Repr c dc) = Repr (f 0 c) $ fmap (*dadc) dc where dadc = runId (snd (df (Id 0) (Id c)))+ lift2 f df (Lift b) (Repr c dc) = Repr (f b c) $ fmap (*dadc) dc where dadc = runId (snd (df (Id b) (Id c)))+ lift2 f df (Repr b db) Zero = Repr (f b 0) $ fmap (dadb*) db where dadb = runId (fst (df (Id b) (Id 0)))+ lift2 f df (Repr b db) (Lift c) = Repr (f b c) $ fmap (dadb*) db where dadb = runId (fst (df (Id b) (Id c)))+ lift2 f df (Repr b db) (Repr c dc) = Repr (f b c) da where+ (Id dadb, Id dadc) = df (Id b) (Id c)+ da = liftR2 productRule db dc+ productRule dbi dci = dadb * dbi + dci * dadc++ lift2_ f _ Zero Zero = Lift (f 0 0)+ lift2_ f _ Zero (Lift c) = Lift (f 0 c)+ lift2_ f _ (Lift b) Zero = Lift (f b 0)+ lift2_ f _ (Lift b) (Lift c) = Lift (f b c)+ lift2_ f df Zero (Repr c dc) = Repr a $ fmap (*dadc) dc where+ a = f 0 c+ (_, Id dadc) = df (Id a) (Id 0) (Id c)+ lift2_ f df (Lift b) (Repr c dc) = Repr a $ fmap (*dadc) dc where+ a = f b c+ (_, Id dadc) = df (Id a) (Id b) (Id c)+ lift2_ f df (Repr b db) Zero = Repr a $ fmap (dadb*) db where+ a = f b 0+ (Id dadb, _) = df (Id a) (Id b) (Id 0)+ lift2_ f df (Repr b db) (Lift c) = Repr a $ fmap (dadb*) db where+ a = f b c+ (Id dadb, _) = df (Id a) (Id b) (Id c)+ lift2_ f df (Repr b db) (Repr c dc) = Repr a $ liftR2 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++mul :: (Representable f, Num a) => Repr f a -> Repr f a -> Repr f a+mul = lift2 (*) (\x y -> (y, x))++#define BODY1(x) (Representable f, x) =>+#define BODY2(x,y) (Representable f, x, y) =>+#define HEAD (Repr f a)+#include "instances.h"
src/Numeric/AD/Internal/Doctest.hs view
@@ -2,7 +2,7 @@ ----------------------------------------------------------------------------- -- |--- Copyright : (c) Edward Kmett 2010-2019+-- Copyright : (c) Edward Kmett 2010-2021 -- License : BSD3 -- Maintainer : ekmett@gmail.com -- Stability : experimental
src/Numeric/AD/Internal/Forward.hs view
@@ -4,14 +4,13 @@ {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE Rank2Types #-}-{-# LANGUAGE TemplateHaskell #-} {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE UndecidableInstances #-} {-# OPTIONS_HADDOCK not-home #-} ----------------------------------------------------------------------------- -- |--- Copyright : (c) Edward Kmett 2010-2015+-- Copyright : (c) Edward Kmett 2010-2021 -- License : BSD3 -- Maintainer : ekmett@gmail.com -- Stability : experimental@@ -38,25 +37,16 @@ import Control.Monad (join)-#if __GLASGOW_HASKELL__ < 710-import Control.Applicative hiding ((<**>))-import Data.Foldable (Foldable, toList)-import Data.Traversable (Traversable, mapAccumL)-#else import Data.Foldable (toList) import Data.Traversable (mapAccumL)-#endif import Data.Data import Data.Number.Erf+import Numeric 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 = Forward !a a@@ -98,6 +88,10 @@ isKnownZero Zero = True isKnownZero _ = False + asKnownConstant Zero = Just 0+ asKnownConstant (Lift a) = Just a+ asKnownConstant _ = Nothing+ isKnownConstant Forward{} = False isKnownConstant _ = True @@ -121,12 +115,6 @@ Lift a <+> Forward b db = Forward (a + b) db Lift a <+> Lift b = Lift (a + b) -(<**>) :: Floating a => Forward a -> Forward a -> Forward 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 * log xi)) x y- instance Num a => Jacobian (Forward a) where type D (Forward a) = Id a @@ -181,7 +169,7 @@ (Id dadb, Id dadc) = df (Id a) (Id b) (Id c) da = dadb * db + dc * dadc -#define HEAD Forward a+#define HEAD (Forward a) #include "instances.h" bind :: (Traversable f, Num a) => (f (Forward a) -> b) -> f a -> f b
src/Numeric/AD/Internal/Forward/Double.hs view
@@ -4,14 +4,13 @@ {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE MultiParamTypeClasses #-}-{-# LANGUAGE TemplateHaskell #-} {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE UndecidableInstances #-} {-# OPTIONS_HADDOCK not-home #-} ----------------------------------------------------------------------------- ---- |----- Copyright : (c) Edward Kmett 2010-2015+---- Copyright : (c) Edward Kmett 2010-2021 ---- License : BSD3 ---- Maintainer : ekmett@gmail.com ---- Stability : experimental@@ -34,17 +33,11 @@ , transposeWith ) where -#if __GLASGOW_HASKELL__ < 710-import Control.Applicative hiding ((<**>))-import Data.Foldable (Foldable, toList)-import Data.Traversable (Traversable, mapAccumL)-#else import Data.Foldable (toList) import Data.Traversable (mapAccumL)-#endif import Control.Monad (join)-import Data.Function (on) import Data.Number.Erf+import Numeric import Numeric.AD.Internal.Combinators import Numeric.AD.Internal.Identity import Numeric.AD.Jacobian@@ -75,6 +68,9 @@ isKnownZero (ForwardDouble 0 0) = True isKnownZero _ = False + asKnownConstant (ForwardDouble x 0) = Just x+ asKnownConstant _ = Nothing+ isKnownConstant (ForwardDouble _ 0) = True isKnownConstant _ = False @@ -109,98 +105,11 @@ (Id dadb, Id dadc) = df (Id a) (Id b) (Id c) da = dadb * db + dc * dadc -instance Eq ForwardDouble where- (==) = on (==) primal--instance Ord ForwardDouble where- compare = on compare primal--instance Num ForwardDouble 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 where- fromRational 0 = zero- fromRational r = auto (fromRational r)- x / y = x * recip y- recip = lift1_ recip (const . negate . join (*))--instance Floating ForwardDouble 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 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 where- toRational = toRational . primal--instance RealFloat ForwardDouble 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 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 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 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)+#define HEAD ForwardDouble+#define BODY1(x)+#define BODY2(x,y)+#define NO_Bounded+#include "instances.h" bind :: Traversable f => (f ForwardDouble -> b) -> f Double -> f b bind f as = snd $ mapAccumL outer (0 :: Int) as where@@ -230,3 +139,6 @@ transposeWith f as = snd . mapAccumL go xss0 where go xss b = (tail <$> xss, f b (head <$> xss)) xss0 = toList <$> as++mul :: ForwardDouble -> ForwardDouble -> ForwardDouble+mul = lift2 (*) (\x y -> (y, x))
src/Numeric/AD/Internal/Identity.hs view
@@ -9,7 +9,7 @@ ----------------------------------------------------------------------------- -- |--- Copyright : (c) Edward Kmett 2010-2015+-- Copyright : (c) Edward Kmett 2010-2021 -- License : BSD3 -- Maintainer : ekmett@gmail.com -- Stability : experimental@@ -25,9 +25,6 @@ ) where import Data.Data (Data)-#if !(MIN_VERSION_base(4,8,0))-import Data.Monoid (Monoid(..))-#endif #if !(MIN_VERSION_base(4,11,0)) import Data.Semigroup (Semigroup(..)) #endif@@ -36,7 +33,11 @@ import Numeric.AD.Mode newtype Id a = Id { runId :: a } deriving- (Eq, Ord, Show, Enum, Bounded, Num, Real, Fractional, Floating, RealFrac, RealFloat, Semigroup, Monoid, Data, Typeable, Erf, InvErf)+ ( Eq, Ord, Show, Enum, Bounded+ , Num, Real, Fractional, Floating+ , RealFrac, RealFloat, Semigroup+ , Monoid, Data, Typeable, Erf, InvErf+ ) probe :: a -> Id a probe = Id
src/Numeric/AD/Internal/Kahn.hs view
@@ -12,7 +12,7 @@ ----------------------------------------------------------------------------- -- |--- Copyright : (c) Edward Kmett 2010-2015+-- Copyright : (c) Edward Kmett 2010-2021 -- License : BSD3 -- Maintainer : ekmett@gmail.com -- Stability : experimental@@ -48,15 +48,10 @@ , varId ) where -#if __GLASGOW_HASKELL__ < 710-import Prelude hiding (mapM)-import Control.Applicative (Applicative(..),(<$>))-import Data.Traversable (Traversable, mapM)-#endif import Control.Monad.ST import Control.Monad hiding (mapM) import Control.Monad.Trans.State-import Data.List (foldl')+import qualified Data.List as List (foldl') import Data.Array.ST import Data.Array import Data.IntMap (IntMap, fromListWith, findWithDefault)@@ -67,6 +62,7 @@ import System.IO.Unsafe (unsafePerformIO) import Data.Data (Data) import Data.Typeable (Typeable)+import Numeric import Numeric.AD.Internal.Combinators import Numeric.AD.Internal.Identity import Numeric.AD.Jacobian@@ -99,12 +95,17 @@ isKnownZero (Kahn Zero) = True isKnownZero _ = False + asKnownConstant (Kahn Zero) = Just 0+ asKnownConstant (Kahn (Lift n)) = Just n+ asKnownConstant _ = Nothing+ isKnownConstant (Kahn Zero) = True isKnownConstant (Kahn (Lift _)) = True isKnownConstant _ = False 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@@ -112,12 +113,6 @@ (<+>) :: Num a => Kahn a -> Kahn a -> Kahn a (<+>) = binary (+) 1 1 -(<**>) :: Floating a => Kahn a -> Kahn a -> Kahn 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 * xi ** (yi - 1), z * log xi)) x y- primal :: Num a => Kahn a -> a primal (Kahn Zero) = 0 primal (Kahn (Lift a)) = a@@ -162,7 +157,7 @@ mul :: Num a => Kahn a -> Kahn a -> Kahn a mul = lift2 (*) (\x y -> (y, x)) -#define HEAD Kahn a+#define HEAD (Kahn a) #include <instances.h> derivative :: Num a => Kahn a -> a@@ -228,12 +223,12 @@ 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 ((a,_):as) = List.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.
+ src/Numeric/AD/Internal/Kahn/Double.hs view
@@ -0,0 +1,9 @@+{-# LANGUAGE CPP #-}++#define MODULE \+module Numeric.AD.Internal.Kahn.Double+#define IMPORTS+#define AD_EXPORT KahnDouble(..)+#define AD_TYPE KahnDouble+#define SCALAR_TYPE Double+#include <internal_kahn.h>
+ src/Numeric/AD/Internal/Kahn/Float.hs view
@@ -0,0 +1,9 @@+{-# LANGUAGE CPP #-}++#define MODULE \+module Numeric.AD.Internal.Kahn.Float+#define IMPORTS+#define AD_EXPORT KahnFloat(..)+#define AD_TYPE KahnFloat+#define SCALAR_TYPE Float+#include <internal_kahn.h>
src/Numeric/AD/Internal/On.hs view
@@ -1,4 +1,3 @@-{-# LANGUAGE CPP #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE StandaloneDeriving #-} {-# LANGUAGE TypeFamilies #-}@@ -10,7 +9,7 @@ ----------------------------------------------------------------------------- -- |--- Copyright : (c) Edward Kmett 2010-2015+-- Copyright : (c) Edward Kmett 2010-2021 -- License : BSD3 -- Maintainer : ekmett@gmail.com -- Stability : experimental@@ -26,9 +25,6 @@ import Data.Data import Numeric.AD.Mode -#ifdef HLINT-#endif- ------------------------------------------------------------------------------ -- On ------------------------------------------------------------------------------@@ -41,8 +37,11 @@ , InvErf, RealFloat, Typeable ) -instance (Mode t, Mode (Scalar t)) => Mode (On t) where+instance (Mode t, Mode (Scalar t), Num (Scalar (Scalar t))) => Mode (On t) where type Scalar (On t) = Scalar (Scalar t) auto = On . auto . auto+ isKnownZero (On n) = isKnownZero n+ asKnownConstant (On n) = asKnownConstant n >>= asKnownConstant+ isKnownConstant (On n) = maybe False isKnownConstant (asKnownConstant n) a *^ On b = On (auto a *^ b) On a ^* b = On (a ^* auto b)
src/Numeric/AD/Internal/Or.hs view
@@ -1,16 +1,15 @@ {-# LANGUAGE CPP #-} {-# LANGUAGE GADTs #-} {-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE UndecidableInstances #-}-#if __GLASGOW_HASKELL__ >= 707 {-# LANGUAGE DeriveDataTypeable #-}-#endif {-# OPTIONS_HADDOCK not-home #-} ----------------------------------------------------------------------------- -- |--- Copyright : (c) Edward Kmett 2015+-- Copyright : (c) Edward Kmett 2021 -- License : BSD3 -- Maintainer : ekmett@gmail.com -- Stability : experimental@@ -28,13 +27,8 @@ , binary ) where -#if __GLASGOW_HASKELL__ < 710-import Control.Applicative-#endif import Data.Number.Erf-#if __GLASGOW_HASKELL__ >= 707 import Data.Typeable-#endif import Numeric.AD.Mode runL :: Or F a b -> a@@ -58,9 +52,6 @@ binary :: (a -> a -> a) -> (b -> b -> b) -> Or s a b -> Or s a b -> Or s a b binary f _ (L a) (L b) = L (f a b) binary _ g (R a) (R b) = R (g a b)-#if __GLASGOW_HASKELL__ < 800-binary _ _ _ _ = impossible-#endif data F data T@@ -79,29 +70,16 @@ data Or s a b where L :: a -> Or F a b R :: b -> Or T a b-#if __GLASGOW_HASKELL__ >= 707 deriving Typeable #endif-#endif -#if __GLASGOW_HASKELL__ < 800-impossible :: a-impossible = error "Numeric.AD.Internal.Or: impossible case"-#endif- instance (Eq a, Eq b) => Eq (Or s a b) where L a == L b = a == b R a == R b = a == b-#if __GLASGOW_HASKELL__ < 800- _ == _ = impossible-#endif instance (Ord a, Ord b) => Ord (Or s a b) where L a `compare` L b = compare a b R a `compare` R b = compare a b-#if __GLASGOW_HASKELL__ < 800- _ `compare` _ = impossible-#endif instance (Enum a, Enum b, Chosen s) => Enum (Or s a b) where pred = unary pred pred@@ -112,19 +90,10 @@ enumFrom (R a) = R <$> enumFrom a enumFromThen (L a) (L b) = L <$> enumFromThen a b enumFromThen (R a) (R b) = R <$> enumFromThen a b-#if __GLASGOW_HASKELL__ < 800- enumFromThen _ _ = impossible-#endif enumFromTo (L a) (L b) = L <$> enumFromTo a b enumFromTo (R a) (R b) = R <$> enumFromTo a b-#if __GLASGOW_HASKELL__ < 800- enumFromTo _ _ = impossible-#endif enumFromThenTo (L a) (L b) (L c) = L <$> enumFromThenTo a b c enumFromThenTo (R a) (R b) (R c) = R <$> enumFromThenTo a b c-#if __GLASGOW_HASKELL__ < 800- enumFromThenTo _ _ _ = impossible-#endif instance (Bounded a, Bounded b, Chosen s) => Bounded (Or s a b) where maxBound = choose maxBound maxBound@@ -209,6 +178,8 @@ type Scalar (Or s a b) = Scalar a auto = choose <$> auto <*> auto isKnownConstant = chosen isKnownConstant isKnownConstant+ asKnownConstant (L a) = asKnownConstant a+ asKnownConstant (R b) = asKnownConstant b isKnownZero = chosen isKnownZero isKnownZero x *^ L a = L (x *^ a) x *^ R a = R (x *^ a)
src/Numeric/AD/Internal/Reverse.hs view
@@ -6,7 +6,6 @@ {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE Rank2Types #-} {-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE TemplateHaskell #-} {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE UndecidableInstances #-} {-# OPTIONS_GHC -fno-full-laziness #-}@@ -14,7 +13,7 @@ ----------------------------------------------------------------------------- -- |--- Copyright : (c) Edward Kmett 2012-2015+-- Copyright : (c) Edward Kmett 2012-2021 -- License : BSD3 -- Maintainer : ekmett@gmail.com -- Stability : experimental@@ -37,6 +36,7 @@ , Head(..) , Cells(..) , reifyTape+ , reifyTypeableTape , partials , partialArrayOf , partialMapOf@@ -64,12 +64,9 @@ 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 import Numeric.AD.Internal.Combinators import Numeric.AD.Internal.Identity import Numeric.AD.Jacobian@@ -78,17 +75,11 @@ 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 -> a -> Cells -> Cells Binary :: {-# UNPACK #-} !Int -> {-# UNPACK #-} !Int -> a -> a -> Cells -> Cells-#endif dropCells :: Int -> Cells -> Cells dropCells 0 xs = xs@@ -128,13 +119,11 @@ 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 deriving (Show, Typeable)-#endif instance (Reifies s Tape, Num a) => Mode (Reverse s a) where type Scalar (Reverse s a) = a@@ -142,6 +131,10 @@ isKnownZero Zero = True isKnownZero _ = False + asKnownConstant Zero = Just 0+ asKnownConstant (Lift n) = Just n+ asKnownConstant _ = Nothing+ isKnownConstant Reverse{} = False isKnownConstant _ = True @@ -154,12 +147,6 @@ (<+>) :: (Reifies s Tape, Num a) => Reverse s a -> Reverse s a -> Reverse s a (<+>) = binary (+) 1 1 -(<**>) :: (Reifies s Tape, Floating a) => Reverse s a -> Reverse s a -> Reverse s 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 * xi ** (yi - 1), z * log xi)) x y- primal :: Num a => Reverse s a -> a primal Zero = 0 primal (Lift a) = a@@ -168,7 +155,7 @@ instance (Reifies s Tape, Num a) => Jacobian (Reverse s a) where type D (Reverse s a) = Id a - unary f _ (Zero) = Lift (f 0)+ 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 @@ -202,9 +189,9 @@ mul :: (Reifies s Tape, Num a) => Reverse s a -> Reverse s a -> Reverse s a mul = lift2 (*) (\x y -> (y, x)) -#define BODY1(x) (Reifies s Tape,x)-#define BODY2(x,y) (Reifies s Tape,x,y)-#define HEAD Reverse s a+#define BODY1(x) (Reifies s Tape,x) =>+#define BODY2(x,y) (Reifies s Tape,x,y) =>+#define HEAD (Reverse s a) #include "instances.h" -- | Helper that extracts the derivative of a chain when the chain was constructed with 1 variable.@@ -264,6 +251,13 @@ h <- newIORef (Head vs Nil) return (reify (Tape h) k) {-# NOINLINE reifyTape #-}++-- | Construct a tape that starts with @n@ variables.+reifyTypeableTape :: Int -> (forall s. (Typeable s, Reifies s Tape) => Proxy s -> r) -> r+reifyTypeableTape vs k = unsafePerformIO $ do+ h <- newIORef (Head vs Nil)+ return (reifyTypeable (Tape h) k)+{-# NOINLINE reifyTypeableTape #-} var :: a -> Int -> Reverse s a var a v = Reverse v a
+ src/Numeric/AD/Internal/Reverse/Double.hs view
@@ -0,0 +1,335 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE BangPatterns #-}+{-# 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 #-}++#ifdef AD_FFI+{-# LANGUAGE ForeignFunctionInterface #-}+#endif++module Numeric.AD.Internal.Reverse.Double+ ( ReverseDouble(..)+ , Tape(..)+ , reifyTape+ , reifyTypeableTape+ , partials+ , partialArrayOf+ , partialMapOf+ , derivativeOf+ , derivativeOf'+ , bind+ , unbind+ , unbindMap+ , unbindWith+ , unbindMapWithDefault+ , var+ , varId+ , primal+ ) where++#ifdef AD_FFI+import Foreign.Ptr+import Foreign.ForeignPtr+import Foreign.C.Types+import qualified Foreign.Marshal.Array as MA+import qualified Foreign.Marshal.Alloc as MA+#else+import Control.Monad.ST+import Data.Array.ST+import Data.Array.Unsafe as Unsafe+import Data.IORef+import Unsafe.Coerce+#endif++import Data.Functor+import Control.Monad hiding (mapM)+import Control.Monad.Trans.State+import Data.Array+import Data.IntMap (IntMap, fromDistinctAscList, findWithDefault)+import Data.Number.Erf+import Data.Proxy+import Data.Reflection+import Data.Traversable (mapM)+import Data.Typeable+import Numeric+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)++#ifdef AD_FFI++newtype Tape = Tape { getTape :: ForeignPtr Tape }++foreign import ccall unsafe "tape_alloc" c_tape_alloc :: CInt -> CInt -> IO (Ptr Tape)+foreign import ccall unsafe "tape_push" c_tape_push :: Ptr Tape -> CInt -> CInt -> Double -> Double -> IO Int+foreign import ccall unsafe "tape_backPropagate" c_tape_backPropagate :: Ptr Tape -> CInt -> Ptr Double -> IO ()+foreign import ccall unsafe "tape_variables" c_tape_variables :: Ptr Tape -> IO CInt+foreign import ccall unsafe "&tape_free" c_ref_tape_free :: FinalizerPtr Tape++pushTape :: Reifies s Tape => p s -> Int -> Int -> Double -> Double -> IO Int+pushTape p i1 i2 d1 d2 = do+ withForeignPtr (getTape (reflect p)) $ \tape -> + c_tape_push tape (fromIntegral i1) (fromIntegral i2) d1 d2+{-# INLINE pushTape #-}++-- | 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 _) = unsafePerformIO $+ withForeignPtr (getTape (reflect (Proxy :: Proxy s))) $ \tape -> do+ l <- fromIntegral <$> c_tape_variables tape+ arr <- MA.mallocArray l+ c_tape_backPropagate tape (fromIntegral k) arr+ ps <- MA.peekArray l arr+ MA.free arr+ return ps+{-# INLINE partials #-}++newTape :: Int -> IO Tape+newTape vs = do+ p <- c_tape_alloc (fromIntegral vs) (4 * 1024)+ Tape <$> newForeignPtr c_ref_tape_free p++-- | Construct a tape that starts with @n@ variables.+reifyTape :: Int -> (forall s. Reifies s Tape => Proxy s -> r) -> r+reifyTape vs k = unsafePerformIO $ fmap (\t -> reify t k) (newTape vs)+{-# NOINLINE reifyTape #-}++-- | Construct a tape that starts with @n@ variables.+reifyTypeableTape :: Int -> (forall s. (Reifies s Tape, Typeable s) => Proxy s -> r) -> r+reifyTypeableTape vs k = unsafePerformIO $ fmap (\t -> reifyTypeable t k) (newTape vs)+{-# NOINLINE reifyTypeableTape #-}++-- | 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 (pushTape (Proxy :: Proxy s) i (-1) di 0.0)) $! 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 (pushTape (Proxy :: Proxy s) i j di dj)) $! f b c+{-# INLINE binarily #-}++#else++data Cells where+ Nil :: Cells+ Unary :: {-# UNPACK #-} !Int -> {-# UNPACK #-} !Double -> !Cells -> Cells+ Binary :: {-# UNPACK #-} !Int -> {-# UNPACK #-} !Int -> {-# UNPACK #-} !Double -> {-# UNPACK #-} !Double -> !Cells -> Cells++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 }++-- | 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)++-- | 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 #-}++-- | Construct a tape that starts with @n@ variables.+reifyTypeableTape :: Int -> (forall s. (Reifies s Tape, Typeable s) => Proxy s -> r) -> r+reifyTypeableTape vs k = unsafePerformIO $ do+ h <- newIORef (Head vs Nil)+ return (reifyTypeable (Tape h) k)+{-# NOINLINE reifyTypeableTape #-}++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 #-}++#endif++data ReverseDouble s where+ Zero :: ReverseDouble s+ Lift :: {-# UNPACK #-} !Double -> ReverseDouble s+ ReverseDouble :: {-# UNPACK #-} !Int -> {-# UNPACK #-} !Double -> ReverseDouble s+ deriving (Show, Typeable)++instance Reifies s Tape => Mode (ReverseDouble s) where+ type Scalar (ReverseDouble s) = Double++ isKnownZero Zero = True+ isKnownZero (Lift 0) = 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++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)++mul :: Reifies s Tape => ReverseDouble s -> ReverseDouble s -> ReverseDouble s+mul = lift2 (*) (\x y -> (y, x))++#define BODY1(x) Reifies s Tape =>+#define BODY2(x,y) Reifies s Tape =>+#define HEAD (ReverseDouble s)+#define NO_Bounded+#include "instances.h"++-- | 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' #-}+++-- | 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 #-}++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
src/Numeric/AD/Internal/Sparse.hs view
@@ -12,7 +12,7 @@ {-# OPTIONS_HADDOCK not-home #-} ----------------------------------------------------------------------------- -- |--- Copyright : (c) Edward Kmett 2010-2015+-- Copyright : (c) Edward Kmett 2010-2021 -- License : BSD3 -- Maintainer : ekmett@gmail.com -- Stability : experimental@@ -44,35 +44,20 @@ ) where import Prelude hiding (lookup)-#if __GLASGOW_HASKELL__ < 710-import Control.Applicative hiding ((<**>))-#endif import Control.Comonad.Cofree-import Control.Monad (join)+import Control.Monad (join, guard) import Data.Data-import Data.IntMap (IntMap, unionWith, findWithDefault, toAscList, singleton, insertWith, lookup)+import Data.IntMap (IntMap, unionWith, findWithDefault, singleton, lookup) import qualified Data.IntMap as IntMap import Data.Number.Erf import Data.Traversable import Data.Typeable ()+import Numeric import Numeric.AD.Internal.Combinators+import Numeric.AD.Internal.Sparse.Common import Numeric.AD.Jacobian import Numeric.AD.Mode -newtype Monomial = Monomial (IntMap Int)--emptyMonomial :: Monomial-emptyMonomial = Monomial IntMap.empty-{-# INLINE emptyMonomial #-}--addToMonomial :: Int -> Monomial -> Monomial-addToMonomial k (Monomial m) = Monomial (insertWith (+) k 1 m)-{-# INLINE addToMonomial #-}--indices :: Monomial -> [Int]-indices (Monomial as) = uncurry (flip replicate) `concatMap` toAscList as-{-# INLINE indices #-}- -- | We only store partials in sorted order, so the map contained in a partial -- will only contain partials with equal or greater keys to that of the map in -- which it was found. This should be key for efficiently computing sparse hessians.@@ -92,12 +77,8 @@ apply f = f . vars {-# INLINE apply #-} -skeleton :: Traversable f => f a -> f Int-skeleton = snd . mapAccumL (\ !n _ -> (n + 1, n)) 0-{-# INLINE skeleton #-}- d :: (Traversable f, Num a) => f b -> Sparse a -> f a-d fs (Zero) = 0 <$ fs+d fs Zero = 0 <$ fs d fs (Sparse _ da) = snd $ mapAccumL (\ !n _ -> (n + 1, maybe 0 primal $ lookup n da)) 0 fs {-# INLINE d #-} @@ -108,7 +89,7 @@ ds :: (Traversable f, Num a) => f b -> Sparse a -> Cofree f a ds fs Zero = r where r = 0 :< (r <$ fs)-ds fs (as@(Sparse a _)) = a :< (go emptyMonomial <$> fns) where+ds fs as@(Sparse a _) = a :< (go emptyMonomial <$> fns) where fns = skeleton fs -- go :: Monomial -> Int -> Cofree f a go ix i = partial (indices ix') as :< (go ix' <$> fns) where@@ -139,17 +120,16 @@ primal (Sparse a _) = a primal Zero = 0 -(<**>) :: Floating a => Sparse a -> Sparse a -> Sparse a-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 * xi ** (yi - 1), z * log xi)) x y- instance Num a => Mode (Sparse a) where type Scalar (Sparse a) = a auto a = Sparse a IntMap.empty zero = Zero+ isKnownZero Zero = True+ isKnownZero _ = False+ isKnownConstant Zero = True+ isKnownConstant (Sparse _ m) = null m+ asKnownConstant Zero = Just 0+ asKnownConstant (Sparse a m) = a <$ guard (null m) Zero ^* _ = Zero Sparse a as ^* b = Sparse (a * b) $ fmap (^* b) as _ *^ Zero = Zero@@ -176,7 +156,7 @@ lift1_ f _ Zero = auto (f 0) lift1_ f df b@(Sparse pb bs) = a where- a = Sparse (f pb) $ IntMap.map ((df a b) *) bs+ a = Sparse (f pb) $ IntMap.map (df a b *) bs binary f _ _ Zero Zero = auto (f 0 0) binary f _ dadc Zero (Sparse pc dc) = Sparse (f 0 pc) $ IntMap.map (dadc *) dc@@ -199,8 +179,7 @@ a = Sparse (f pb pc) da da = unionWith (<+>) (IntMap.map (dadb *) db) (IntMap.map (dadc *) dc) --#define HEAD Sparse a+#define HEAD (Sparse a) #include "instances.h" class Num a => Grad i o o' a | i -> a o o', o -> a i o', o' -> a i o where@@ -250,29 +229,6 @@ isZero :: Sparse a -> Bool isZero Zero = True isZero _ = False---- |--- The value of the derivative of (f*g) of order mi is------ @--- 'sum' [a * 'primal' ('partialS' ('indices' b) f) * 'primal' ('partialS' ('indices' c) g) | (a,b,c) <- 'terms' mi ]--- @------ It is a bit more complicated in 'mul' below, since we build the whole tree of--- derivatives and want to prune the tree with 'Zero's as much as possible.--- The number of terms in the sum for order mi as of differentiation has--- @'sum' ('map' (+1) as)@ terms, so this is *much* more efficient--- than the naive recursive differentiation with @2^'sum' as@ terms.--- The coefficients @a@, which collect equivalent derivatives, are suitable products--- of binomial coefficients.-terms :: Monomial -> [(Integer,Monomial,Monomial)]-terms (Monomial m) = t (toAscList m) where- t [] = [(1,emptyMonomial,emptyMonomial)]- t ((k,a):ts) = concatMap (f (t ts)) (zip (bins!!a) [0..a]) where- f ps (b,i) = map (\(w,Monomial mf,Monomial mg) -> (w*b,Monomial (IntMap.insert k i mf), Monomial (IntMap.insert k (a-i) mg))) ps- bins = iterate next [1]- next xs@(_:ts) = 1 : zipWith (+) xs ts ++ [1]- next [] = error "impossible" mul :: Num a => Sparse a -> Sparse a -> Sparse a mul Zero _ = Zero
+ src/Numeric/AD/Internal/Sparse/Common.hs view
@@ -0,0 +1,54 @@+{-# LANGUAGE BangPatterns #-}+{-# OPTIONS_GHC -fno-warn-name-shadowing #-}+{-# OPTIONS_HADDOCK not-home #-}+-----------------------------------------------------------------------------+-- |+-- Copyright : (c) Edward Kmett 2010-2021+-- License : BSD3+-- Maintainer : ekmett@gmail.com+-- Stability : experimental+-- Portability : GHC only+--+-- common guts for Sparse.Double and Sparse mode+--+-- Handle with care.+-----------------------------------------------------------------------------+module Numeric.AD.Internal.Sparse.Common+ ( Monomial(..)+ , emptyMonomial+ , addToMonomial+ , indices+ , skeleton+ , terms+ ) where++import Data.IntMap (IntMap, toAscList, insertWith)+import qualified Data.IntMap as IntMap+import Data.Traversable++newtype Monomial = Monomial (IntMap Int)++emptyMonomial :: Monomial+emptyMonomial = Monomial IntMap.empty+{-# INLINE emptyMonomial #-}++addToMonomial :: Int -> Monomial -> Monomial+addToMonomial k (Monomial m) = Monomial (insertWith (+) k 1 m)+{-# INLINE addToMonomial #-}++indices :: Monomial -> [Int]+indices (Monomial as) = uncurry (flip replicate) `concatMap` toAscList as+{-# INLINE indices #-}++skeleton :: Traversable f => f a -> f Int+skeleton = snd . mapAccumL (\ !n _ -> (n + 1, n)) 0+{-# INLINE skeleton #-}++terms :: Monomial -> [(Integer,Monomial,Monomial)]+terms (Monomial m) = t (toAscList m) where+ t [] = [(1,emptyMonomial,emptyMonomial)]+ t ((k,a):ts) = concatMap (f (t ts)) (zip (bins!!a) [0..a]) where+ f ps (b,i) = map (\(w,Monomial mf,Monomial mg) -> (w*b,Monomial (IntMap.insert k i mf), Monomial (IntMap.insert k (a-i) mg))) ps+ bins = iterate next [1]+ next xs@(_:ts) = 1 : zipWith (+) xs ts ++ [1]+ next [] = error "impossible"
+ src/Numeric/AD/Internal/Sparse/Double.hs view
@@ -0,0 +1,258 @@+{-# 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 #-}+-----------------------------------------------------------------------------+-- |+-- Copyright : (c) Edward Kmett 2010-2021+-- 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.Sparse.Double+ ( Monomial(..)+ , emptyMonomial+ , addToMonomial+ , indices+ , SparseDouble(..)+ , apply+ , vars+ , d, d', ds+ , skeleton+ , spartial+ , partial+ , vgrad+ , vgrad'+ , vgrads+ , Grad(..)+ , Grads(..)+ , terms+ , primal+ ) where++import Prelude hiding (lookup)+import Control.Comonad.Cofree+import Control.Monad (join, guard)+import Data.Data+import Data.IntMap (IntMap, unionWith, findWithDefault, singleton, lookup)+import qualified Data.IntMap as IntMap+import Data.Number.Erf+import Data.Traversable+import Data.Typeable ()+import Numeric+import Numeric.AD.Internal.Combinators+import Numeric.AD.Internal.Sparse.Common+import Numeric.AD.Jacobian+import Numeric.AD.Mode++-- | We only store partials in sorted order, so the map contained in a partial+-- will only contain partials with equal or greater keys to that of the map in+-- which it was found. This should be key for efficiently computing sparse hessians.+-- there are only @n + k - 1@ choose @k@ distinct nth partial derivatives of a+-- function with k inputs.+data SparseDouble+ = Sparse {-# UNPACK #-} !Double (IntMap SparseDouble)+ | Zero+ deriving (Show, Data, Typeable)++vars :: Traversable f => f Double -> f SparseDouble+vars = snd . mapAccumL var 0 where+ var !n a = (n + 1, Sparse a $ singleton n $ auto 1)+{-# INLINE vars #-}++apply :: Traversable f => (f SparseDouble -> b) -> f Double -> b+apply f = f . vars+{-# INLINE apply #-}++d :: Traversable f => f b -> SparseDouble -> f Double+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 => f Double -> SparseDouble -> (Double, f Double)+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 => f b -> SparseDouble -> Cofree f Double+ds fs Zero = r where r = 0 :< (r <$ fs)+ds fs as@(Sparse a _) = a :< (go emptyMonomial <$> fns) where+ fns = skeleton fs+ -- go :: Monomial -> Int -> Cofree f a+ go ix i = partial (indices ix') as :< (go ix' <$> fns) where+ ix' = addToMonomial i ix+{-# INLINE ds #-}++partialS :: [Int] -> SparseDouble -> SparseDouble+partialS [] a = a+partialS (n:ns) (Sparse _ da) = partialS ns $ findWithDefault Zero n da+partialS _ Zero = Zero+{-# INLINE partialS #-}++partial :: [Int] -> SparseDouble -> Double+partial [] (Sparse a _) = a+partial (n:ns) (Sparse _ da) = partial ns $ findWithDefault (auto 0) n da+partial _ Zero = 0+{-# INLINE partial #-}++spartial :: [Int] -> SparseDouble -> Maybe Double+spartial [] (Sparse a _) = Just a+spartial (n:ns) (Sparse _ da) = do+ a' <- lookup n da+ spartial ns a'+spartial _ Zero = Nothing+{-# INLINE spartial #-}++primal :: SparseDouble -> Double+primal (Sparse a _) = a+primal Zero = 0++instance Mode SparseDouble where+ type Scalar SparseDouble = Double++ auto a = Sparse a IntMap.empty++ zero = Zero++ isKnownZero Zero = True+ isKnownZero (Sparse 0 m) = null m+ isKnownZero _ = False++ isKnownConstant Zero = True+ isKnownConstant (Sparse _ m) = null m++ asKnownConstant Zero = Just 0+ asKnownConstant (Sparse a m) = a <$ guard (null m)++ 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++infixr 6 <+>++(<+>) :: SparseDouble -> SparseDouble -> SparseDouble+Zero <+> a = a+a <+> Zero = a+Sparse a as <+> Sparse b bs = Sparse (a + b) $ unionWith (<+>) as bs++-- The instances for Jacobian for Sparse and Tower are almost identical;+-- could easily be made exactly equal by small changes.+instance Jacobian SparseDouble where+ type D SparseDouble = SparseDouble+ unary f _ Zero = auto (f 0)+ unary f dadb (Sparse pb bs) = Sparse (f pb) $ IntMap.map (* dadb) bs++ lift1 f _ Zero = auto (f 0)+ lift1 f df b@(Sparse pb bs) = Sparse (f pb) $ IntMap.map (* df b) bs++ lift1_ f _ Zero = auto (f 0)+ lift1_ f df b@(Sparse pb bs) = a where+ a = Sparse (f pb) $ IntMap.map (df a b *) bs++ binary f _ _ Zero Zero = auto (f 0 0)+ binary f _ dadc Zero (Sparse pc dc) = Sparse (f 0 pc) $ IntMap.map (dadc *) dc+ binary f dadb _ (Sparse pb db) Zero = Sparse (f pb 0 ) $ IntMap.map (dadb *) db+ binary f dadb dadc (Sparse pb db) (Sparse pc dc) = Sparse (f pb pc) $+ unionWith (<+>) (IntMap.map (dadb *) db) (IntMap.map (dadc *) dc)++ lift2 f _ Zero Zero = auto (f 0 0)+ lift2 f df Zero c@(Sparse pc dc) = Sparse (f 0 pc) $ IntMap.map (dadc *) dc where dadc = snd (df zero c)+ lift2 f df b@(Sparse pb db) Zero = Sparse (f pb 0) $ IntMap.map (* 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 (<+>) (IntMap.map (dadb *) db) (IntMap.map (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) (IntMap.map (fst (df a b zero) *) db)+ lift2_ f df Zero c@(Sparse pc dc) = a where a = Sparse (f 0 pc) (IntMap.map (* 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 (<+>) (IntMap.map (dadb *) db) (IntMap.map (dadc *) dc)++#define HEAD SparseDouble+#define BODY1(x)+#define BODY2(x,y)+#define NO_Bounded+#include "instances.h"++class Grad i o o' | i -> o o', o -> i o', o' -> i o where+ pack :: i -> [SparseDouble] -> SparseDouble+ unpack :: ([Double] -> [Double]) -> o+ unpack' :: ([Double] -> (Double, [Double])) -> o'++instance Grad SparseDouble [Double] (Double, [Double]) where+ pack i _ = i+ unpack f = f []+ unpack' f = f []++instance Grad i o o' => Grad (SparseDouble -> i) (Double -> o) (Double -> o') 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' => i -> o+vgrad i = unpack (unsafeGrad (pack i)) where+ unsafeGrad f as = d as $ apply f as+{-# INLINE vgrad #-}++vgrad' :: Grad i o o' => i -> o'+vgrad' i = unpack' (unsafeGrad' (pack i)) where+ unsafeGrad' f as = d' as $ apply f as+{-# INLINE vgrad' #-}++class Grads i o | i -> o, o -> i where+ packs :: i -> [SparseDouble] -> SparseDouble+ unpacks :: ([Double] -> Cofree [] Double) -> o++instance Grads SparseDouble (Cofree [] Double) where+ packs i _ = i+ unpacks f = f []++instance Grads i o => Grads (SparseDouble -> i) (Double -> o) 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 => i -> o+vgrads i = unpacks (unsafeGrads (packs i)) where+ unsafeGrads f as = ds as $ apply f as+{-# INLINE vgrads #-}++isZero :: SparseDouble -> Bool+isZero Zero = True+isZero _ = False++mul :: SparseDouble -> SparseDouble -> SparseDouble+mul Zero _ = Zero+mul _ Zero = Zero+mul f@(Sparse _ am) g@(Sparse _ bm) = Sparse (primal f * primal g) (derivs 0 emptyMonomial) where+ derivs v mi = IntMap.unions (map fn [v..kMax]) where+ fn w+ | and zs = IntMap.empty+ | otherwise = IntMap.singleton w (Sparse (sum ds) (derivs w mi'))+ where+ mi' = addToMonomial w mi+ (zs,ds) = unzip (map derVal (terms mi'))+ derVal (bin,mif,mig) = (isZero fder || isZero gder, fromIntegral bin * primal fder * primal gder) where+ fder = partialS (indices mif) f+ gder = partialS (indices mig) g+ kMax = maybe (-1) (fst.fst) (IntMap.maxViewWithKey am) `max` maybe (-1) (fst.fst) (IntMap.maxViewWithKey bm)
src/Numeric/AD/Internal/Tower.hs view
@@ -11,7 +11,7 @@ ----------------------------------------------------------------------------- -- |--- Copyright : (c) Edward Kmett 2010-2015+-- Copyright : (c) Edward Kmett 2010-2021 -- License : BSD3 -- Maintainer : ekmett@gmail.com -- Stability : experimental@@ -35,14 +35,12 @@ ) where import Prelude hiding (all, sum)-#if __GLASGOW_HASKELL__ < 710-import Control.Applicative hiding ((<**>))-#endif import Control.Monad (join) import Data.Foldable import Data.Data (Data) import Data.Number.Erf import Data.Typeable (Typeable)+import Numeric import Numeric.AD.Internal.Combinators import Numeric.AD.Jacobian import Numeric.AD.Mode@@ -60,7 +58,7 @@ {-# INLINE zeroPad #-} zeroPadF :: (Functor f, Num a) => [f a] -> [f a]-zeroPadF fxs@(fx:_) = fxs ++ repeat (const 0 <$> fx)+zeroPadF fxs@(fx:_) = fxs ++ repeat (0 <$ fx) zeroPadF _ = error "zeroPadF :: empty list" {-# INLINE zeroPadF #-} @@ -121,6 +119,13 @@ instance Num a => Mode (Tower a) where type Scalar (Tower a) = a auto a = Tower [a]+ isKnownZero (Tower xs) = null xs+ asKnownConstant (Tower []) = Just 0+ asKnownConstant (Tower [a]) = Just a+ asKnownConstant Tower {} = Nothing+ isKnownConstant (Tower []) = True+ isKnownConstant (Tower [_]) = True+ isKnownConstant Tower {} = False zero = Tower [] a *^ Tower bs = Tower (map (a*) bs) Tower as ^* b = Tower (map (*b) as)@@ -158,12 +163,6 @@ a = bundle a0 da (dadb, dadc) = df a b c -(<**>) :: Floating a => Tower a -> Tower a -> Tower a-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- -- mul xs ys = [ sum [xs!!j * ys!!(k-j)*bin k j | j <- [0..k]] | k <- [0..] ] -- adapted for efficiency and to handle finite lists xs, ys mul:: Num a => Tower a -> Tower a -> Tower a@@ -181,5 +180,5 @@ next' xs = zipWith (+) xs (tail xs) ++ [1] -- end part of next row in Pascal's triangle sumProd3 as bs cs = sum (zipWith3 (\x y z -> x*y*z) as bs cs) -#define HEAD Tower a+#define HEAD (Tower a) #include <instances.h>
+ src/Numeric/AD/Internal/Tower/Double.hs view
@@ -0,0 +1,271 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-# OPTIONS_GHC -fno-warn-name-shadowing #-}+{-# OPTIONS_HADDOCK not-home #-}++-----------------------------------------------------------------------------+-- |+-- Copyright : (c) Edward Kmett 2010-2021+-- License : BSD3+-- Maintainer : ekmett@gmail.com+-- Stability : experimental+-- Portability : GHC only+--+-----------------------------------------------------------------------------++#ifndef MIN_VERSION_base+#define MIN_VERSION_base(x,y,z) 1+#endif++module Numeric.AD.Internal.Tower.Double+ ( TowerDouble(..)+ , List(..)+ , zeroPad+ , zeroPadF+ , transposePadF+ , d, dl+ , d', dl'+ , withD+ , tangents+ , bundle+ , apply+ , getADTower+ , tower+ ) where++import Prelude hiding (all, sum)+import Control.Monad (join)+import Data.Foldable+import Data.Data (Data)+import Data.Number.Erf+import Data.Typeable (Typeable)+import Numeric+import Numeric.AD.Internal.Combinators+import Numeric.AD.Jacobian+import Numeric.AD.Mode+import Text.Read+import GHC.Exts as Exts (IsList(..))+#if !(MIN_VERSION_base(4,11,0))+import Data.Semigroup (Semigroup(..))+#endif++-- spine lazy, value strict list of doubles+data List+ = Nil+ | {-# UNPACK #-} !Double :! List+ deriving (Eq,Ord,Typeable,Data)++infixr 5 :!+++instance Semigroup List where+ Nil <> xs = xs+ (x :! xs) <> ys = x :! (xs <> ys)++instance Monoid List where+ mempty = Nil+ mappend = (<>)++instance IsList List where+ type Item List = Double+ toList Nil = []+ toList (a :! as) = a : Exts.toList as+ fromList [] = Nil+ fromList (a : as) = a :! Exts.fromList as++instance Show List where+ showsPrec d = showsPrec d . Exts.toList++instance Read List where+ readPrec = Exts.fromList <$> step readPrec++lmap :: (Double -> Double) -> List -> List+lmap f (a :! as) = f a :! lmap f as+lmap _ Nil = Nil+++-- | @Tower@ is an AD 'Mode' that calculates a tangent tower by forward AD, and provides fast 'diffsUU', 'diffsUF'+newtype TowerDouble = Tower { getTower :: List }+ deriving (Data, Typeable)++instance Show TowerDouble where+ showsPrec n (Tower as) = showParen (n > 10) $ showString "Tower " . showsPrec 11 as++-- Local combinators++zeroPad :: Num a => [a] -> [a]+zeroPad xs = xs ++ repeat 0+{-# INLINE zeroPad #-}++zeroPadF :: (Functor f, Num a) => [f a] -> [f a]+zeroPadF fxs@(fx:_) = fxs ++ repeat (0 <$ fx)+zeroPadF _ = error "zeroPadF :: empty list"+{-# INLINE zeroPadF #-}++lnull :: List -> Bool+lnull Nil = True+lnull _ = False++transposePadF :: (Foldable f, Functor f) => Double -> f List -> [f Double]+transposePadF pad fx+ | all lnull fx = []+ | otherwise = fmap headPad fx : transposePadF pad (drop1 <$> fx)+ where+ headPad Nil = pad+ headPad (x :! _) = x+ drop1 (_ :! xs) = xs+ drop1 xs = xs++d :: Num a => [a] -> a+d (_:da:_) = da+d _ = 0+{-# INLINE d #-}++dl :: List -> Double+dl (_ :! da :! _) = da+dl _ = 0+{-# INLINE dl #-}++d' :: Num a => [a] -> (a, a)+d' (a:da:_) = (a, da)+d' (a:_) = (a, 0)+d' _ = (0, 0)+{-# INLINE d' #-}++dl' :: List -> (Double, Double)+dl' (a:!da:!_) = (a, da)+dl' (a:!_) = (a, 0)+dl' _ = (0, 0)+{-# INLINE dl' #-}++tangents :: TowerDouble -> TowerDouble+tangents (Tower Nil) = Tower Nil+tangents (Tower (_ :! xs)) = Tower xs+{-# INLINE tangents #-}++truncated :: TowerDouble -> Bool+truncated (Tower Nil) = True+truncated _ = False+{-# INLINE truncated #-}++bundle :: Double -> TowerDouble -> TowerDouble+bundle a (Tower as) = Tower (a :! as)+{-# INLINE bundle #-}++withD :: (Double, Double) -> TowerDouble+withD (a, da) = Tower (a :! da :! Nil)+{-# INLINE withD #-}++apply :: (TowerDouble -> b) -> Double -> b+apply f a = f (Tower (a :! 1 :! Nil))+{-# INLINE apply #-}++getADTower :: TowerDouble -> [Double]+getADTower = Exts.toList . getTower+{-# INLINE getADTower #-}++tower :: [Double] -> TowerDouble+tower = Tower . Exts.fromList++primal :: TowerDouble -> Double+primal (Tower (x:!_)) = x+primal _ = 0++instance Mode TowerDouble where+ type Scalar TowerDouble = Double++ auto a = Tower (a :! Nil)++ isKnownZero (Tower Nil) = True+ isKnownZero (Tower (0 :! Nil)) = True+ isKnownZero _ = False++ asKnownConstant (Tower Nil) = Just 0+ asKnownConstant (Tower (a :! Nil)) = Just a+ asKnownConstant Tower {} = Nothing++ isKnownConstant (Tower Nil) = True+ isKnownConstant (Tower (_ :! Nil)) = True+ isKnownConstant Tower {} = False++ zero = Tower Nil++ a *^ Tower bs = Tower (lmap (a*) bs)++ Tower as ^* b = Tower (lmap (*b) as)++ Tower as ^/ b = Tower (lmap (/b) as)++infixr 6 <+>++(<+>) :: TowerDouble -> TowerDouble -> TowerDouble+Tower Nil <+> bs = bs+as <+> Tower Nil = as+Tower (a:!as) <+> Tower (b:!bs) = Tower (c:!cs) where+ c = a + b+ Tower cs = Tower as <+> Tower bs++instance Jacobian TowerDouble where+ type D TowerDouble = TowerDouble+ unary f dadb b = bundle (f (primal b)) (tangents b * dadb)+ lift1 f df b = bundle (f (primal b)) (tangents b * df b)+ lift1_ f df b = a where+ a = bundle (f (primal b)) (tangents b * df a b)++ binary f dadb dadc b c = bundle (f (primal b) (primal c)) (tangents b * dadb + tangents c * dadc)+ lift2 f df b c = bundle (f (primal b) (primal c)) 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++lzipWith :: (Double -> Double -> Double) -> List -> List -> List+lzipWith f (a :! as) (b :! bs) = f a b :! lzipWith f as bs+lzipWith _ _ _ = Nil++lsumProd3 :: List -> List -> List -> Double+lsumProd3 as0 bs0 cs0 = go as0 bs0 cs0 0 where+ go (a :! as) (b :! bs) (c :! cs) !acc = go as bs cs (a*b*c + acc)+ go _ _ _ acc = acc;++ltail :: List -> List+ltail (_ :! as) = as+ltail _ = error "ltail"++-- mul xs ys = [ sum [xs!!j * ys!!(k-j)*bin k j | j <- [0..k]] | k <- [0..] ]+-- adapted for efficiency and to handle finite lists xs, ys+mul:: TowerDouble -> TowerDouble -> TowerDouble+mul (Tower Nil) _ = Tower Nil+mul (Tower (a :! as)) (Tower bs) = Tower (convs' (1 :! Nil) (a :! Nil) as bs)+ where convs' _ _ _ Nil = Nil+ convs' ps ars as bs = lsumProd3 ps ars bs :!+ case as of+ Nil -> convs'' (next' ps) ars bs+ a:!as -> convs' (next ps) (a:!ars) as bs+ convs'' _ _ Nil = undefined -- convs'' never called with last argument empty+ convs'' _ _ (_:! Nil) = Nil+ convs'' ps ars (_:!bs) = lsumProd3 ps ars bs :! convs'' (next' ps) ars bs+ next xs = 1 :! lzipWith (+) xs (ltail xs) <> (1 :! Nil) -- next row in Pascal's triangle+ next' xs = lzipWith (+) xs (ltail xs) <> (1 :! Nil) -- end part of next row in Pascal's triangle++#define HEAD TowerDouble+#define BODY1(x)+#define BODY2(x,y)+#define NO_Bounded+#include <instances.h>
src/Numeric/AD/Internal/Type.hs view
@@ -3,7 +3,7 @@ {-# LANGUAGE GeneralizedNewtypeDeriving #-} ----------------------------------------------------------------------------- ---- |----- Copyright : (c) Edward Kmett 2010-2015+---- Copyright : (c) Edward Kmett 2010-2021 ---- License : BSD3 ---- Maintainer : ekmett@gmail.com ---- Stability : experimental@@ -24,6 +24,7 @@ instance Mode a => Mode (AD s a) where type Scalar (AD s a) = Scalar a isKnownConstant = isKnownConstant . runAD+ asKnownConstant = asKnownConstant . runAD isKnownZero = isKnownZero . runAD zero = AD zero auto = AD . auto
src/Numeric/AD/Jacobian.hs view
@@ -1,4 +1,3 @@-{-# LANGUAGE CPP #-} {-# LANGUAGE DefaultSignatures #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-}@@ -11,7 +10,7 @@ {-# LANGUAGE UndecidableInstances #-} ----------------------------------------------------------------------------- -- |--- Copyright : (c) Edward Kmett 2010-2015+-- Copyright : (c) Edward Kmett 2010-2021 -- License : BSD3 -- Maintainer : ekmett@gmail.com -- Stability : experimental
src/Numeric/AD/Jet.hs view
@@ -1,14 +1,10 @@-{-# LANGUAGE CPP #-} {-# LANGUAGE TypeOperators #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE FlexibleContexts #-}-#if __GLASGOW_HASKELL__ >= 707 {-# LANGUAGE DeriveDataTypeable #-}-{-# LANGUAGE StandaloneDeriving #-}-#endif ----------------------------------------------------------------------------- -- |--- Copyright : (c) Edward Kmett 2010-2015+-- Copyright : (c) Edward Kmett 2010-2021 -- License : BSD3 -- Maintainer : ekmett@gmail.com -- Stability : experimental@@ -20,23 +16,14 @@ , headJet , tailJet , jet+ , unjet ) where -#ifndef MIN_VERSION_base-#define MIN_VERSION_base(x,y,z) 1-#endif--#if __GLASGOW_HASKELL__ < 710-import Control.Applicative-import Data.Foldable-import Data.Traversable-import Data.Monoid-#endif-+import Data.Functor.Rep 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 --@@ -44,6 +31,7 @@ -- -- > a :- f a :- f (f a) :- f (f (f a)) :- ... data Jet f a = a :- Jet f (f a)+ deriving Typeable -- | Used to sidestep the need for UndecidableInstances. newtype Showable = Showable (Int -> String -> String)@@ -86,19 +74,8 @@ 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+unjet :: Representable f => Jet f a -> Cofree f a+unjet (a :- as) = a :< (unjet <$> undist as) where+ undist :: Representable f => Jet f (f a) -> f (Jet f a)+ undist (fa :- fas) = tabulate $ \i -> index fa i :- index (undist fas) i -jetTyCon :: TyCon-#if MIN_VERSION_base(4,4,0)-jetTyCon = mkTyCon3 "ad" "Numeric.AD.Internal.Jet" "Jet"-#else-jetTyCon = mkTyCon "Numeric.AD.Internal.Jet.Jet"-#endif-{-# NOINLINE jetTyCon #-}-#endif
src/Numeric/AD/Mode.hs view
@@ -1,16 +1,19 @@-{-# LANGUAGE CPP #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE DefaultSignatures #-} {-# LANGUAGE FunctionalDependencies #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE PatternGuards #-}+{-# LANGUAGE PatternSynonyms #-}+{-# LANGUAGE ViewPatterns #-} {-# LANGUAGE Rank2Types #-} {-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-} {-# LANGUAGE UndecidableInstances #-} ----------------------------------------------------------------------------- -- |--- Copyright : (c) Edward Kmett 2010-2015+-- Copyright : (c) Edward Kmett 2010-2021 -- License : BSD3 -- Maintainer : ekmett@gmail.com -- Stability : experimental@@ -22,6 +25,8 @@ ( -- * AD modes Mode(..)+ , pattern KnownZero+ , pattern Auto ) where import Numeric.Natural@@ -36,16 +41,23 @@ class (Num t, Num (Scalar t)) => Mode t where type Scalar t+ type Scalar t = t+ -- | allowed to return False for items with a zero derivative, but we'll give more NaNs than strictly necessary isKnownConstant :: t -> Bool isKnownConstant _ = False + asKnownConstant :: t -> Maybe (Scalar t)+ asKnownConstant _ = Nothing+ -- | allowed to return False for zero, but we give more NaN's than strictly necessary isKnownZero :: t -> Bool isKnownZero _ = False -- | Embed a constant auto :: Scalar t -> t+ default auto :: (Scalar t ~ t) => Scalar t -> t+ auto = id -- | Scalar-vector multiplication (*^) :: Scalar t -> t -> t@@ -64,114 +76,106 @@ zero :: t zero = auto 0 +pattern KnownZero :: Mode s => s+pattern KnownZero <- (isKnownZero -> True) where+ KnownZero = zero++pattern Auto :: Mode s => Scalar s -> s+pattern Auto n <- (asKnownConstant -> Just n) where+ Auto n = auto n+ instance Mode Double where- type Scalar Double = Double isKnownConstant _ = True+ asKnownConstant = Just isKnownZero x = 0 == x- auto = id (^/) = (/) instance Mode Float where- type Scalar Float = Float isKnownConstant _ = True+ asKnownConstant = Just isKnownZero x = 0 == x- auto = id (^/) = (/) instance Mode Int where- type Scalar Int = Int isKnownConstant _ = True+ asKnownConstant = Just isKnownZero x = 0 == x- auto = id (^/) = (/) instance Mode Integer where- type Scalar Integer = Integer isKnownConstant _ = True+ asKnownConstant = Just isKnownZero x = 0 == x- auto = id (^/) = (/) instance Mode Int8 where- type Scalar Int8 = Int8 isKnownConstant _ = True+ asKnownConstant = Just isKnownZero x = 0 == x- auto = id (^/) = (/) instance Mode Int16 where- type Scalar Int16 = Int16 isKnownConstant _ = True+ asKnownConstant = Just isKnownZero x = 0 == x- auto = id (^/) = (/) instance Mode Int32 where- type Scalar Int32 = Int32 isKnownConstant _ = True+ asKnownConstant = Just isKnownZero x = 0 == x- auto = id (^/) = (/) instance Mode Int64 where- type Scalar Int64 = Int64 isKnownConstant _ = True+ asKnownConstant = Just isKnownZero x = 0 == x- auto = id (^/) = (/) instance Mode Natural where- type Scalar Natural = Natural isKnownConstant _ = True+ asKnownConstant = Just isKnownZero x = 0 == x- auto = id (^/) = (/) instance Mode Word where- type Scalar Word = Word isKnownConstant _ = True+ asKnownConstant = Just isKnownZero x = 0 == x- auto = id (^/) = (/) instance Mode Word8 where- type Scalar Word8 = Word8 isKnownConstant _ = True+ asKnownConstant = Just isKnownZero x = 0 == x- auto = id (^/) = (/) instance Mode Word16 where- type Scalar Word16 = Word16 isKnownConstant _ = True+ asKnownConstant = Just isKnownZero x = 0 == x- auto = id (^/) = (/) instance Mode Word32 where- type Scalar Word32 = Word32 isKnownConstant _ = True+ asKnownConstant = Just isKnownZero x = 0 == x- auto = id (^/) = (/) instance Mode Word64 where- type Scalar Word64 = Word64 isKnownConstant _ = True+ asKnownConstant = Just isKnownZero x = 0 == x- auto = id (^/) = (/) instance RealFloat a => Mode (Complex a) where- type Scalar (Complex a) = Complex a isKnownConstant _ = True+ asKnownConstant = Just isKnownZero x = 0 == x- auto = id (^/) = (/) instance Integral a => Mode (Ratio a) where- type Scalar (Ratio a) = Ratio a isKnownConstant _ = True+ asKnownConstant = Just isKnownZero x = 0 == x- auto = id (^/) = (/)
+ src/Numeric/AD/Mode/Dense.hs view
@@ -0,0 +1,68 @@+{-# LANGUAGE Rank2Types #-}+-----------------------------------------------------------------------------+-- |+-- Copyright : (c) Edward Kmett 2010-2021+-- License : BSD3+-- Maintainer : ekmett@gmail.com+-- Stability : experimental+-- Portability : GHC only+--+-- First order dense forward mode using 'Traversable' functors+--+-----------------------------------------------------------------------------++module Numeric.AD.Mode.Dense+ ( AD, Dense, auto+ -- * Dense Gradients+ , grad+ , grad'+ , gradWith+ , gradWith'++ -- * Dense Jacobians (synonyms)+ , jacobian+ , jacobian'+ , jacobianWith+ , jacobianWith'+ ) where++import Numeric.AD.Internal.Dense (Dense)+import qualified Numeric.AD.Rank1.Dense as Rank1+import Numeric.AD.Internal.Type+import Numeric.AD.Mode++-- | The 'grad' function calculates the gradient of a non-scalar-to-scalar function with dense-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. f (AD s (Dense f a)) -> AD s (Dense f a)) -> f a -> f a+grad f = Rank1.grad (runAD.f.fmap AD)+{-# INLINE grad #-}++grad' :: (Traversable f, Num a) => (forall s. f (AD s (Dense f a)) -> AD s (Dense f a)) -> f a -> (a, f a)+grad' f = Rank1.grad' (runAD.f.fmap AD)+{-# INLINE grad' #-}++gradWith :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. f (AD s (Dense f a)) -> AD s (Dense f a)) -> f a -> f b+gradWith g f = Rank1.gradWith g (runAD.f.fmap AD)+{-# INLINE gradWith #-}++gradWith' :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. f (AD s (Dense f a)) -> AD s (Dense f a)) -> f a -> (a, f b)+gradWith' g f = Rank1.gradWith' g (runAD.f.fmap AD)+{-# INLINE gradWith' #-}++jacobian :: (Traversable f, Functor g, Num a) => (forall s. f (AD s (Dense f a)) -> g (AD s (Dense f a))) -> f a -> g (f a)+jacobian f = Rank1.jacobian (fmap runAD.f.fmap AD)+{-# INLINE jacobian #-}++jacobian' :: (Traversable f, Functor g, Num a) => (forall s. f (AD s (Dense f a)) -> g (AD s (Dense f a))) -> f a -> g (a, f a)+jacobian' f = Rank1.jacobian' (fmap runAD.f.fmap AD)+{-# INLINE jacobian' #-}++jacobianWith :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. f (AD s (Dense f a)) -> g (AD s (Dense f a))) -> f a -> g (f b)+jacobianWith g f = Rank1.jacobianWith g (fmap runAD.f.fmap AD)+{-# INLINE jacobianWith #-}++jacobianWith' :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. f (AD s (Dense f a)) -> g (AD s (Dense f a))) -> f a -> g (a, f b)+jacobianWith' g f = Rank1.jacobianWith' g (fmap runAD.f.fmap AD)+{-# INLINE jacobianWith' #-}
+ src/Numeric/AD/Mode/Dense/Representable.hs view
@@ -0,0 +1,79 @@+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE FlexibleContexts #-}+-----------------------------------------------------------------------------+-- |+-- Copyright : (c) Edward Kmett 2010-2021+-- License : BSD3+-- Maintainer : ekmett@gmail.com+-- Stability : experimental+-- Portability : GHC only+--+-- First order dense forward mode using 'Representable' functors+--+-----------------------------------------------------------------------------++module Numeric.AD.Mode.Dense.Representable+ ( AD, Repr, auto+ -- * Dense Gradients+ , grad+ , grad'+ , gradWith+ , gradWith'++ -- * Dense Jacobians (synonyms)+ , jacobian+ , jacobian'+ , jacobianWith+ , jacobianWith'+ ) where++import Data.Functor.Rep+import Numeric.AD.Internal.Dense.Representable (Repr)+import qualified Numeric.AD.Rank1.Dense.Representable as Rank1+import Numeric.AD.Internal.Type+import Numeric.AD.Mode++-- $setup+-- >>> :set -XDeriveGeneric -XDeriveFunctor+-- >>> import GHC.Generics (Generic1)+-- >>> import Data.Distributive (Distributive (..))+-- >>> import Data.Functor.Rep (Representable, distributeRep)+-- >>> data V3 a = V3 a a a deriving (Generic1, Functor, Show)+-- >>> instance Representable V3; instance Distributive V3 where distribute = distributeRep++-- | The 'grad' function calculates the gradient of a non-scalar-to-scalar function with dense-mode AD in a single pass.+--+-- >>> grad (\(V3 x y z) -> x*y+z) (V3 1 2 3)+-- V3 2 1 1+--+grad :: (Representable f, Eq (Rep f), Num a) => (forall s. f (AD s (Repr f a)) -> AD s (Repr f a)) -> f a -> f a+grad f = Rank1.grad (runAD.f.fmap AD)+{-# INLINE grad #-}++grad' :: (Representable f, Eq (Rep f), Num a) => (forall s. f (AD s (Repr f a)) -> AD s (Repr f a)) -> f a -> (a, f a)+grad' f = Rank1.grad' (runAD.f.fmap AD)+{-# INLINE grad' #-}++gradWith :: (Representable f, Eq (Rep f), Num a) => (a -> a -> b) -> (forall s. f (AD s (Repr f a)) -> AD s (Repr f a)) -> f a -> f b+gradWith g f = Rank1.gradWith g (runAD.f.fmap AD)+{-# INLINE gradWith #-}++gradWith' :: (Representable f, Eq (Rep f), Num a) => (a -> a -> b) -> (forall s. f (AD s (Repr f a)) -> AD s (Repr f a)) -> f a -> (a, f b)+gradWith' g f = Rank1.gradWith' g (runAD.f.fmap AD)+{-# INLINE gradWith' #-}++jacobian :: (Representable f, Eq (Rep f), Functor g, Num a) => (forall s. f (AD s (Repr f a)) -> g (AD s (Repr f a))) -> f a -> g (f a)+jacobian f = Rank1.jacobian (fmap runAD.f.fmap AD)+{-# INLINE jacobian #-}++jacobian' :: (Representable f, Eq (Rep f), Functor g, Num a) => (forall s. f (AD s (Repr f a)) -> g (AD s (Repr f a))) -> f a -> g (a, f a)+jacobian' f = Rank1.jacobian' (fmap runAD.f.fmap AD)+{-# INLINE jacobian' #-}++jacobianWith :: (Representable f, Eq (Rep f), Functor g, Num a) => (a -> a -> b) -> (forall s. f (AD s (Repr f a)) -> g (AD s (Repr f a))) -> f a -> g (f b)+jacobianWith g f = Rank1.jacobianWith g (fmap runAD.f.fmap AD)+{-# INLINE jacobianWith #-}++jacobianWith' :: (Representable f, Eq (Rep f), Functor g, Num a) => (a -> a -> b) -> (forall s. f (AD s (Repr f a)) -> g (AD s (Repr f a))) -> f a -> g (a, f b)+jacobianWith' g f = Rank1.jacobianWith' g (fmap runAD.f.fmap AD)+{-# INLINE jacobianWith' #-}
src/Numeric/AD/Mode/Forward.hs view
@@ -1,8 +1,7 @@-{-# LANGUAGE CPP #-} {-# LANGUAGE Rank2Types #-} ----------------------------------------------------------------------------- -- |--- Copyright : (c) Edward Kmett 2010-2015+-- Copyright : (c) Edward Kmett 2010-2021 -- License : BSD3 -- Maintainer : ekmett@gmail.com -- Stability : experimental@@ -44,9 +43,6 @@ , duF' ) where -#if __GLASGOW_HASKELL__ < 710-import Data.Traversable (Traversable)-#endif import Numeric.AD.Internal.Forward import Numeric.AD.Internal.On import Numeric.AD.Internal.Type@@ -150,21 +146,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 +168,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)])
src/Numeric/AD/Mode/Forward/Double.hs view
@@ -1,8 +1,7 @@-{-# LANGUAGE CPP #-} {-# LANGUAGE Rank2Types #-} ----------------------------------------------------------------------------- -- |--- Copyright : (c) Edward Kmett 2010-2015+-- Copyright : (c) Edward Kmett 2010-2021 -- License : BSD3 -- Maintainer : ekmett@gmail.com -- Stability : experimental@@ -41,9 +40,6 @@ , duF' ) where -#if __GLASGOW_HASKELL__ < 710-import Data.Traversable (Traversable)-#endif import Numeric.AD.Internal.Type (AD(AD), runAD) import Numeric.AD.Internal.Forward.Double (ForwardDouble) import qualified Numeric.AD.Rank1.Forward.Double as Rank1@@ -144,21 +140,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 +162,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)])
src/Numeric/AD/Mode/Kahn.hs view
@@ -1,14 +1,12 @@-{-# LANGUAGE CPP #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE FunctionalDependencies #-} {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE Rank2Types #-} {-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE TemplateHaskell #-} {-# LANGUAGE UndecidableInstances #-} ----------------------------------------------------------------------------- -- |--- Copyright : (c) Edward Kmett 2010-2015+-- Copyright : (c) Edward Kmett 2010-2021 -- License : BSD3 -- Maintainer : ekmett@gmail.com -- Stability : experimental@@ -46,9 +44,6 @@ , diffF' ) where -#if __GLASGOW_HASKELL__ < 710-import Data.Traversable (Traversable)-#endif import Numeric.AD.Internal.Kahn (Kahn) import Numeric.AD.Internal.On import Numeric.AD.Internal.Type (AD(..))
+ src/Numeric/AD/Mode/Kahn/Double.hs view
@@ -0,0 +1,240 @@+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE UndecidableInstances #-}+-----------------------------------------------------------------------------+-- |+-- Copyright : (c) Edward Kmett 2010-2021+-- License : BSD3+-- Maintainer : ekmett@gmail.com+-- Stability : experimental+-- Portability : GHC only+--+-- This module provides reverse-mode Automatic Differentiation using post-hoc linear time+-- topological sorting.+--+-- For reverse mode AD we use 'System.Mem.StableName.StableName' to recover sharing information from+-- the tape to avoid combinatorial explosion, and thus run asymptotically faster+-- than it could without such sharing information, but the use of side-effects+-- contained herein is benign.+--+-----------------------------------------------------------------------------++module Numeric.AD.Mode.Kahn.Double+ ( AD, Kahn, KahnDouble, auto+ -- * Gradient+ , grad+ , grad'+ , gradWith+ , gradWith'+ -- * Jacobian+ , jacobian+ , jacobian'+ , jacobianWith+ , jacobianWith'+ -- * Hessian+ , hessian+ , hessianF+ -- * Derivatives+ , diff+ , diff'+ , diffF+ , diffF'+ ) where++import Numeric.AD.Internal.Kahn (Kahn)+import Numeric.AD.Internal.Kahn.Double (KahnDouble)+import Numeric.AD.Internal.On+import Numeric.AD.Internal.Type (AD(..))+import Numeric.AD.Mode+import qualified Numeric.AD.Rank1.Kahn.Double as Rank1++-- $setup+--+-- >>> import Numeric.AD.Internal.Doctest++-- | 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.0,1.0,1.0]+--+-- >>> grad (\[x,y] -> x**y) [0,2]+-- [0.0,NaN]+grad+ :: Traversable f+ => (forall s. f (AD s KahnDouble) -> AD s KahnDouble)+ -> f Double+ -> f Double+grad f = Rank1.grad (runAD.f.fmap AD)+{-# INLINE grad #-}++-- | 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+ => (forall s. f (AD s KahnDouble) -> AD s KahnDouble)+ -> f Double+ -> (Double, f Double)+grad' f = Rank1.grad' (runAD.f.fmap AD)+{-# INLINE grad' #-}++-- | @'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@.+--+-- @+-- 'grad' = 'gradWith' (\_ dx -> dx)+-- 'id' = 'gradWith' const+-- @+--+--+gradWith+ :: Traversable f+ => (Double -> Double -> b)+ -> (forall s. f (AD s KahnDouble) -> AD s KahnDouble)+ -> f Double -> f b+gradWith g f = Rank1.gradWith g (runAD.f.fmap AD)+{-# INLINE gradWith #-}++-- | @'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+ => (Double -> Double -> b)+ -> (forall s. f (AD s KahnDouble) -> AD s KahnDouble)+ -> f Double -> (Double, f b)+gradWith' g f = Rank1.gradWith' g (runAD.f.fmap AD)+{-# INLINE gradWith' #-}++-- | 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.0,1.0],[1.0,0.0],[1.0,2.0]]+--+-- >>> 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)+ => (forall s. f (AD s KahnDouble) -> g (AD s KahnDouble))+ -> f Double -> g (f Double)+jacobian f = Rank1.jacobian (fmap runAD.f.fmap AD)+{-# INLINE jacobian #-}++-- | 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)+ => (forall s. f (AD s KahnDouble) -> g (AD s KahnDouble))+ -> f Double -> g (Double, f Double)+jacobian' f = Rank1.jacobian' (fmap runAD.f.fmap AD)+{-# INLINE jacobian' #-}++-- | '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@.+--+-- @+-- 'jacobian' = 'jacobianWith' (\_ dx -> dx)+-- 'jacobianWith' 'const' = (\f x -> 'const' x '<$>' f x)+-- @+jacobianWith+ :: (Traversable f, Functor g)+ => (Double -> Double -> b)+ -> (forall s. f (AD s KahnDouble) -> g (AD s KahnDouble))+ -> f Double -> g (f b)+jacobianWith g f = Rank1.jacobianWith g (fmap runAD.f.fmap AD)+{-# INLINE jacobianWith #-}++-- | '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)+ => (Double -> Double -> b) -> (forall s. f (AD s KahnDouble) -> g (AD s KahnDouble))+ -> f Double -> g (Double, f b)+jacobianWith' g f = Rank1.jacobianWith' g (fmap runAD.f.fmap AD)+{-# INLINE jacobianWith' #-}++-- | Compute the derivative of a function.+--+-- >>> diff sin 0+-- 1.0+--+-- >>> cos 0+-- 1.0+diff+ :: (forall s. AD s KahnDouble -> AD s KahnDouble) -> Double -> Double+diff f = Rank1.diff (runAD.f.AD)+{-# INLINE diff #-}++-- | The 'diff'' function calculates the value and derivative, as a+-- pair, of a scalar-to-scalar function.+--+--+-- >>> diff' sin 0+-- (0.0,1.0)+diff'+ :: (forall s. AD s KahnDouble -> AD s KahnDouble)+ -> Double+ -> (Double, Double)+diff' f = Rank1.diff' (runAD.f.AD)+{-# INLINE diff' #-}++-- | Compute the derivatives of a function that returns a vector with regards to its single input.+--+-- >>> diffF (\a -> [sin a, cos a]) 0+-- [1.0,0.0]+diffF+ :: Functor f+ => (forall s. AD s KahnDouble -> f (AD s KahnDouble))+ -> Double+ -> f Double+diffF f = Rank1.diffF (fmap runAD.f.AD)+{-# INLINE diffF #-}++-- | Compute the derivatives of a function that returns a vector with regards to its single input+-- as well as the primal answer.+--+-- >>> diffF' (\a -> [sin a, cos a]) 0+-- [(0.0,1.0),(1.0,0.0)]+diffF'+ :: Functor f+ => (forall s. AD s KahnDouble -> f (AD s KahnDouble))+ -> Double+ -> f (Double, Double)+diffF' f = Rank1.diffF' (fmap runAD.f.AD)+{-# INLINE diffF' #-}++-- | 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.0,1.0],[1.0,0.0]]+hessian+ :: Traversable f+ => (forall s. f (AD s (On (Kahn KahnDouble))) -> AD s (On (Kahn KahnDouble)))+ -> f Double -> f (f Double)+hessian f = Rank1.hessian (runAD.f.fmap AD)++-- | 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+ :: (Traversable f, Functor g)+ => (forall s. f (AD s (On (Kahn KahnDouble))) -> g (AD s (On (Kahn KahnDouble))))+ -> f Double -> g (f (f Double))+hessianF f = Rank1.hessianF (fmap runAD.f.fmap AD)
src/Numeric/AD/Mode/Reverse.hs view
@@ -1,6 +1,4 @@-{-# LANGUAGE CPP #-} {-# LANGUAGE Rank2Types #-}-{-# LANGUAGE TemplateHaskell #-} {-# LANGUAGE BangPatterns #-} {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE FunctionalDependencies #-}@@ -9,7 +7,7 @@ {-# LANGUAGE ScopedTypeVariables #-} ----------------------------------------------------------------------------- -- |--- Copyright : (c) Edward Kmett 2010-2015+-- Copyright : (c) Edward Kmett 2010-2021 -- License : BSD3 -- Maintainer : ekmett@gmail.com -- Stability : experimental@@ -45,10 +43,7 @@ , diffF' ) where -#if __GLASGOW_HASKELL__ < 710-import Data.Functor ((<$>))-import Data.Traversable (Traversable)-#endif+import Data.Typeable import Data.Functor.Compose import Data.Reflection (Reifies) import Numeric.AD.Internal.On@@ -67,8 +62,12 @@ -- -- >>> grad (\[x,y] -> x**y) [0,2] -- [0.0,NaN]-grad :: (Traversable f, Num a) => (forall s. Reifies s Tape => f (Reverse s a) -> Reverse s a) -> f a -> f a-grad f as = reifyTape (snd bds) $ \p -> unbind vs $! partialArrayOf p bds $! f vs where+grad+ :: (Traversable f, Num a)+ => (forall s. (Reifies s Tape, Typeable s) => f (Reverse s a) -> Reverse s a)+ -> f a+ -> f a+grad f as = reifyTypeableTape (snd bds) $ \p -> unbind vs $! partialArrayOf p bds $! f vs where (vs, bds) = bind as {-# INLINE grad #-} @@ -76,8 +75,12 @@ -- -- >>> grad' (\[x,y,z] -> x*y+z) [1,2,3] -- (5,[2,1,1])-grad' :: (Traversable f, Num a) => (forall s. Reifies s Tape => f (Reverse s a) -> Reverse s a) -> f a -> (a, f a)-grad' f as = reifyTape (snd bds) $ \p -> case f vs of+grad'+ :: (Traversable f, Num a)+ => (forall s. (Reifies s Tape, Typeable s) => f (Reverse s a) -> Reverse s a)+ -> f a+ -> (a, f a)+grad' f as = reifyTypeableTape (snd bds) $ \p -> case f vs of r -> (primal r, unbind vs $! partialArrayOf p bds $! r) where (vs, bds) = bind as {-# INLINE grad' #-}@@ -89,8 +92,13 @@ -- 'grad' == 'gradWith' (\_ dx -> dx) -- 'id' == 'gradWith' 'const' -- @-gradWith :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Reifies s Tape => f (Reverse s a) -> Reverse s a) -> f a -> f b-gradWith g f as = reifyTape (snd bds) $ \p -> unbindWith g vs $! partialArrayOf p bds $! f vs+gradWith+ :: (Traversable f, Num a)+ => (a -> a -> b)+ -> (forall s. (Reifies s Tape, Typeable s) => f (Reverse s a) -> Reverse s a)+ -> f a+ -> f b+gradWith g f as = reifyTypeableTape (snd bds) $ \p -> unbindWith g vs $! partialArrayOf p bds $! f vs where (vs,bds) = bind as {-# INLINE gradWith #-} @@ -100,8 +108,13 @@ -- @ -- 'grad'' == 'gradWith'' (\_ dx -> dx) -- @-gradWith' :: (Traversable f, Num a) => (a -> a -> b) -> (forall s. Reifies s Tape => f (Reverse s a) -> Reverse s a) -> f a -> (a, f b)-gradWith' g f as = reifyTape (snd bds) $ \p -> case f vs of+gradWith'+ :: (Traversable f, Num a)+ => (a -> a -> b)+ -> (forall s. (Reifies s Tape, Typeable s) => f (Reverse s a) -> Reverse s a)+ -> f a+ -> (a, f b)+gradWith' g f as = reifyTypeableTape (snd bds) $ \p -> case f vs of r -> (primal r, unbindWith g vs $! partialArrayOf p bds $! r) where (vs, bds) = bind as {-# INLINE gradWith' #-}@@ -110,8 +123,12 @@ -- -- >>> jacobian (\[x,y] -> [y,x,x*y]) [2,1] -- [[0,1],[1,0],[1,2]]-jacobian :: (Traversable f, Functor g, Num a) => (forall s. Reifies s Tape => f (Reverse s a) -> g (Reverse s a)) -> f a -> g (f a)-jacobian f as = reifyTape (snd bds) $ \p -> unbind vs . partialArrayOf p bds <$> f vs where+jacobian+ :: (Traversable f, Functor g, Num a)+ => (forall s. (Reifies s Tape, Typeable s) => f (Reverse s a) -> g (Reverse s a))+ -> f a+ -> g (f a)+jacobian f as = reifyTypeableTape (snd bds) $ \p -> unbind vs . partialArrayOf p bds <$> f vs where (vs, bds) = bind as {-# INLINE jacobian #-} @@ -121,8 +138,12 @@ -- -- >>> 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. Reifies s Tape => f (Reverse s a) -> g (Reverse s a)) -> f a -> g (a, f a)-jacobian' f as = reifyTape (snd bds) $ \p ->+jacobian'+ :: (Traversable f, Functor g, Num a)+ => (forall s. (Reifies s Tape, Typeable s) => f (Reverse s a) -> g (Reverse s a))+ -> f a+ -> g (a, f a)+jacobian' f as = reifyTypeableTape (snd bds) $ \p -> let row a = (primal a, unbind vs $! partialArrayOf p bds $! a) in row <$> f vs where (vs, bds) = bind as@@ -136,8 +157,13 @@ -- 'jacobian' == 'jacobianWith' (\_ dx -> dx) -- 'jacobianWith' 'const' == (\f x -> 'const' x '<$>' f x) -- @-jacobianWith :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. Reifies s Tape => f (Reverse s a) -> g (Reverse s a)) -> f a -> g (f b)-jacobianWith g f as = reifyTape (snd bds) $ \p -> unbindWith g vs . partialArrayOf p bds <$> f vs where+jacobianWith+ :: (Traversable f, Functor g, Num a)+ => (a -> a -> b)+ -> (forall s. (Reifies s Tape, Typeable s) => f (Reverse s a) -> g (Reverse s a))+ -> f a+ -> g (f b)+jacobianWith g f as = reifyTypeableTape (snd bds) $ \p -> unbindWith g vs . partialArrayOf p bds <$> f vs where (vs, bds) = bind as {-# INLINE jacobianWith #-} @@ -148,8 +174,13 @@ -- -- @'jacobian'' == 'jacobianWith'' (\_ dx -> dx)@ ---jacobianWith' :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (forall s. Reifies s Tape => f (Reverse s a) -> g (Reverse s a)) -> f a -> g (a, f b)-jacobianWith' g f as = reifyTape (snd bds) $ \p ->+jacobianWith'+ :: (Traversable f, Functor g, Num a)+ => (a -> a -> b)+ -> (forall s. (Reifies s Tape, Typeable s) => f (Reverse s a) -> g (Reverse s a))+ -> f a+ -> g (a, f b)+jacobianWith' g f as = reifyTypeableTape (snd bds) $ \p -> let row a = (primal a, unbindWith g vs $! partialArrayOf p bds $! a) in row <$> f vs where (vs, bds) = bind as@@ -159,8 +190,12 @@ -- -- >>> diff sin 0 -- 1.0-diff :: Num a => (forall s. Reifies s Tape => Reverse s a -> Reverse s a) -> a -> a-diff f a = reifyTape 1 $ \p -> derivativeOf p $! f (var a 0)+diff+ :: Num a+ => (forall s. (Reifies s Tape, Typeable s) => Reverse s a -> Reverse s a)+ -> a+ -> a+diff f a = reifyTypeableTape 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.@@ -170,25 +205,37 @@ -- -- >>> diff' exp 0 -- (1.0,1.0)-diff' :: Num a => (forall s. Reifies s Tape => Reverse s a -> Reverse s a) -> a -> (a, a)-diff' f a = reifyTape 1 $ \p -> derivativeOf' p $! f (var a 0)+diff'+ :: Num a+ => (forall s. (Reifies s Tape, Typeable s) => Reverse s a -> Reverse s a)+ -> a+ -> (a, a)+diff' f a = reifyTypeableTape 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, Num a) => (forall s. Reifies s Tape => Reverse s a -> f (Reverse s a)) -> a -> f a-diffF f a = reifyTape 1 $ \p -> derivativeOf p <$> f (var a 0)++diffF+ :: (Functor f, Num a)+ => (forall s. (Reifies s Tape, Typeable s) => Reverse s a -> f (Reverse s a))+ -> a+ -> f a+diffF f a = reifyTypeableTape 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, Num a) => (forall s. Reifies s Tape => Reverse s a -> f (Reverse s a)) -> a -> f (a, a)-diffF' f a = reifyTape 1 $ \p -> derivativeOf' p <$> f (var a 0)+diffF'+ :: (Functor f, Num a)+ => (forall s. (Reifies s Tape, Typeable s) => Reverse s a -> f (Reverse s a))+ -> a+ -> f (a, a)+diffF' f a = reifyTypeableTape 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.@@ -197,7 +244,14 @@ -- -- >>> hessian (\[x,y] -> x*y) [1,2] -- [[0,1],[1,0]]-hessian :: (Traversable f, Num a) => (forall s s'. (Reifies s Tape, Reifies s' Tape) => f (On (Reverse s (Reverse s' a))) -> (On (Reverse s (Reverse s' a)))) -> f a -> f (f a)+hessian+ :: (Traversable f, Num a)+ => ( forall s s'.+ (Reifies s Tape, Typeable s, Reifies s' Tape, Typeable s') =>+ f (On (Reverse s (Reverse s' a))) -> On (Reverse s (Reverse s' a))+ )+ -> f a+ -> f (f a) hessian f = jacobian (grad (off . f . fmap On)) {-# INLINE hessian #-} @@ -207,6 +261,13 @@ -- -- >>> hessianF (\[x,y] -> [x*y,x+y,exp x*cos y]) [1,2 :: RDouble] -- [[[0.0,1.0],[1.0,0.0]],[[0.0,0.0],[0.0,0.0]],[[-1.131204383757,-2.471726672005],[-2.471726672005,1.131204383757]]]-hessianF :: (Traversable f, Functor g, Num a) => (forall s s'. (Reifies s Tape, Reifies s' Tape) => f (On (Reverse s (Reverse s' a))) -> g (On (Reverse s (Reverse s' a)))) -> f a -> g (f (f a))+hessianF+ :: (Traversable f, Functor g, Num a)+ => (forall s s'. + (Reifies s Tape, Typeable s, Reifies s' Tape, Typeable s') =>+ f (On (Reverse s (Reverse s' a))) -> g (On (Reverse s (Reverse s' a)))+ )+ -> f a+ -> g (f (f a)) hessianF f = getCompose . jacobian (Compose . jacobian (fmap off . f . fmap On)) {-# INLINE hessianF #-}
+ src/Numeric/AD/Mode/Reverse/Double.hs view
@@ -0,0 +1,274 @@+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE ScopedTypeVariables #-}+-----------------------------------------------------------------------------+-- |+-- Copyright : (c) Edward Kmett 2010-2021+-- 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 to be unboxed.+--+-----------------------------------------------------------------------------++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++import Data.Typeable+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, Typeable s) => f (ReverseDouble s) -> ReverseDouble s)+ -> f Double+ -> f Double+grad f as = reifyTypeableTape (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, Typeable s) => f (ReverseDouble s) -> ReverseDouble s)+ -> f Double+ -> (Double, f Double)+grad' f as = reifyTypeableTape (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, Typeable s) => f (ReverseDouble s) -> ReverseDouble s)+ -> f Double+ -> f b+gradWith g f as = reifyTypeableTape (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, Typeable s) => f (ReverseDouble s) -> ReverseDouble s)+ -> f Double+ -> (Double, f b)+gradWith' g f as = reifyTypeableTape (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, Typeable s) => f (ReverseDouble s) -> g (ReverseDouble s))+ -> f Double+ -> g (f Double)+jacobian f as = reifyTypeableTape (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, Typeable s) => f (ReverseDouble s) -> g (ReverseDouble s))+ -> f Double+ -> g (Double, f Double)+jacobian' f as = reifyTypeableTape (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, Typeable s) => f (ReverseDouble s) -> g (ReverseDouble s))+ -> f Double+ -> g (f b)+jacobianWith g f as = reifyTypeableTape (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, Typeable s) => f (ReverseDouble s) -> g (ReverseDouble s))+ -> f Double+ -> g (Double, f b)+jacobianWith' g f as = reifyTypeableTape (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, Typeable s) => ReverseDouble s -> ReverseDouble s)+ -> Double+ -> Double+diff f a = reifyTypeableTape 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, Typeable s) => ReverseDouble s -> ReverseDouble s)+ -> Double+ -> (Double, Double)+diff' f a = reifyTypeableTape 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, Typeable s) => ReverseDouble s -> f (ReverseDouble s))+ -> Double+ -> f Double+diffF f a = reifyTypeableTape 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, Typeable s) => ReverseDouble s -> f (ReverseDouble s))+ -> Double+ -> f (Double, Double)+diffF' f a = reifyTypeableTape 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, Typeable s, Reifies s' Tape, Typeable s') =>+ 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, Typeable s, Reifies s' Tape, Typeable s') =>+ 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 #-}
src/Numeric/AD/Mode/Sparse.hs view
@@ -1,8 +1,7 @@-{-# LANGUAGE CPP #-} {-# LANGUAGE Rank2Types #-} ----------------------------------------------------------------------------- -- |--- Copyright : (c) Edward Kmett 2010-2015+-- Copyright : (c) Edward Kmett 2010-2021 -- License : BSD3 -- Maintainer : ekmett@gmail.com -- Stability : experimental@@ -37,9 +36,6 @@ ) where import Control.Comonad.Cofree (Cofree)-#if __GLASGOW_HASKELL__ < 710-import Data.Traversable (Traversable)-#endif import Numeric.AD.Internal.Sparse (Sparse) import qualified Numeric.AD.Rank1.Sparse as Rank1 import Numeric.AD.Internal.Type
+ src/Numeric/AD/Mode/Sparse/Double.hs view
@@ -0,0 +1,162 @@+{-# LANGUAGE Rank2Types #-}+-----------------------------------------------------------------------------+-- |+-- Copyright : (c) Edward Kmett 2010-2021+-- License : BSD3+-- Maintainer : ekmett@gmail.com+-- Stability : experimental+-- Portability : GHC only+--+-- Higher order derivatives via a \"dual number tower\".+--+-----------------------------------------------------------------------------++module Numeric.AD.Mode.Sparse.Double+ ( AD, SparseDouble, auto+ -- * Sparse Gradients+ , grad+ , grad'+ , grads+ , gradWith+ , gradWith'++ -- * Sparse Jacobians (synonyms)+ , jacobian+ , jacobian'+ , jacobianWith+ , jacobianWith'+ , jacobians++ -- * Sparse Hessians+ , hessian+ , hessian'++ , hessianF+ , hessianF'+ ) where++import Control.Comonad.Cofree (Cofree)+import Numeric.AD.Internal.Sparse.Double (SparseDouble)+import qualified Numeric.AD.Rank1.Sparse.Double as Rank1+import Numeric.AD.Internal.Type+import Numeric.AD.Mode++-- | The 'grad' function calculates the gradient of a non-scalar-to-scalar function with sparse-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. f (AD s SparseDouble) -> AD s SparseDouble)+ -> f Double+ -> f Double+grad f = Rank1.grad (runAD.f.fmap AD)+{-# INLINE grad #-}++grad'+ :: Traversable f+ => (forall s. f (AD s SparseDouble) -> AD s SparseDouble)+ -> f Double+ -> (Double, f Double)+grad' f = Rank1.grad' (runAD.f.fmap AD)+{-# INLINE grad' #-}++gradWith+ :: Traversable f+ => (Double -> Double -> b)+ -> (forall s. f (AD s SparseDouble) -> AD s SparseDouble)+ -> f Double+ -> f b+gradWith g f = Rank1.gradWith g (runAD.f.fmap AD)+{-# INLINE gradWith #-}++gradWith'+ :: Traversable f+ => (Double -> Double -> b)+ -> (forall s. f (AD s SparseDouble) -> AD s SparseDouble)+ -> f Double+ -> (Double, f b)+gradWith' g f = Rank1.gradWith' g (runAD.f.fmap AD)+{-# INLINE gradWith' #-}++jacobian+ :: (Traversable f, Functor g)+ => (forall s. f (AD s SparseDouble) -> g (AD s SparseDouble))+ -> f Double -> g (f Double)+jacobian f = Rank1.jacobian (fmap runAD.f.fmap AD)+{-# INLINE jacobian #-}++jacobian'+ :: (Traversable f, Functor g)+ => (forall s. f (AD s SparseDouble) -> g (AD s SparseDouble))+ -> f Double+ -> g (Double, f Double)+jacobian' f = Rank1.jacobian' (fmap runAD.f.fmap AD)+{-# INLINE jacobian' #-}++jacobianWith+ :: (Traversable f, Functor g)+ => (Double -> Double -> b)+ -> (forall s. f (AD s SparseDouble) -> g (AD s SparseDouble))+ -> f Double+ -> g (f b)+jacobianWith g f = Rank1.jacobianWith g (fmap runAD.f.fmap AD)+{-# INLINE jacobianWith #-}++jacobianWith'+ :: (Traversable f, Functor g)+ => (Double -> Double -> b)+ -> (forall s. f (AD s SparseDouble) -> g (AD s SparseDouble))+ -> f Double+ -> g (Double, f b)+jacobianWith' g f = Rank1.jacobianWith' g (fmap runAD.f.fmap AD)+{-# INLINE jacobianWith' #-}++grads+ :: Traversable f+ => (forall s. f (AD s SparseDouble) -> AD s SparseDouble)+ -> f Double -> Cofree f Double+grads f = Rank1.grads (runAD.f.fmap AD)+{-# INLINE grads #-}++jacobians+ :: (Traversable f, Functor g)+ => (forall s. f (AD s SparseDouble) -> g (AD s SparseDouble))+ -> f Double+ -> g (Cofree f Double)+jacobians f = Rank1.jacobians (fmap runAD.f.fmap AD)+{-# INLINE jacobians #-}++hessian+ :: Traversable f+ => (forall s. f (AD s SparseDouble) -> AD s SparseDouble)+ -> f Double+ -> f (f Double)+hessian f = Rank1.hessian (runAD.f.fmap AD)+{-# INLINE hessian #-}++hessian'+ :: Traversable f+ => (forall s. f (AD s SparseDouble) -> AD s SparseDouble)+ -> f Double -> (Double, f (Double, f Double))+hessian' f = Rank1.hessian' (runAD.f.fmap AD)+{-# INLINE hessian' #-}++hessianF+ :: (Traversable f, Functor g)+ => (forall s. f (AD s SparseDouble) -> g (AD s SparseDouble))+ -> f Double -> g (f (f Double))+hessianF f = Rank1.hessianF (fmap runAD.f.fmap AD)+{-# INLINE hessianF #-}++hessianF'+ :: (Traversable f, Functor g)+ => (forall s. f (AD s SparseDouble) -> g (AD s SparseDouble))+ -> f Double+ -> g (Double, f (Double, f Double))+hessianF' f = Rank1.hessianF' (fmap runAD.f.fmap AD)+{-# INLINE hessianF' #-}
src/Numeric/AD/Mode/Tower.hs view
@@ -1,7 +1,7 @@ {-# LANGUAGE Rank2Types #-} ----------------------------------------------------------------------------- -- |--- Copyright : (c) Edward Kmett 2010-2015+-- Copyright : (c) Edward Kmett 2010-2021 -- License : BSD3 -- Maintainer : ekmett@gmail.com -- Stability : experimental@@ -42,7 +42,6 @@ import Numeric.AD.Internal.Tower (Tower) import Numeric.AD.Internal.Type (AD(..)) import Numeric.AD.Mode- diffs :: Num a => (forall s. AD s (Tower a) -> AD s (Tower a)) -> a -> [a] diffs f = Rank1.diffs (runAD.f.AD)
+ src/Numeric/AD/Mode/Tower/Double.hs view
@@ -0,0 +1,115 @@+{-# LANGUAGE Rank2Types #-}+-----------------------------------------------------------------------------+-- |+-- Copyright : (c) Edward Kmett 2010-2021+-- License : BSD3+-- Maintainer : ekmett@gmail.com+-- Stability : experimental+-- Portability : GHC only+--+-- Higher order derivatives via a \"dual number tower\".+--+-----------------------------------------------------------------------------+module Numeric.AD.Mode.Tower.Double+ ( AD+ , TowerDouble+ , auto+ -- * Taylor Series+ , taylor+ , taylor0+ -- * Maclaurin Series+ , maclaurin+ , maclaurin0+ -- * Derivatives+ , diff -- first derivative of (Double -> a)+ , diff' -- answer and first derivative of (Double -> a)+ , diffs -- answer and all derivatives of (Double -> a)+ , diffs0 -- zero padded derivatives of (Double -> a)+ , diffsF -- answer and all derivatives of (Double -> f a)+ , diffs0F -- zero padded derivatives of (Double -> f a)+ -- * Directional Derivatives+ , du -- directional derivative of (Double -> a)+ , du' -- answer and directional derivative of (Double -> a)+ , dus -- answer and all directional derivatives of (Double -> a)+ , dus0 -- answer and all zero padded directional derivatives of (Double -> a)+ , duF -- directional derivative of (Double -> f a)+ , duF' -- answer and directional derivative of (Double -> f a)+ , dusF -- answer and all directional derivatives of (Double -> f a)+ , dus0F -- answer and all zero padded directional derivatives of (Double -> a)+ ) where++import qualified Numeric.AD.Rank1.Tower.Double as Rank1+import Numeric.AD.Internal.Tower.Double (TowerDouble)+import Numeric.AD.Internal.Type (AD(..))+import Numeric.AD.Mode++diffs :: (forall s. AD s TowerDouble -> AD s TowerDouble) -> Double -> [Double]+diffs f = Rank1.diffs (runAD.f.AD)+{-# INLINE diffs #-}++diffs0 :: (forall s. AD s TowerDouble -> AD s TowerDouble) -> Double -> [Double]+diffs0 f = Rank1.diffs0 (runAD.f.AD)+{-# INLINE diffs0 #-}++diffsF :: Functor f => (forall s. AD s TowerDouble -> f (AD s TowerDouble)) -> Double -> f [Double]+diffsF f = Rank1.diffsF (fmap runAD.f.AD)+{-# INLINE diffsF #-}++diffs0F :: Functor f => (forall s. AD s TowerDouble -> f (AD s TowerDouble)) -> Double -> f [Double]+diffs0F f = Rank1.diffs0F (fmap runAD.f.AD)+{-# INLINE diffs0F #-}++taylor :: (forall s. AD s TowerDouble -> AD s TowerDouble) -> Double -> Double -> [Double]+taylor f = Rank1.taylor (runAD.f.AD)++taylor0 :: (forall s. AD s TowerDouble -> AD s TowerDouble) -> Double -> Double -> [Double]+taylor0 f = Rank1.taylor0 (runAD.f.AD)+{-# INLINE taylor0 #-}++maclaurin :: (forall s. AD s TowerDouble -> AD s TowerDouble) -> Double -> [Double]+maclaurin f = Rank1.maclaurin (runAD.f.AD)+{-# INLINE maclaurin #-}++maclaurin0 :: (forall s. AD s TowerDouble -> AD s TowerDouble) -> Double -> [Double]+maclaurin0 f = Rank1.maclaurin0 (runAD.f.AD)+{-# INLINE maclaurin0 #-}++diff :: (forall s. AD s TowerDouble -> AD s TowerDouble) -> Double -> Double+diff f = Rank1.diff (runAD.f.AD)+{-# INLINE diff #-}++diff' :: (forall s. AD s TowerDouble -> AD s TowerDouble) -> Double -> (Double, Double)+diff' f = Rank1.diff' (runAD.f.AD)+{-# INLINE diff' #-}++du :: Functor f => (forall s. f (AD s TowerDouble) -> AD s TowerDouble) -> f (Double, Double) -> Double+du f = Rank1.du (runAD.f. fmap AD)+{-# INLINE du #-}++du' :: Functor f => (forall s. f (AD s TowerDouble) -> AD s TowerDouble) -> f (Double, Double) -> (Double, Double)+du' f = Rank1.du' (runAD.f.fmap AD)+{-# INLINE du' #-}++duF :: (Functor f, Functor g) => (forall s. f (AD s TowerDouble) -> g (AD s TowerDouble)) -> f (Double, Double) -> g Double+duF f = Rank1.duF (fmap runAD.f.fmap AD)+{-# INLINE duF #-}++duF' :: (Functor f, Functor g) => (forall s. f (AD s TowerDouble) -> g (AD s TowerDouble)) -> f (Double, Double) -> g (Double, Double)+duF' f = Rank1.duF' (fmap runAD.f.fmap AD)+{-# INLINE duF' #-}++dus :: Functor f => (forall s. f (AD s TowerDouble) -> AD s TowerDouble) -> f [Double] -> [Double]+dus f = Rank1.dus (runAD.f.fmap AD)+{-# INLINE dus #-}++dus0 :: Functor f => (forall s. f (AD s TowerDouble) -> AD s TowerDouble) -> f [Double] -> [Double]+dus0 f = Rank1.dus0 (runAD.f.fmap AD)+{-# INLINE dus0 #-}++dusF :: (Functor f, Functor g) => (forall s. f (AD s TowerDouble) -> g (AD s TowerDouble)) -> f [Double] -> g [Double]+dusF f = Rank1.dusF (fmap runAD.f.fmap AD)+{-# INLINE dusF #-}++dus0F :: (Functor f, Functor g) => (forall s. f (AD s TowerDouble) -> g (AD s TowerDouble)) -> f [Double] -> g [Double]+dus0F f = Rank1.dus0F (fmap runAD.f.fmap AD)+{-# INLINE dus0F #-}
src/Numeric/AD/Newton.hs view
@@ -1,4 +1,3 @@-{-# LANGUAGE CPP #-} {-# LANGUAGE RankNTypes #-} {-# LANGUAGE GADTs #-} {-# LANGUAGE BangPatterns #-}@@ -11,7 +10,7 @@ {-# LANGUAGE ParallelListComp #-} ----------------------------------------------------------------------------- -- |--- Copyright : (c) Edward Kmett 2010-2015+-- Copyright : (c) Edward Kmett 2010-2021 -- License : BSD3 -- Maintainer : ekmett@gmail.com -- Stability : experimental@@ -38,11 +37,7 @@ , stochasticGradientDescent ) where -#if __GLASGOW_HASKELL__ < 710-import Data.Foldable (Foldable, all, sum)-#else import Data.Foldable (all, sum)-#endif import Data.Reflection (Reifies) import Data.Traversable import Numeric.AD.Internal.Combinators@@ -261,8 +256,8 @@ -- It uses reverse mode automatic differentiation to compute the gradient -- The learning rate is constant through out, and is set to 0.001 stochasticGradientDescent :: (Traversable f, Fractional a, Ord a)- => (forall s. Reifies s Tape => f (Scalar a) -> f (Reverse s a) -> Reverse s a)- -> [f (Scalar a)]+ => (forall s. Reifies s Tape => e -> f (Reverse s a) -> Reverse s a)+ -> [e] -> f a -> [f a] stochasticGradientDescent errorSingle d0 x0 = go xgx0 0.001 dLeft
src/Numeric/AD/Newton/Double.hs view
@@ -5,7 +5,7 @@ {-# LANGUAGE TypeFamilies #-} ----------------------------------------------------------------------------- -- |--- Copyright : (c) Edward Kmett 2015+-- Copyright : (c) Edward Kmett 2015-2021 -- License : BSD3 -- Maintainer : ekmett@gmail.com -- Stability : experimental@@ -17,13 +17,9 @@ ( -- * Newton's Method (Forward AD) findZero- , findZeroNoEq , inverse- , inverseNoEq , fixedPoint- , fixedPointNoEq , extremum- , extremumNoEq -- * Gradient Ascent/Descent (Reverse AD) , conjugateGradientDescent , conjugateGradientAscent@@ -38,7 +34,7 @@ import Numeric.AD.Internal.Or import Numeric.AD.Internal.Type (AD(..)) import Numeric.AD.Mode-import Numeric.AD.Rank1.Kahn as Kahn (Kahn, grad)+import Numeric.AD.Rank1.Kahn.Double as Kahn (KahnDouble, grad) import qualified Numeric.AD.Rank1.Newton.Double as Rank1 import Prelude hiding (all, mapM, sum) @@ -55,12 +51,6 @@ findZero f = Rank1.findZero (runAD.f.AD) {-# INLINE findZero #-} --- | The 'findZeroNoEq' function behaves the same as 'findZero' except that it--- doesn't truncate the list once the results become constant.-findZeroNoEq :: (forall s. AD s ForwardDouble -> AD s ForwardDouble) -> Double -> [Double]-findZeroNoEq f = Rank1.findZeroNoEq (runAD.f.AD)-{-# INLINE findZeroNoEq #-}- -- | The 'inverse' function inverts a scalar function using -- Newton's method; its output is a stream of increasingly accurate -- results. (Modulo the usual caveats.) If the stream becomes@@ -74,12 +64,6 @@ inverse f = Rank1.inverse (runAD.f.AD) {-# INLINE inverse #-} --- | The 'inverseNoEq' function behaves the same as 'inverse' except that it--- doesn't truncate the list once the results become constant.-inverseNoEq :: (forall s. AD s ForwardDouble -> AD s ForwardDouble) -> Double -> Double -> [Double]-inverseNoEq f = Rank1.inverseNoEq (runAD.f.AD)-{-# INLINE inverseNoEq #-}- -- | The 'fixedPoint' function find a fixedpoint of a scalar -- function using Newton's method; its output is a stream of -- increasingly accurate results. (Modulo the usual caveats.)@@ -93,12 +77,6 @@ fixedPoint f = Rank1.fixedPoint (runAD.f.AD) {-# INLINE fixedPoint #-} --- | The 'fixedPointNoEq' function behaves the same as 'fixedPoint' except that--- doesn't truncate the list once the results become constant.-fixedPointNoEq :: (forall s. AD s ForwardDouble -> AD s ForwardDouble) -> Double -> [Double]-fixedPointNoEq f = Rank1.fixedPointNoEq (runAD.f.AD)-{-# INLINE fixedPointNoEq #-}- -- | The 'extremum' function finds an extremum of a scalar -- function using Newton's method; produces a stream of increasingly -- accurate results. (Modulo the usual caveats.) If the stream@@ -110,12 +88,6 @@ extremum f = Rank1.extremum (runAD.f.AD) {-# INLINE extremum #-} --- | The 'extremumNoEq' function behaves the same as 'extremum' except that it--- doesn't truncate the list once the results become constant.-extremumNoEq :: (forall s. AD s (On (Forward ForwardDouble)) -> AD s (On (Forward ForwardDouble))) -> Double -> [Double]-extremumNoEq f = Rank1.extremumNoEq (runAD.f.AD)-{-# INLINE extremumNoEq #-}- -- | 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@@ -126,7 +98,7 @@ -- True conjugateGradientDescent :: Traversable f- => (forall s. Chosen s => f (Or s (On (Forward ForwardDouble)) (Kahn Double)) -> Or s (On (Forward ForwardDouble)) (Kahn Double))+ => (forall s. Chosen s => f (Or s (On (Forward ForwardDouble)) KahnDouble) -> Or s (On (Forward ForwardDouble)) KahnDouble) -> f Double -> [f Double] conjugateGradientDescent f = conjugateGradientAscent (negate . f) {-# INLINE conjugateGradientDescent #-}@@ -140,7 +112,7 @@ -- | Perform a conjugate gradient ascent using reverse mode automatic differentiation to compute the gradient. conjugateGradientAscent :: Traversable f- => (forall s. Chosen s => f (Or s (On (Forward ForwardDouble)) (Kahn Double)) -> Or s (On (Forward ForwardDouble)) (Kahn Double))+ => (forall s. Chosen s => f (Or s (On (Forward ForwardDouble)) KahnDouble) -> Or s (On (Forward ForwardDouble)) KahnDouble) -> f Double -> [f Double] conjugateGradientAscent f x0 = takeWhile (all (\a -> a == a)) (go x0 d0 d0 delta0) where
+ src/Numeric/AD/Rank1/Dense.hs view
@@ -0,0 +1,102 @@+-----------------------------------------------------------------------------+-- |+-- Copyright : (c) Edward Kmett 2010-2021+-- License : BSD3+-- Maintainer : ekmett@gmail.com+-- Stability : experimental+-- Portability : GHC only+--+-----------------------------------------------------------------------------++module Numeric.AD.Rank1.Dense+ ( Dense+ , auto+ -- * Sparse Gradients+ , grad+ , grad'+ , gradWith+ , gradWith'++ -- * Sparse Jacobians (synonyms)+ , jacobian+ , jacobian'+ , jacobianWith+ , jacobianWith'++ ) where++import Numeric.AD.Internal.Dense+import Numeric.AD.Internal.Combinators+import Numeric.AD.Mode++second :: (a -> b) -> (c, a) -> (c, b)+second g (a,b) = (a, g b)+{-# INLINE second #-}++grad+ :: (Traversable f, Num a)+ => (f (Dense f a) -> Dense f a)+ -> f a+ -> f a+grad f as = ds (0 <$ as) $ apply f as+{-# INLINE grad #-}++grad'+ :: (Traversable f, Num a)+ => (f (Dense f a) -> Dense f a)+ -> f a+ -> (a, f a)+grad' f as = ds' (0 <$ as) $ apply f as+{-# INLINE grad' #-}++gradWith+ :: (Traversable f, Num a)+ => (a -> a -> b)+ -> (f (Dense f a) -> Dense f a)+ -> f a+ -> f b+gradWith g f as = zipWithT g as $ grad f as+{-# INLINE gradWith #-}++gradWith'+ :: (Traversable f, Num a)+ => (a -> a -> b)+ -> (f (Dense f a) -> Dense f a)+ -> f a+ -> (a, f b)+gradWith' g f as = second (zipWithT g as) $ grad' f as+{-# INLINE gradWith' #-}++jacobian+ :: (Traversable f, Functor g, Num a)+ => (f (Dense f a) -> g (Dense f a))+ -> f a+ -> g (f a)+jacobian f as = ds (0 <$ as) <$> apply f as+{-# INLINE jacobian #-}++jacobian'+ :: (Traversable f, Functor g, Num a)+ => (f (Dense f a) -> g (Dense f a))+ -> f a+ -> g (a, f a)+jacobian' f as = ds' (0 <$ as) <$> apply f as+{-# INLINE jacobian' #-}++jacobianWith+ :: (Traversable f, Functor g, Num a)+ => (a -> a -> b)+ -> (f (Dense f a) -> g (Dense f a))+ -> f a+ -> g (f b)+jacobianWith g f as = zipWithT g as <$> jacobian f as+{-# INLINE jacobianWith #-}++jacobianWith'+ :: (Traversable f, Functor g, Num a)+ => (a -> a -> b)+ -> (f (Dense f a) -> g (Dense f a))+ -> f a+ -> g (a, f b)+jacobianWith' g f as = second (zipWithT g as) <$> jacobian' f as+{-# INLINE jacobianWith' #-}
+ src/Numeric/AD/Rank1/Dense/Representable.hs view
@@ -0,0 +1,103 @@+{-# LANGUAGE FlexibleContexts #-}+-----------------------------------------------------------------------------+-- |+-- Copyright : (c) Edward Kmett 2010-2021+-- License : BSD3+-- Maintainer : ekmett@gmail.com+-- Stability : experimental+-- Portability : GHC only+--+-- Dense forward mode automatic differentiation with representable functors.+--+-----------------------------------------------------------------------------++module Numeric.AD.Rank1.Dense.Representable+ ( Repr+ , auto+ -- * Sparse Gradients+ , grad+ , grad'+ , gradWith+ , gradWith'+ -- * Sparse Jacobians (synonyms)+ , jacobian+ , jacobian'+ , jacobianWith+ , jacobianWith'+ ) where++import Data.Functor.Rep+import Numeric.AD.Internal.Dense.Representable+import Numeric.AD.Mode++second :: (a -> b) -> (c, a) -> (c, b)+second g (a,b) = (a, g b)+{-# INLINE second #-}++grad+ :: (Representable f, Eq (Rep f), Num a)+ => (f (Repr f a) -> Repr f a)+ -> f a+ -> f a+grad f as = ds (pureRep 0) $ apply f as+{-# INLINE grad #-}++grad'+ :: (Representable f, Eq (Rep f), Num a)+ => (f (Repr f a) -> Repr f a)+ -> f a+ -> (a, f a)+grad' f as = ds' (pureRep 0) $ apply f as+{-# INLINE grad' #-}++gradWith+ :: (Representable f, Eq (Rep f), Num a)+ => (a -> a -> b)+ -> (f (Repr f a) -> Repr f a)+ -> f a+ -> f b+gradWith g f as = liftR2 g as $ grad f as+{-# INLINE gradWith #-}++gradWith'+ :: (Representable f, Eq (Rep f), Num a)+ => (a -> a -> b)+ -> (f (Repr f a) -> Repr f a)+ -> f a+ -> (a, f b)+gradWith' g f as = second (liftR2 g as) $ grad' f as+{-# INLINE gradWith' #-}++jacobian+ :: (Representable f, Eq (Rep f), Functor g, Num a)+ => (f (Repr f a) -> g (Repr f a))+ -> f a+ -> g (f a)+jacobian f as = ds (0 <$ as) <$> apply f as+{-# INLINE jacobian #-}++jacobian'+ :: (Representable f, Eq (Rep f), Functor g, Num a)+ => (f (Repr f a) -> g (Repr f a))+ -> f a+ -> g (a, f a)+jacobian' f as = ds' (0 <$ as) <$> apply f as+{-# INLINE jacobian' #-}++jacobianWith+ :: (Representable f, Eq (Rep f), Functor g, Num a)+ => (a -> a -> b)+ -> (f (Repr f a) -> g (Repr f a))+ -> f a+ -> g (f b)+jacobianWith g f as = liftR2 g as <$> jacobian f as+{-# INLINE jacobianWith #-}++jacobianWith'+ :: (Representable f, Eq (Rep f), Functor g, Num a)+ => (a -> a -> b)+ -> (f (Repr f a) -> g (Repr f a))+ -> f a+ -> g (a, f b)+jacobianWith' g f as = second (liftR2 g as) <$> jacobian' f as+{-# INLINE jacobianWith' #-}
src/Numeric/AD/Rank1/Forward.hs view
@@ -1,7 +1,6 @@-{-# LANGUAGE CPP #-} ----------------------------------------------------------------------------- -- |--- Copyright : (c) Edward Kmett 2010-2015+-- Copyright : (c) Edward Kmett 2010-2021 -- License : BSD3 -- Maintainer : ekmett@gmail.com -- Stability : experimental@@ -42,32 +41,43 @@ , duF' ) where -#if __GLASGOW_HASKELL__ < 710-import Data.Traversable (Traversable)-import Control.Applicative-#endif import Numeric.AD.Internal.Forward import Numeric.AD.Internal.On import Numeric.AD.Mode - -- | Compute the directional derivative of a function given a zipped up 'Functor' of the input values and their derivatives-du :: (Functor f, Num a) => (f (Forward a) -> Forward a) -> f (a, a) -> a+du+ :: (Functor f, Num a)+ => (f (Forward a) -> Forward a)+ -> 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) => (f (Forward a) -> Forward a) -> f (a, a) -> (a, a)+du'+ :: (Functor f, Num a)+ => (f (Forward a) -> Forward a)+ -> 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) => (f (Forward a) -> g (Forward a)) -> f (a, a) -> g a+duF+ :: (Functor f, Functor g, Num a)+ => (f (Forward a) -> g (Forward a))+ -> 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) => (f (Forward a) -> g (Forward a)) -> f (a, a) -> g (a, a)+duF'+ :: (Functor f, Functor g, Num a)+ => (f (Forward a) -> g (Forward a))+ -> f (a, a)+ -> g (a, a) duF' f = fmap unbundle . f . fmap (uncurry bundle) {-# INLINE duF' #-} @@ -75,7 +85,11 @@ -- -- >>> diff sin 0 -- 1.0-diff :: Num a => (Forward a -> Forward a) -> a -> a+diff+ :: Num a+ => (Forward a -> Forward a)+ -> a+ -> a diff f a = tangent $ apply f a {-# INLINE diff #-} @@ -92,7 +106,11 @@ -- >>> diff' exp 0 -- (1.0,1.0) -diff' :: Num a => (Forward a -> Forward a) -> a -> (a, a)+diff'+ :: Num a+ => (Forward a -> Forward a)+ -> a+ -> (a, a) diff' f a = unbundle $ apply f a {-# INLINE diff' #-} @@ -100,7 +118,11 @@ -- -- >>> diffF (\a -> [sin a, cos a]) 0 -- [1.0,-0.0]-diffF :: (Functor f, Num a) => (Forward a -> f (Forward a)) -> a -> f a+diffF+ :: (Functor f, Num a)+ => (Forward a -> f (Forward a))+ -> a+ -> f a diffF f a = tangent <$> apply f a {-# INLINE diffF #-} @@ -108,50 +130,78 @@ -- -- >>> diffF' (\a -> [sin a, cos a]) 0 -- [(0.0,1.0),(1.0,-0.0)]-diffF' :: (Functor f, Num a) => (Forward a -> f (Forward a)) -> a -> f (a, a)+diffF'+ :: (Functor f, Num a)+ => (Forward a -> f (Forward a))+ -> a+ -> f (a, a) diffF' f a = unbundle <$> apply f a {-# INLINE diffF' #-} -- | A fast, simple, transposed Jacobian computed with forward-mode AD.-jacobianT :: (Traversable f, Functor g, Num a) => (f (Forward a) -> g (Forward a)) -> f a -> f (g a)+jacobianT+ :: (Traversable f, Functor g, Num a)+ => (f (Forward a) -> g (Forward a))+ -> 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) -> (f (Forward a) -> g (Forward a)) -> f a -> f (g b)+jacobianWithT+ :: (Traversable f, Functor g, Num a)+ => (a -> a -> b)+ -> (f (Forward a) -> g (Forward a))+ -> 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) => (f (Forward a) -> g (Forward a)) -> f a -> g (f a)+jacobian+ :: (Traversable f, Traversable g, Num a)+ => (f (Forward a) -> g (Forward a))+ -> 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) -> (f (Forward a) -> g (Forward a)) -> f a -> g (f b)+jacobianWith+ :: (Traversable f, Traversable g, Num a)+ => (a -> a -> b)+ -> (f (Forward a) -> g (Forward a))+ -> f a+ -> g (f b) jacobianWith g f as = transposeWith (const id) t p where (p, t) = bindWith' g' f as g' a ga = g a . tangent <$> ga {-# INLINE jacobianWith #-} -- | Compute the Jacobian using 'Forward' mode 'AD' along with the actual answer.-jacobian' :: (Traversable f, Traversable g, Num a) => (f (Forward a) -> g (Forward a)) -> f a -> g (a, f a)+jacobian'+ :: (Traversable f, Traversable g, Num a)+ => (f (Forward a) -> g (Forward a))+ -> f a+ -> g (a, f a) jacobian' f as = transposeWith row t p where (p, t) = bind' f as row x as' = (primal x, tangent <$> as') {-# INLINE jacobian' #-} -- | 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) -> (f (Forward a) -> g (Forward a)) -> f a -> g (a, f b)+jacobianWith'+ :: (Traversable f, Traversable g, Num a)+ => (a -> a -> b)+ -> (f (Forward a) -> g (Forward 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')@@ -160,44 +210,69 @@ -- | 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) => (f (Forward a) -> Forward a) -> f a -> f a+-- 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.-grad' :: (Traversable f, Num a) => (f (Forward a) -> Forward a) -> f a -> (a, f a)+-- 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 {-# 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) -> (f (Forward a) -> Forward a) -> f a -> f b+-- 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 #-} -- | 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)])-gradWith' :: (Traversable f, Num a) => (a -> a -> b) -> (f (Forward a) -> Forward a) -> f a -> (a, f b)+gradWith'+ :: (Traversable f, Num a)+ => (a -> a -> b)+ -> (f (Forward a) -> Forward a)+ -> 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) => (f (On (Forward (Forward a))) -> On (Forward (Forward a))) -> f (a, a) -> f a+hessianProduct+ :: (Traversable f, Num a)+ => (f (On (Forward (Forward a))) -> On (Forward (Forward a)))+ -> 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) => (f (On (Forward (Forward a))) -> On (Forward (Forward a))) -> f (a, a) -> f (a, a)+hessianProduct'+ :: (Traversable f, Num a)+ => (f (On (Forward (Forward a))) -> On (Forward (Forward a)))+ -> f (a, a) -> f (a, a) hessianProduct' f = duF' $ grad $ off . f . fmap On {-# INLINE hessianProduct' #-}
src/Numeric/AD/Rank1/Forward/Double.hs view
@@ -1,4 +1,13 @@-{-# LANGUAGE CPP #-}+-----------------------------------------------------------------------------+-- |+-- Copyright : (c) Edward Kmett 2010-2021+-- License : BSD3+-- Maintainer : ekmett@gmail.com+-- Stability : experimental+-- Portability : GHC only+--+-----------------------------------------------------------------------------+ module Numeric.AD.Rank1.Forward.Double ( ForwardDouble -- * Gradient@@ -26,30 +35,42 @@ , duF' ) where -#if __GLASGOW_HASKELL__ < 710-import Control.Applicative-import Data.Traversable (Traversable)-#endif 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 => (f ForwardDouble -> ForwardDouble) -> f (Double, Double) -> Double+du+ :: Functor f+ => (f ForwardDouble -> ForwardDouble)+ -> 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 => (f ForwardDouble -> ForwardDouble) -> f (Double, Double) -> (Double, Double)+du'+ :: Functor f+ => (f ForwardDouble -> ForwardDouble)+ -> 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) => (f ForwardDouble -> g ForwardDouble) -> f (Double, Double) -> g Double+duF+ :: (Functor f, Functor g)+ => (f ForwardDouble -> g ForwardDouble)+ -> 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) => (f ForwardDouble -> g ForwardDouble) -> f (Double, Double) -> g (Double, Double)+duF'+ :: (Functor f, Functor g)+ => (f ForwardDouble -> g ForwardDouble)+ -> f (Double, Double)+ -> g (Double, Double) duF' f = fmap unbundle . f . fmap (uncurry bundle) {-# INLINE duF' #-} @@ -57,7 +78,10 @@ -- -- >>> diff sin 0 -- 1.0-diff :: (ForwardDouble -> ForwardDouble) -> Double -> Double+diff+ :: (ForwardDouble -> ForwardDouble)+ -> Double+ -> Double diff f a = tangent $ apply f a {-# INLINE diff #-} @@ -73,7 +97,10 @@ -- -- >>> diff' exp 0 -- (1.0,1.0)-diff' :: (ForwardDouble -> ForwardDouble) -> Double -> (Double, Double)+diff'+ :: (ForwardDouble -> ForwardDouble)+ -> Double+ -> (Double, Double) diff' f a = unbundle $ apply f a {-# INLINE diff' #-} @@ -81,7 +108,11 @@ -- -- >>> diffF (\a -> [sin a, cos a]) 0 -- [1.0,-0.0]-diffF :: Functor f => (ForwardDouble -> f ForwardDouble) -> Double -> f Double+diffF+ :: Functor f+ => (ForwardDouble -> f ForwardDouble)+ -> Double+ -> f Double diffF f a = tangent <$> apply f a {-# INLINE diffF #-} @@ -89,48 +120,78 @@ -- -- >>> diffF' (\a -> [sin a, cos a]) 0 -- [(0.0,1.0),(1.0,-0.0)]-diffF' :: Functor f => (ForwardDouble -> f ForwardDouble) -> Double -> f (Double, Double)+diffF'+ :: Functor f+ => (ForwardDouble -> f ForwardDouble)+ -> 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) => (f ForwardDouble -> g ForwardDouble) -> f Double -> f (g Double)+jacobianT+ :: (Traversable f, Functor g)+ => (f ForwardDouble -> g ForwardDouble)+ -> 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) -> (f ForwardDouble -> g ForwardDouble) -> f Double -> f (g b)+jacobianWithT+ :: (Traversable f, Functor g)+ => (Double -> Double -> b)+ -> (f ForwardDouble -> g ForwardDouble)+ -> 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) => (f ForwardDouble -> g ForwardDouble) -> f Double -> g (f Double)+jacobian+ :: (Traversable f, Traversable g)+ => (f ForwardDouble -> g ForwardDouble)+ -> 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) -> (f ForwardDouble -> g ForwardDouble) -> f Double -> g (f b)+jacobianWith+ :: (Traversable f, Traversable g)+ => (Double -> Double -> b)+ -> (f ForwardDouble -> g ForwardDouble)+ -> 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) => (f ForwardDouble -> g ForwardDouble) -> f Double -> g (Double, f Double)+jacobian'+ :: (Traversable f, Traversable g)+ => (f ForwardDouble -> g ForwardDouble)+ -> 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) -> (f ForwardDouble -> g ForwardDouble) -> f Double -> g (Double, f b)+jacobianWith'+ :: (Traversable f, Traversable g)+ => (Double -> Double -> b)+ -> (f ForwardDouble -> g ForwardDouble)+ -> 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')@@ -139,15 +200,23 @@ -- | 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 => (f ForwardDouble -> ForwardDouble) -> f Double -> f Double+-- 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.-grad' :: Traversable f => (f ForwardDouble -> ForwardDouble) -> f Double -> (Double, f Double)+-- 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 (b, bs) = bind' f as@@ -155,18 +224,28 @@ -- | 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) -> (f ForwardDouble -> ForwardDouble) -> f Double -> f b+-- 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 #-} -- | 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)])-gradWith' :: Traversable f => (Double -> Double -> b) -> (f ForwardDouble -> ForwardDouble) -> f Double -> (Double, f b)+gradWith'+ :: Traversable f+ => (Double -> Double -> b)+ -> (f ForwardDouble -> ForwardDouble)+ -> f Double+ -> (Double, f b) gradWith' g f as = (primal $ f (auto <$> as), bindWith g (tangent . f) as) {-# INLINE gradWith' #-}
src/Numeric/AD/Rank1/Halley.hs view
@@ -1,7 +1,7 @@ {-# LANGUAGE CPP #-} ----------------------------------------------------------------------------- -- |--- Copyright : (c) Edward Kmett 2010-2015+-- Copyright : (c) Edward Kmett 2010-2021 -- License : BSD3 -- Maintainer : ekmett@gmail.com -- Stability : experimental@@ -61,7 +61,9 @@ findZeroNoEq :: Fractional a => (Tower a -> Tower a) -> a -> [a] findZeroNoEq f = iterate go where go x = xn where- (y:y':y'':_) = diffs0 f x+ (y,y',y'') = case diffs0 f x of+ (z:z':z'':_) -> (z,z',z'')+ _ -> error "findZeroNoEq: Impossible (diffs0 should produce an infinite list)" xn = x - 2*y*y'/(2*y'*y'-y*y'') -- 9.606671960457536 bits error -- = x - recip (y'/y - y''/ y') -- "improved error" = 6.640625e-2 bits -- = x - y' / (y'/y/y' - y''/2) -- "improved error" = 1.4
+ src/Numeric/AD/Rank1/Halley/Double.hs view
@@ -0,0 +1,122 @@+{-# LANGUAGE CPP #-}+-----------------------------------------------------------------------------+-- |+-- Copyright : (c) Edward Kmett 2010-2021+-- License : BSD3+-- Maintainer : ekmett@gmail.com+-- Stability : experimental+-- Portability : GHC only+--+-- Root finding using Halley's rational method (the second in+-- the class of Householder methods). Assumes the function is three+-- times continuously differentiable and converges cubically when+-- progress can be made.+--+-----------------------------------------------------------------------------++module Numeric.AD.Rank1.Halley.Double+ (+ -- * Halley's Method (Tower AD)+ findZero+ , inverse+ , fixedPoint+ , extremum+ ) where++import Prelude hiding (all)+import Numeric.AD.Internal.Forward (Forward)+import Numeric.AD.Internal.On+import Numeric.AD.Internal.Tower.Double (TowerDouble)+import Numeric.AD.Mode+import Numeric.AD.Rank1.Tower.Double (diffs0)+import Numeric.AD.Rank1.Forward (diff)+import Numeric.AD.Internal.Combinators (takeWhileDifferent)++-- $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+-- ("it converges"), no further elements are returned.+--+-- Examples:+--+-- >>> take 10 $ findZero (\x->x^2-4) 1+-- [1.0,1.8571428571428572,1.9997967892704736,1.9999999999994755,2.0]+findZero :: (TowerDouble -> TowerDouble) -> Double -> [Double]+findZero f = takeWhileDifferent . findZeroNoEq f+{-# INLINE findZero #-}++-- | The 'findZeroNoEq' function behaves the same as 'findZero' except that it+-- doesn't truncate the list once the results become constant. This means it+-- can be used with types without an 'Eq' instance.+findZeroNoEq :: (TowerDouble -> TowerDouble) -> Double -> [Double]+findZeroNoEq f = iterate go where+ go x = xn where+ (y,y',y'') = case diffs0 f x of+ (z:z':z'':_) -> (z,z',z'')+ _ -> error "findZeroNoEq: Impossible (diffs0 should produce an infinite list)"+ xn = x - 2*y*y'/(2*y'*y'-y*y'') -- 9.606671960457536 bits error+ -- = x - recip (y'/y - y''/ y') -- "improved error" = 6.640625e-2 bits+ -- = x - y' / (y'/y/y' - y''/2) -- "improved error" = 1.4+#ifdef HERBIE+{-# ANN findZeroNoEq "NoHerbie" #-}+#endif+{-# INLINE findZeroNoEq #-}++-- | The 'inverse' function inverts 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+-- ("it converges"), no further elements are returned.+--+-- 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 :: (TowerDouble -> TowerDouble) -> Double -> Double -> [Double]+inverse f x0 = takeWhileDifferent . inverseNoEq f x0+{-# INLINE inverse #-}++-- | The 'inverseNoEq' function behaves the same as 'inverse' except that it+-- doesn't truncate the list once the results become constant. This means it+-- can be used with types without an 'Eq' instance.+inverseNoEq :: (TowerDouble -> TowerDouble) -> Double -> Double -> [Double]+inverseNoEq f x0 y = findZeroNoEq (\x -> f x - auto y) x0+{-# INLINE inverseNoEq #-}++-- | The 'fixedPoint' function find a fixedpoint of a scalar+-- function using Halley's method; its output is a stream of+-- increasingly accurate results. (Modulo the usual caveats.)+--+-- If the stream becomes constant ("it converges"), no further+-- elements are returned.+--+-- >>> last $ take 10 $ fixedPoint cos 1+-- 0.7390851332151607+fixedPoint :: (TowerDouble -> TowerDouble) -> Double -> [Double]+fixedPoint f = takeWhileDifferent . fixedPointNoEq f+{-# INLINE fixedPoint #-}++-- | The 'fixedPointNoEq' function behaves the same as 'fixedPoint' except that+-- it doesn't truncate the list once the results become constant. This means it+-- can be used with types without an 'Eq' instance.+fixedPointNoEq :: (TowerDouble -> TowerDouble) -> Double -> [Double]+fixedPointNoEq f = findZeroNoEq (\x -> f x - x)+{-# INLINE fixedPointNoEq #-}++-- | 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+-- constant ("it converges"), no further elements are returned.+--+-- >>> take 10 $ extremum cos 1+-- [1.0,0.29616942658570555,4.59979519460002e-3,1.6220740159042513e-8,0.0]+extremum :: (On (Forward TowerDouble) -> On (Forward TowerDouble)) -> Double -> [Double]+extremum f = takeWhileDifferent . extremumNoEq f+{-# INLINE extremum #-}++-- | The 'extremumNoEq' function behaves the same as 'extremum' except that it+-- doesn't truncate the list once the results become constant. This means it+-- can be used with types without an 'Eq' instance.+extremumNoEq :: (On (Forward TowerDouble) -> On (Forward TowerDouble)) -> Double -> [Double]+extremumNoEq f = findZeroNoEq (diff (off . f . On))+{-# INLINE extremumNoEq #-}
src/Numeric/AD/Rank1/Kahn.hs view
@@ -1,223 +1,20 @@ {-# LANGUAGE CPP #-}-{-# LANGUAGE FlexibleInstances #-}-{-# LANGUAGE FunctionalDependencies #-}-{-# LANGUAGE MultiParamTypeClasses #-}-{-# LANGUAGE Rank2Types #-}-{-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE TemplateHaskell #-}-{-# LANGUAGE UndecidableInstances #-}--------------------------------------------------------------------------------- |--- Copyright : (c) Edward Kmett 2010-2015--- License : BSD3--- Maintainer : ekmett@gmail.com--- Stability : experimental--- Portability : GHC only------ This module provides reverse-mode Automatic Differentiation using post-hoc linear time--- topological sorting.------ For reverse mode AD we use 'System.Mem.StableName.StableName' to recover sharing information from--- the tape to avoid combinatorial explosion, and thus run asymptotically faster--- than it could without such sharing information, but the use of side-effects--- contained herein is benign.--------------------------------------------------------------------------------- +#define MODULE \ module Numeric.AD.Rank1.Kahn- ( Kahn- , auto- -- * 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 -#if __GLASGOW_HASKELL__ < 710-import Control.Applicative ((<$>))-import Data.Traversable (Traversable)-#endif-import Data.Functor.Compose-import Numeric.AD.Internal.On-import Numeric.AD.Internal.Kahn-import Numeric.AD.Mode---- $setup------ >>> import Numeric.AD.Internal.Doctest---- | 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) => (f (Kahn a) -> Kahn a) -> f a -> f a-grad f as = unbind vs (partialArray bds $ f vs) where- (vs,bds) = bind as-{-# INLINE grad #-}---- | The 'grad'' function calculates the result and gradient of a non-scalar-to-scalar function with 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) => (f (Kahn a) -> Kahn a) -> f a -> (a, f a)-grad' f as = (primal r, unbind vs $ partialArray bds r) where- (vs, bds) = bind as- r = f vs-{-# INLINE grad' #-}---- | @'grad' g f@ function calculates the gradient of a non-scalar-to-scalar function @f@ with kahn-mode AD in a single pass.--- The gradient is combined element-wise with the argument using the function @g@.------ @--- 'grad' = 'gradWith' (\_ dx -> dx)--- 'id' = 'gradWith' const--- @-------gradWith :: (Traversable f, Num a) => (a -> a -> b) -> (f (Kahn a) -> Kahn a) -> f a -> f b-gradWith g f as = unbindWith g vs (partialArray bds $ f vs) where- (vs,bds) = bind as-{-# INLINE gradWith #-}---- | @'grad'' g f@ calculates the result and gradient of a non-scalar-to-scalar function @f@ with 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) -> (f (Kahn a) -> Kahn a) -> f a -> (a, f b)-gradWith' g f as = (primal r, unbindWith g vs $ partialArray bds r) where- (vs, bds) = bind as- r = f vs-{-# INLINE gradWith' #-}---- | The 'jacobian' function calculates the jacobian of a non-scalar-to-non-scalar function with 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) => (f (Kahn a) -> g (Kahn a)) -> f a -> g (f a)-jacobian f as = unbind vs . partialArray bds <$> f vs where- (vs, bds) = bind as-{-# INLINE jacobian #-}---- | The 'jacobian'' function calculates both the result and the Jacobian of a nonscalar-to-nonscalar function, using @m@ invocations of 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) => (f (Kahn a) -> g (Kahn a)) -> f a -> g (a, f a)-jacobian' f as = row <$> f vs where- (vs, bds) = bind as- row a = (primal a, unbind vs (partialArray bds a))-{-# INLINE jacobian' #-}---- | 'jacobianWith g f' calculates the Jacobian of a non-scalar-to-non-scalar function @f@ with 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@.------ @--- 'jacobian' = 'jacobianWith' (\_ dx -> dx)--- 'jacobianWith' 'const' = (\f x -> 'const' x '<$>' f x)--- @-jacobianWith :: (Traversable f, Functor g, Num a) => (a -> a -> b) -> (f (Kahn a) -> g (Kahn a)) -> f a -> g (f b)-jacobianWith g f as = unbindWith g vs . partialArray bds <$> f vs where- (vs, bds) = bind as-{-# INLINE jacobianWith #-}---- | 'jacobianWith' g f' calculates both the result and the Jacobian of a nonscalar-to-nonscalar function @f@, using @m@ invocations of 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) -> (f (Kahn a) -> g (Kahn a)) -> f a -> g (a, f b)-jacobianWith' g f as = row <$> f vs where- (vs, bds) = bind as- row a = (primal a, unbindWith g vs (partialArray bds a))-{-# INLINE jacobianWith' #-}---- | Compute the derivative of a function.------ >>> diff sin 0--- 1.0------ >>> cos 0--- 1.0-diff :: Num a => (Kahn a -> Kahn a) -> a -> a-diff f a = derivative $ f (var a 0)-{-# INLINE diff #-}---- | The 'diff'' function calculates the value and derivative, as a--- pair, of a scalar-to-scalar function.--------- >>> diff' sin 0--- (0.0,1.0)-diff' :: Num a => (Kahn a -> Kahn a) -> a -> (a, a)-diff' f a = derivative' $ f (var a 0)-{-# INLINE diff' #-}---- | Compute the derivatives of a function that returns a vector with regards to its single input.------ >>> diffF (\a -> [sin a, cos a]) 0--- [1.0,0.0]-diffF :: (Functor f, Num a) => (Kahn a -> f (Kahn a)) -> a -> f a-diffF f a = derivative <$> f (var a 0)-{-# INLINE diffF #-}---- | Compute the derivatives of a function that returns a vector with regards to its single input--- as well as the primal answer.------ >>> diffF' (\a -> [sin a, cos a]) 0--- [(0.0,1.0),(1.0,0.0)]-diffF' :: (Functor f, Num a) => (Kahn a -> f (Kahn a)) -> a -> f (a, a)-diffF' f a = derivative' <$> f (var a 0)-{-# INLINE diffF' #-}-+#define AD_EXPORT Kahn --- | 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) => (f (On (Kahn (Kahn a))) -> (On (Kahn (Kahn a)))) -> f a -> f (f a)-hessian f = jacobian (grad (off . f . fmap On))+#define IMPORTS \+import Numeric.AD.Internal.Kahn --- | 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 :: RDouble]--- [[[0.0,1.0],[1.0,0.0]],[[0.0,0.0],[0.0,0.0]],[[-1.131204383757,-2.471726672005],[-2.471726672005,1.131204383757]]]-hessianF :: (Traversable f, Functor g, Num a) => (f (On (Kahn (Kahn a))) -> g (On (Kahn (Kahn a)))) -> f a -> g (f (f a))-hessianF f = getCompose . jacobian (Compose . jacobian (fmap off . f . fmap On))+#define UNBINDWITH unbindWith+#define JACOBIAN jacobian+#define GRAD grad --- $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!+#define AD_TYPE (Kahn a)+#define SCALAR_TYPE a+#define BASE0_1(x) x =>+#define BASE1_1(x,y) (x,y)+#define BASE2_1(x,y,z) (x,y,z)+#include "rank1_kahn.h"
+ src/Numeric/AD/Rank1/Kahn/Double.hs view
@@ -0,0 +1,22 @@+{-# LANGUAGE CPP #-}++#define MODULE \+module Numeric.AD.Rank1.Kahn.Double++#define AD_EXPORT KahnDouble++#define IMPORTS \+import Numeric.AD.Internal.Kahn (Kahn); \+import qualified Numeric.AD.Rank1.Kahn as Kahn; \+import Numeric.AD.Internal.Kahn.Double++#define UNBINDWITH unbindWithUArray+#define GRAD Kahn.grad+#define JACOBIAN Kahn.jacobian++#define AD_TYPE KahnDouble+#define SCALAR_TYPE Double+#define BASE0_1(x)+#define BASE1_1(x,y) x+#define BASE2_1(x,y,z) (x,y)+#include "rank1_kahn.h"
+ src/Numeric/AD/Rank1/Kahn/Float.hs view
@@ -0,0 +1,22 @@+{-# LANGUAGE CPP #-}++#define MODULE \+module Numeric.AD.Rank1.Kahn.Float++#define AD_EXPORT KahnFloat++#define IMPORTS \+import Numeric.AD.Internal.Kahn (Kahn); \+import qualified Numeric.AD.Rank1.Kahn as Kahn; \+import Numeric.AD.Internal.Kahn.Float++#define UNBINDWITH unbindWithUArray+#define GRAD Kahn.grad+#define JACOBIAN Kahn.jacobian++#define AD_TYPE KahnFloat+#define SCALAR_TYPE Float+#define BASE0_1(x)+#define BASE1_1(x,y) x+#define BASE2_1(x,y,z) (x,y)+#include "rank1_kahn.h"
src/Numeric/AD/Rank1/Newton.hs view
@@ -1,11 +1,10 @@-{-# LANGUAGE CPP #-} {-# LANGUAGE BangPatterns #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TypeFamilies #-} ----------------------------------------------------------------------------- -- |--- Copyright : (c) Edward Kmett 2010-2015+-- Copyright : (c) Edward Kmett 2010-2021 -- License : BSD3 -- Maintainer : ekmett@gmail.com -- Stability : experimental@@ -31,9 +30,6 @@ import Prelude hiding (all, mapM) import Data.Foldable (all)-#if __GLASGOW_HASKELL__ < 710-import Data.Traversable-#endif import Numeric.AD.Mode import Numeric.AD.Rank1.Forward (Forward, diff, diff') import Numeric.AD.Rank1.Kahn as Kahn (Kahn, gradWith')
src/Numeric/AD/Rank1/Newton/Double.hs view
@@ -1,11 +1,9 @@-{-# LANGUAGE CPP #-}-{-# LANGUAGE BangPatterns #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TypeFamilies #-} ----------------------------------------------------------------------------- -- |--- Copyright : (c) Edward Kmett 2015+-- Copyright : (c) Edward Kmett 2015-2021 -- License : BSD3 -- Maintainer : ekmett@gmail.com -- Stability : experimental@@ -17,13 +15,9 @@ ( -- * Newton's Method (Forward) findZero- , findZeroNoEq , inverse- , inverseNoEq , fixedPoint- , fixedPointNoEq , extremum- , extremumNoEq ) where import Prelude hiding (all, mapM)@@ -50,7 +44,7 @@ -- | The 'findZeroNoEq' function behaves the same as 'findZero' except that it -- doesn't truncate the list once the results become constant. findZeroNoEq :: (ForwardDouble -> ForwardDouble) -> Double -> [Double]-findZeroNoEq f = iterate go where+findZeroNoEq f = takeWhileDifferent . iterate go where go x = xn where (y,y') = diff' f x xn = x - y/y'
src/Numeric/AD/Rank1/Sparse.hs view
@@ -1,7 +1,6 @@-{-# LANGUAGE CPP #-} ----------------------------------------------------------------------------- -- |--- Copyright : (c) Edward Kmett 2010-2015+-- Copyright : (c) Edward Kmett 2010-2021 -- License : BSD3 -- Maintainer : ekmett@gmail.com -- Stability : experimental@@ -46,9 +45,6 @@ ) where import Control.Comonad-#if __GLASGOW_HASKELL__ < 710-import Data.Traversable-#endif import Control.Comonad.Cofree import Numeric.AD.Jet import Numeric.AD.Internal.Sparse
+ src/Numeric/AD/Rank1/Sparse/Double.hs view
@@ -0,0 +1,192 @@+-----------------------------------------------------------------------------+-- |+-- Copyright : (c) Edward Kmett 2010-2021+-- License : BSD3+-- Maintainer : ekmett@gmail.com+-- Stability : experimental+-- Portability : GHC only+--+-- Higher order derivatives via a \"dual number tower\".+--+-----------------------------------------------------------------------------++module Numeric.AD.Rank1.Sparse.Double+ ( SparseDouble+ , auto+ -- * Sparse Gradients+ , grad+ , grad'+ , gradWith+ , gradWith'+ -- * Variadic Gradients+ -- $vgrad+ , Grad+ , vgrad+ -- * Higher-Order Gradients+ , grads+ -- * Variadic Higher-Order Gradients+ , Grads+ , vgrads++ -- * Sparse Jacobians (synonyms)+ , jacobian+ , jacobian'+ , jacobianWith+ , jacobianWith'+ , jacobians++ -- * Sparse Hessians+ , hessian+ , hessian'++ , hessianF+ , hessianF'++ ) where++import Control.Comonad+import Control.Comonad.Cofree+import Numeric.AD.Jet+import Numeric.AD.Internal.Sparse.Double+import Numeric.AD.Internal.Combinators+import Numeric.AD.Mode++second :: (a -> b) -> (c, a) -> (c, b)+second g (a,b) = (a, g b)+{-# INLINE second #-}++grad+ :: Traversable f+ => (f SparseDouble -> SparseDouble)+ -> f Double -> f Double+grad f as = d as $ apply f as+{-# INLINE grad #-}++grad'+ :: Traversable f+ => (f SparseDouble -> SparseDouble)+ -> f Double -> (Double, f Double)+grad' f as = d' as $ apply f as+{-# INLINE grad' #-}++gradWith+ :: Traversable f+ => (Double -> Double -> b)+ -> (f SparseDouble -> SparseDouble)+ -> f Double+ -> f b+gradWith g f as = zipWithT g as $ grad f as+{-# INLINE gradWith #-}++gradWith'+ :: Traversable f+ => (Double -> Double -> b)+ -> (f SparseDouble -> SparseDouble)+ -> f Double+ -> (Double, f b)+gradWith' g f as = second (zipWithT g as) $ grad' f as+{-# INLINE gradWith' #-}++jacobian+ :: (Traversable f, Functor g)+ => (f SparseDouble -> g SparseDouble)+ -> f Double+ -> g (f Double)+jacobian f as = d as <$> apply f as+{-# INLINE jacobian #-}++jacobian'+ :: (Traversable f, Functor g)+ => (f SparseDouble -> g SparseDouble)+ -> f Double+ -> g (Double, f Double)+jacobian' f as = d' as <$> apply f as+{-# INLINE jacobian' #-}++jacobianWith+ :: (Traversable f, Functor g)+ => (Double -> Double -> b)+ -> (f SparseDouble -> g SparseDouble)+ -> f Double+ -> g (f b)+jacobianWith g f as = zipWithT g as <$> jacobian f as+{-# INLINE jacobianWith #-}++jacobianWith'+ :: (Traversable f, Functor g)+ => (Double -> Double -> b)+ -> (f SparseDouble -> g SparseDouble)+ -> f Double+ -> g (Double, f b)+jacobianWith' g f as = second (zipWithT g as) <$> jacobian' f as+{-# INLINE jacobianWith' #-}++grads+ :: Traversable f+ => (f SparseDouble -> SparseDouble)+ -> f Double+ -> Cofree f Double+grads f as = ds as $ apply f as+{-# INLINE grads #-}++jacobians+ :: (Traversable f, Functor g)+ => (f SparseDouble -> g SparseDouble)+ -> f Double+ -> g (Cofree f Double)+jacobians f as = ds as <$> apply f as+{-# INLINE jacobians #-}++d2 :: Functor f+ => Cofree f a+ -> f (f a)+d2 = headJet . tailJet . tailJet . jet+{-# INLINE d2 #-}++d2'+ :: Functor f+ => Cofree f a+ -> (a, f (a, f a))+d2' (a :< as) = (a, fmap (\(da :< das) -> (da, extract <$> das)) as)+{-# INLINE d2' #-}++hessian+ :: Traversable f+ => (f SparseDouble -> SparseDouble)+ -> f Double+ -> f (f Double)+hessian f as = d2 $ grads f as+{-# INLINE hessian #-}++hessian'+ :: Traversable f+ => (f SparseDouble -> SparseDouble)+ -> f Double+ -> (Double, f (Double, f Double))+hessian' f as = d2' $ grads f as+{-# INLINE hessian' #-}++hessianF+ :: (Traversable f, Functor g)+ => (f SparseDouble -> g SparseDouble)+ -> f Double+ -> g (f (f Double))+hessianF f as = d2 <$> jacobians f as+{-# INLINE hessianF #-}++hessianF'+ :: (Traversable f, Functor g)+ => (f SparseDouble -> g SparseDouble)+ -> f Double+ -> g (Double, f (Double, f Double))+hessianF' f as = d2' <$> jacobians f as+{-# INLINE hessianF' #-}++-- $vgrad+--+-- Variadic combinators for variadic mixed-mode automatic differentiation.+--+-- Unfortunately, variadicity comes at the expense of being able to use+-- quantification to avoid sensitivity confusion, so be careful when+-- counting the number of 'auto' calls you use when taking the gradient+-- of a function that takes gradients!
src/Numeric/AD/Rank1/Tower.hs view
@@ -1,9 +1,8 @@-{-# LANGUAGE CPP #-} {-# LANGUAGE Rank2Types #-} {-# LANGUAGE BangPatterns #-} ----------------------------------------------------------------------------- -- |--- Copyright : (c) Edward Kmett 2010-2015+-- Copyright : (c) Edward Kmett 2010-2021 -- License : BSD3 -- Maintainer : ekmett@gmail.com -- Stability : experimental@@ -40,99 +39,169 @@ , dus0F -- answer and all zero padded directional derivatives of (f a -> g a) ) where -#if __GLASGOW_HASKELL__ < 710-import Data.Functor ((<$>))-#endif import Numeric.AD.Internal.Tower import Numeric.AD.Mode -- | Compute the answer and all derivatives of a function @(a -> a)@-diffs :: Num a => (Tower a -> Tower a) -> a -> [a]+diffs+ :: Num a+ => (Tower a -> Tower a)+ -> a+ -> [a] diffs f a = getADTower $ apply f a {-# INLINE diffs #-} -- | Compute the zero-padded derivatives of a function @(a -> a)@-diffs0 :: Num a => (Tower a -> Tower a) -> a -> [a]+diffs0+ :: Num a+ => (Tower a -> Tower a)+ -> a+ -> [a] diffs0 f a = zeroPad (diffs f a) {-# INLINE diffs0 #-} -- | Compute the answer and all derivatives of a function @(a -> f a)@-diffsF :: (Functor f, Num a) => (Tower a -> f (Tower a)) -> a -> f [a]+diffsF+ :: (Functor f, Num a)+ => (Tower a -> f (Tower a))+ -> a+ -> f [a] diffsF f a = getADTower <$> apply f a {-# INLINE diffsF #-} -- | Compute the zero-padded derivatives of a function @(a -> f a)@-diffs0F :: (Functor f, Num a) => (Tower a -> f (Tower a)) -> a -> f [a]-diffs0F f a = (zeroPad . getADTower) <$> apply f a+diffs0F+ :: (Functor f, Num a)+ => (Tower a -> f (Tower a))+ -> a+ -> f [a]+diffs0F f a = zeroPad . getADTower <$> apply f a {-# INLINE diffs0F #-} -- | @taylor f x@ compute the Taylor series of @f@ around @x@.-taylor :: Fractional a => (Tower a -> Tower a) -> a -> a -> [a]+taylor+ :: Fractional a+ => (Tower a -> Tower a)+ -> a+ -> a+ -> [a] taylor f x dx = go 1 1 (diffs f x) where go !n !acc (a:as) = a * acc : go (n + 1) (acc * dx / n) as go _ _ [] = [] -- | @taylor0 f x@ compute the Taylor series of @f@ around @x@, zero-padded.-taylor0 :: Fractional a => (Tower a -> Tower a) -> a -> a -> [a]+taylor0+ :: Fractional a+ => (Tower a -> Tower a)+ -> a+ -> a+ -> [a] taylor0 f x dx = zeroPad (taylor f x dx) {-# INLINE taylor0 #-} -- | @maclaurin f@ compute the Maclaurin series of @f@-maclaurin :: Fractional a => (Tower a -> Tower a) -> a -> [a]+maclaurin+ :: Fractional a+ => (Tower a -> Tower a)+ -> a+ -> [a] maclaurin f = taylor f 0 {-# INLINE maclaurin #-} -- | @maclaurin f@ compute the Maclaurin series of @f@, zero-padded-maclaurin0 :: Fractional a => (Tower a -> Tower a) -> a -> [a]+maclaurin0+ :: Fractional a+ => (Tower a -> Tower a)+ -> a+ -> [a] maclaurin0 f = taylor0 f 0 {-# INLINE maclaurin0 #-} -- | Compute the first derivative of a function @(a -> a)@-diff :: Num a => (Tower a -> Tower a) -> a -> a+diff+ :: Num a+ => (Tower a -> Tower a)+ -> a+ -> a diff f = d . diffs f {-# INLINE diff #-} -- | Compute the answer and first derivative of a function @(a -> a)@-diff' :: Num a => (Tower a -> Tower a) -> a -> (a, a)+diff'+ :: Num a+ => (Tower a -> Tower a)+ -> a+ -> (a, a) diff' f = d' . diffs f {-# INLINE diff' #-} -- | Compute a directional derivative of a function @(f a -> a)@-du :: (Functor f, Num a) => (f (Tower a) -> Tower a) -> f (a, a) -> a+du+ :: (Functor f, Num a)+ => (f (Tower a) -> Tower a)+ -> f (a, a) -> a du f = d . getADTower . f . fmap withD {-# INLINE du #-} -- | Compute the answer and a directional derivative of a function @(f a -> a)@-du' :: (Functor f, Num a) => (f (Tower a) -> Tower a) -> f (a, a) -> (a, a)+du'+ :: (Functor f, Num a)+ => (f (Tower a) -> Tower a)+ -> f (a, a)+ -> (a, a) du' f = d' . getADTower . f . fmap withD {-# INLINE du' #-} -- | Compute a directional derivative of a function @(f a -> g a)@-duF :: (Functor f, Functor g, Num a) => (f (Tower a) -> g (Tower a)) -> f (a, a) -> g a+duF+ :: (Functor f, Functor g, Num a)+ => (f (Tower a) -> g (Tower a))+ -> f (a, a)+ -> g a duF f = fmap (d . getADTower) . f . fmap withD {-# INLINE duF #-} -- | Compute the answer and a directional derivative of a function @(f a -> g a)@-duF' :: (Functor f, Functor g, Num a) => (f (Tower a) -> g (Tower a)) -> f (a, a) -> g (a, a)+duF'+ :: (Functor f, Functor g, Num a)+ => (f (Tower a) -> g (Tower a))+ -> f (a, a)+ -> g (a, a) duF' f = fmap (d' . getADTower) . f . fmap withD {-# INLINE duF' #-} -- | Given a function @(f a -> a)@, and a tower of derivatives, compute the corresponding directional derivatives.-dus :: (Functor f, Num a) => (f (Tower a) -> Tower a) -> f [a] -> [a]+dus+ :: (Functor f, Num a)+ => (f (Tower a) -> Tower a)+ -> f [a]+ -> [a] dus f = getADTower . f . fmap tower {-# INLINE dus #-} -- | Given a function @(f a -> a)@, and a tower of derivatives, compute the corresponding directional derivatives, zero-padded-dus0 :: (Functor f, Num a) => (f (Tower a) -> Tower a) -> f [a] -> [a]+dus0+ :: (Functor f, Num a)+ => (f (Tower a) -> Tower a)+ -> f [a]+ -> [a] dus0 f = zeroPad . getADTower . f . fmap tower {-# INLINE dus0 #-} -- | Given a function @(f a -> g a)@, and a tower of derivatives, compute the corresponding directional derivatives-dusF :: (Functor f, Functor g, Num a) => (f (Tower a) -> g (Tower a)) -> f [a] -> g [a]+dusF+ :: (Functor f, Functor g, Num a)+ => (f (Tower a) -> g (Tower a))+ -> f [a]+ -> g [a] dusF f = fmap getADTower . f . fmap tower {-# INLINE dusF #-} -- | Given a function @(f a -> g a)@, and a tower of derivatives, compute the corresponding directional derivatives, zero-padded-dus0F :: (Functor f, Functor g, Num a) => (f (Tower a) -> g (Tower a)) -> f [a] -> g [a]+dus0F+ :: (Functor f, Functor g, Num a)+ => (f (Tower a) -> g (Tower a))+ -> f [a]+ -> g [a] dus0F f = fmap getADTower . f . fmap tower {-# INLINE dus0F #-}
+ src/Numeric/AD/Rank1/Tower/Double.hs view
@@ -0,0 +1,199 @@+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE BangPatterns #-}+-----------------------------------------------------------------------------+-- |+-- Copyright : (c) Edward Kmett 2010-2021+-- License : BSD3+-- Maintainer : ekmett@gmail.com+-- Stability : experimental+-- Portability : GHC only+--+-- Higher order derivatives via a \"dual number tower\".+--+-----------------------------------------------------------------------------++module Numeric.AD.Rank1.Tower.Double+ ( TowerDouble+ , auto+ -- * 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 (f a -> a)+ , du' -- answer and directional derivative of (f a -> a)+ , dus -- answer and all directional derivatives of (f a -> a)+ , dus0 -- answer and all zero padded directional derivatives of (f a -> a)+ , duF -- directional derivative of (f a -> g a)+ , duF' -- answer and directional derivative of (f a -> g a)+ , dusF -- answer and all directional derivatives of (f a -> g a)+ , dus0F -- answer and all zero padded directional derivatives of (f a -> g a)+ ) where++import Numeric.AD.Internal.Tower.Double+import Numeric.AD.Mode++-- | Compute the answer and all derivatives of a function @(a -> a)@+diffs+ :: (TowerDouble -> TowerDouble)+ -> Double+ -> [Double]+diffs f a = getADTower $ apply f a+{-# INLINE diffs #-}++-- | Compute the zero-padded derivatives of a function @(a -> a)@+diffs0+ :: (TowerDouble -> TowerDouble)+ -> Double+ -> [Double]+diffs0 f a = zeroPad (diffs f a)+{-# INLINE diffs0 #-}++-- | Compute the answer and all derivatives of a function @(a -> f a)@+diffsF+ :: Functor f+ => (TowerDouble -> f TowerDouble)+ -> Double+ -> f [Double]+diffsF f a = getADTower <$> apply f a+{-# INLINE diffsF #-}++-- | Compute the zero-padded derivatives of a function @(a -> f a)@+diffs0F+ :: Functor f+ => (TowerDouble -> f TowerDouble)+ -> Double+ -> f [Double]+diffs0F f a = zeroPad . getADTower <$> apply f a+{-# INLINE diffs0F #-}++-- | @taylor f x@ compute the Taylor series of @f@ around @x@.+taylor+ :: (TowerDouble -> TowerDouble)+ -> Double+ -> Double+ -> [Double]+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 f x@ compute the Taylor series of @f@ around @x@, zero-padded.+taylor0+ :: (TowerDouble -> TowerDouble)+ -> Double+ -> Double+ -> [Double]+taylor0 f x dx = zeroPad (taylor f x dx)+{-# INLINE taylor0 #-}++-- | @maclaurin f@ compute the Maclaurin series of @f@+maclaurin+ :: (TowerDouble -> TowerDouble)+ -> Double+ -> [Double]+maclaurin f = taylor f 0+{-# INLINE maclaurin #-}++-- | @maclaurin f@ compute the Maclaurin series of @f@, zero-padded+maclaurin0+ :: (TowerDouble -> TowerDouble)+ -> Double+ -> [Double]+maclaurin0 f = taylor0 f 0+{-# INLINE maclaurin0 #-}++-- | Compute the first derivative of a function @(a -> a)@+diff+ :: (TowerDouble -> TowerDouble)+ -> Double+ -> Double+diff f = d . diffs f+{-# INLINE diff #-}++-- | Compute the answer and first derivative of a function @(a -> a)@+diff'+ :: (TowerDouble -> TowerDouble)+ -> Double+ -> (Double, Double)+diff' f = d' . diffs f+{-# INLINE diff' #-}++-- | Compute a directional derivative of a function @(f a -> a)@+du+ :: Functor f+ => (f TowerDouble -> TowerDouble)+ -> f (Double, Double) -> Double+du f = d . getADTower . f . fmap withD+{-# INLINE du #-}++-- | Compute the answer and a directional derivative of a function @(f a -> a)@+du'+ :: Functor f+ => (f TowerDouble -> TowerDouble)+ -> f (Double, Double)+ -> (Double, Double)+du' f = d' . getADTower . f . fmap withD+{-# INLINE du' #-}++-- | Compute a directional derivative of a function @(f a -> g a)@+duF+ :: (Functor f, Functor g)+ => (f TowerDouble -> g TowerDouble)+ -> f (Double, Double)+ -> g Double+duF f = fmap (d . getADTower) . f . fmap withD+{-# INLINE duF #-}++-- | Compute the answer and a directional derivative of a function @(f a -> g a)@+duF'+ :: (Functor f, Functor g)+ => (f TowerDouble -> g TowerDouble)+ -> f (Double, Double)+ -> g (Double, Double)+duF' f = fmap (d' . getADTower) . f . fmap withD+{-# INLINE duF' #-}++-- | Given a function @(f a -> a)@, and a tower of derivatives, compute the corresponding directional derivatives.+dus+ :: Functor f+ => (f TowerDouble -> TowerDouble)+ -> f [Double]+ -> [Double]+dus f = getADTower . f . fmap tower+{-# INLINE dus #-}++-- | Given a function @(f a -> a)@, and a tower of derivatives, compute the corresponding directional derivatives, zero-padded+dus0+ :: Functor f+ => (f TowerDouble -> TowerDouble)+ -> f [Double]+ -> [Double]+dus0 f = zeroPad . getADTower . f . fmap tower+{-# INLINE dus0 #-}++-- | Given a function @(f a -> g a)@, and a tower of derivatives, compute the corresponding directional derivatives+dusF+ :: (Functor f, Functor g)+ => (f TowerDouble -> g TowerDouble)+ -> f [Double]+ -> g [Double]+dusF f = fmap getADTower . f . fmap tower+{-# INLINE dusF #-}++-- | Given a function @(f a -> g a)@, and a tower of derivatives, compute the corresponding directional derivatives, zero-padded+dus0F+ :: (Functor f, Functor g)+ => (f TowerDouble -> g TowerDouble)+ -> f [Double]+ -> g [Double]+dus0F f = fmap getADTower . f . fmap tower+{-# INLINE dus0F #-}
+ tests/Regression.hs view
@@ -0,0 +1,158 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE NoMonomorphismRestriction #-}+{-# LANGUAGE RankNTypes #-}++module Main (main) where++import Numeric+import qualified Numeric.AD.Mode.Forward as F+import qualified Numeric.AD.Mode.Forward.Double as FD+import qualified Numeric.AD.Mode.Reverse as R+import qualified Numeric.AD.Mode.Reverse.Double as RD++import Text.Printf+import Test.Tasty+import Test.Tasty.HUnit++type Diff' = (forall a. Floating a => a -> a) -> Double -> (Double, Double)+type Grad = (forall a. Floating a => [a] -> a) -> [Double] -> [Double]+type Jacobian = (forall a. Floating a => [a] -> [a]) -> [Double] -> [[Double]]+type Hessian = (forall a. Floating a => [a] -> a) -> [Double] -> [[Double]]++main :: IO ()+main = defaultMain tests++-- TODO: the forward-double tests are currently failing due to discrepancies between the modes+-- see also https://github.com/ekmett/ad/issues/109 and https://github.com/ekmett/ad/pull/110+tests :: TestTree+tests = testGroup "tests" [+ mode "forward" (\ f -> F.diff' f) (\ f -> F.grad f) (\ f -> F.jacobian f) (\ f -> F.jacobian $ F.grad f),+ --mode "forward-double" (\ f -> FD.diff' f) (\ f -> FD.grad f) (\ f -> FD.jacobian f) (\ f -> FD.jacobian $ F.grad f),+ mode "reverse" (\ f -> R.diff' f) (\ f -> R.grad f) (\ f -> R.jacobian f) (\ f -> R.hessian f),+ mode "reverse-double" (\ f -> RD.diff' f) (\ f -> RD.grad f) (\ f -> RD.jacobian f) (\ f -> RD.hessian f)]++mode :: String -> Diff' -> Grad -> Jacobian -> Hessian -> TestTree+mode name diff grad jacobian hessian = testGroup name [+ basic diff grad jacobian hessian,+ issue97 diff,+ issue104 diff grad,+ issue108 diff]++basic :: Diff' -> Grad -> Jacobian -> Hessian -> TestTree+basic diff grad jacobian hessian = testGroup "basic" [tdiff, tgrad, tjacobian, thessian] where+ tdiff = testCase "diff" $ do+ expect (list eq) [11, 5.5, 3, 3.5, 7, 13.5, 23, 35.5, 51] $ snd . diff p <$> [-2, -1.5, -1, -0.5, 0, 0.5, 1, 1.5, 2]+ expect (list eq) [nan, inf, 1, 0.5, 0.25] $ snd . diff sqrt <$> [-1, 0, 0.25, 1, 4]+ expect (list eq) [1, 0, 1] $ [snd . diff sin, snd . diff cos, snd . diff tan] <*> [0]+ expect (list eq) [-1, 0, 1] $ snd . diff abs <$> [-1, 0, 1]+ expect (list eq) [1, exp 1, inf, 1] $ [snd . diff exp, snd . diff log] <*> [0, 1]+ tgrad = testCase "grad" $ do+ expect (list eq) [2, 1, 1] $ grad f [1, 2, 3]+ expect (list eq) [1, 0.25] $ grad h [2, 8]+ expect (list eq) [0, nan] $ grad power [0, 2]+ tjacobian = testCase "jacobian" $ do+ expect (list $ list eq) [[0, 1], [1, 0], [1, 2]] $ jacobian g [2, 1]+ thessian = testCase "hessian" $ do+ expect (list $ list eq) [[0, 1, 0], [1, 0, 0], [0, 0, 0]] $ hessian f [1, 2, 3]+ expect (list $ list eq) [[0, 0], [0, 0]] $ hessian sum [1, 2]+ expect (list $ list eq) [[0, 1], [1, 0]] $ hessian product [1, 2]+ expect (list $ list eq) [[2, 1], [1, 0]] $ hessian power [1, 2]+ sum [x, y] = x + y+ product [x, y] = x * y+ power [x, y] = x ** y+ f [x, y, z] = x * y + z+ g [x, y] = [y, x, x * y]+ h [x, y] = sqrt $ x * y+ p x = 12 + 7 * x + 5 * x ^ 2 + 2 * x ^ 3++-- Reverse.Double +ffi initializes the tape with a block of size 4096+-- The large term in this function forces the allocation of an additional block+issue97 :: Diff' -> TestTree+issue97 diff = testCase "issue-97" $ expect eq 5000 $ snd $ diff f 0 where f = sum . replicate 5000++issue104 :: Diff' -> Grad -> TestTree+issue104 diff grad = testGroup "issue-104" [inside, outside] where+ inside = testGroup "inside" [tdiff, tgrad] where+ tdiff = testCase "diff" $ do+ expect (list eq) [nan, nan] $ snd . diff (0 `f`) <$> [0, 1]+ expect (list eq) [inf, 0.5] $ snd . diff (1 `f`) <$> [0, 1]+ expect (list eq) [nan, nan] $ snd . diff (`f` 0) <$> [0, 1]+ expect (list eq) [inf, 0.5] $ snd . diff (`f` 1) <$> [0, 1]+ tgrad = testCase "grad" $ do+ expect (list eq) [nan, nan] $ grad (binary f) [0, 0]+ expect (list eq) [nan, inf] $ grad (binary f) [1, 0]+ expect (list eq) [inf, nan] $ grad (binary f) [0, 1]+ expect (list eq) [0.5, 0.5] $ grad (binary f) [1, 1]+ f x y = sqrt $ x * y -- grad f [x, y] = [y / (2 * f x y), x / (2 * f x y)]+ outside = testGroup "outside" [tdiff, tgrad] where+ tdiff = testCase "diff" $ do+ expect (list eq) [nan, 0.0] $ snd . diff (0 `f`) <$> [0, 1]+ expect (list eq) [inf, 0.5] $ snd . diff (1 `f`) <$> [0, 1]+ expect (list eq) [nan, 0.0] $ snd . diff (`f` 0) <$> [0, 1]+ expect (list eq) [inf, 0.5] $ snd . diff (`f` 1) <$> [0, 1]+ tgrad = testCase "grad" $ do+ expect (list eq) [nan, nan] $ grad (binary f) [0, 0]+ expect (list eq) [0.0, inf] $ grad (binary f) [1, 0]+ expect (list eq) [inf, 0.0] $ grad (binary f) [0, 1]+ expect (list eq) [0.5, 0.5] $ grad (binary f) [1, 1]+ f x y = sqrt x * sqrt y -- grad f [x, y] = [sqrt y / 2 sqrt x, sqrt x / 2 sqrt y]+ binary f [x, y] = f x y++issue108 :: Diff' -> TestTree+issue108 diff = testGroup "issue-108" [tlog1p, texpm1, tlog1pexp, tlog1mexp] where+ tlog1p = testCase "log1p" $ do+ equal (-inf, inf) $ diff log1p (-1)+ equal (-1.0000000000000007e-15, 1.000000000000001) $ diff log1p (-1e-15)+ equal (-1e-20, 1) $ diff log1p (-1e-20)+ equal (0, 1) $ diff log1p 0+ equal (1e-20, 1) $ diff log1p 1e-20+ equal (9.999999999999995e-16, 0.9999999999999989) $ diff log1p 1e-15+ equal (0.6931471805599453, 0.5) $ diff log1p 1+ texpm1 = testCase "expm1" $ do+ equal (-0.6321205588285577, 0.36787944117144233) $ diff expm1 (-1)+ equal (-9.999999999999995e-16, 0.999999999999999) $ diff expm1 (-1e-15)+ equal (-1e-20, 1) $ diff expm1 (-1e-20)+ equal (0, 1) $ diff expm1 0+ equal (1e-20, 1) $ diff expm1 1e-20+ equal (1.0000000000000007e-15, 1.000000000000001) $ diff expm1 1e-15+ equal (1.718281828459045, 2.718281828459045) $ diff expm1 1+ tlog1pexp = testCase "log1pexp" $ do+ equal (0, 0) $ diff log1pexp (-1000)+ equal (3.720075976020836e-44, 3.7200759760208356e-44) $ diff log1pexp (-100)+ equal (0.31326168751822286, 0.2689414213699951) $ diff log1pexp (-1)+ equal (0.6931471805599453, 0.5) $ diff log1pexp 0+ equal (1.3132616875182228, 0.7310585786300049) $ diff log1pexp 1+ equal (100, 1) $ diff log1pexp 100+ equal (1000, 1) $ diff log1pexp 1000+ tlog1mexp = testCase "log1mexp" $ do+ equal (-0, -0) $ diff log1mexp (-1000)+-- old versions of base have a faulty implementation of log1mexp, causing this case to fail+-- see also https://gitlab.haskell.org/ghc/ghc/-/issues/17125+#if MIN_VERSION_base(4, 13, 0)+ equal (-3.720075976020836e-44, -3.7200759760208356e-44) $ diff log1mexp (-100)+#endif+ equal (-0.45867514538708193, -0.5819767068693265) $ diff log1mexp (-1)+ equal (-0.9327521295671886, -1.5414940825367982) $ diff log1mexp (-0.5)+ equal (-2.3521684610440907, -9.50833194477505) $ diff log1mexp (-0.1)+ equal (-34.538776394910684, -9.999999999999994e14) $ diff log1mexp (-1e-15)+ equal (-46.051701859880914, -1e20) $ diff log1mexp (-1e-20)+ equal (-inf, -inf) $ diff log1mexp (-0)+ equal = expect $ \ (a, b) (c, d) -> eq a c && eq b d++-- TODO: ideally, we would consider `0` and `-0` to be different+-- however, zero signedness is currently not reliably propagated through some modes+-- see also https://github.com/ekmett/ad/issues/109 and https://github.com/ekmett/ad/pull/110+eq :: Double -> Double -> Bool+eq a b = isNaN a && isNaN b || a == b++list :: (a -> a -> Bool) -> [a] -> [a] -> Bool+list eq as bs = length as == length bs && and (zipWith eq as bs)++expect :: HasCallStack => Show a => (a -> a -> Bool) -> a -> a -> Assertion+expect eq a b = eq a b @? printf "expected %s but got %s" (show a) (show b)++nan :: Double+nan = 0 / 0++inf :: Double+inf = 1 / 0
− tests/doctests.hs
@@ -1,25 +0,0 @@--------------------------------------------------------------------------------- |--- Module : Main (doctests)--- Copyright : (C) 2012-14 Edward Kmett--- License : BSD-style (see the file LICENSE)--- Maintainer : Edward Kmett <ekmett@gmail.com>--- Stability : provisional--- Portability : portable------ This module provides doctests for a project based on the actual versions--- of the packages it was built with. It requires a corresponding Setup.lhs--- to be added to the project-------------------------------------------------------------------------------module Main where--import Build_doctests (flags, pkgs, module_sources)-import Data.Foldable (traverse_)-import Test.DocTest--main :: IO ()-main = do- traverse_ putStrLn args- doctest args- where- args = flags ++ pkgs ++ module_sources
− travis/cabal-apt-install
@@ -1,27 +0,0 @@-#! /bin/bash-set -eu--APT="sudo apt-get -q -y"-CABAL_INSTALL_DEPS="cabal install --only-dependencies --force-reinstall"--$APT update-$APT install dctrl-tools--# Find potential system packages to satisfy cabal dependencies-deps()-{- local M='^\([^ ]\+\)-[0-9.]\+ (.*$'- local G=' -o ( -FPackage -X libghc-\L\1\E-dev )'- local E="$($CABAL_INSTALL_DEPS "$@" --dry-run -v 2> /dev/null \- | sed -ne "s/$M/$G/p" | sort -u)"- grep-aptavail -n -sPackage \( -FNone -X None \) $E | sort -u-}--$APT install $(deps "$@") libghc-quickcheck2-dev # QuickCheck is special-$CABAL_INSTALL_DEPS "$@" # Install the rest via Hackage--if ! $APT install hlint ; then- $APT install $(deps hlint)- cabal install hlint-fi-
− travis/config
@@ -1,16 +0,0 @@--- This provides a custom ~/.cabal/config file for use when hackage is down that should work on unix------ This is particularly useful for travis-ci to get it to stop complaining--- about a broken build when everything is still correct on our end.------ This uses Luite Stegeman's mirror of hackage provided by his 'hdiff' site instead------ To enable this, uncomment the before_script in .travis.yml--remote-repo: hdiff.luite.com:http://hdiff.luite.com/packages/archive-remote-repo-cache: ~/.cabal/packages-world-file: ~/.cabal/world-build-summary: ~/.cabal/logs/build.log-remote-build-reporting: anonymous-install-dirs user-install-dirs global