ad 4.4.1 → 4.5

raw patch · 71 files changed

+4218/−1437 lines, 71 filesdep +adjunctionsdep −directorydep −doctestdep −filepathdep ~arraydep ~basedep ~containerssetup-changedPVP ok

version bump matches the API change (PVP)

Dependencies added: adjunctions

Dependencies removed: directory, doctest, filepath

Dependency ranges changed: array, base, containers, transformers

API changes (from Hackage documentation)

Files

• .gitignore +1/−0
• .hlint.yaml +5/−0
• .travis.yml +0/−164
• CHANGELOG.markdown +33/−0
• README.markdown +1/−1
• Setup.lhs +5/−32
• Warning.hs +0/−5
• ad.cabal +53/−47
• cbits/tape.c +125/−0
• include/instances.h +23/−20
• include/internal_kahn.h +326/−0
• include/rank1_kahn.h +270/−0
• src/Numeric/AD.hs +34/−9
• 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 +12/−19
• src/Numeric/AD/Internal/Dense/Representable.hs +177/−0
• src/Numeric/AD/Internal/Doctest.hs +1/−1
• src/Numeric/AD/Internal/Forward.hs +6/−18
• src/Numeric/AD/Internal/Forward/Double.hs +12/−100
• src/Numeric/AD/Internal/Identity.hs +6/−5
• src/Numeric/AD/Internal/Kahn.hs +7/−13
• src/Numeric/AD/Internal/Kahn/Double.hs +9/−0
• src/Numeric/AD/Internal/Kahn/Float.hs +9/−0
• src/Numeric/AD/Internal/On.hs +4/−5
• src/Numeric/AD/Internal/Or.hs +3/−33
• src/Numeric/AD/Internal/Reverse.hs +17/−23
• src/Numeric/AD/Internal/Reverse/Double.hs +148/−170
• src/Numeric/AD/Internal/Sparse.hs +14/−59
• src/Numeric/AD/Internal/Sparse/Common.hs +54/−0
• src/Numeric/AD/Internal/Sparse/Double.hs +257/−0
• src/Numeric/AD/Internal/Tower.hs +10/−12
• src/Numeric/AD/Internal/Tower/Double.hs +270/−0
• src/Numeric/AD/Internal/Type.hs +2/−1
• src/Numeric/AD/Jacobian.hs +1/−2
• src/Numeric/AD/Jet.hs +8/−31
• src/Numeric/AD/Mode.hs +37/−34
• src/Numeric/AD/Mode/Dense.hs +68/−0
• src/Numeric/AD/Mode/Dense/Representable.hs +79/−0
• src/Numeric/AD/Mode/Forward.hs +1/−5
• src/Numeric/AD/Mode/Forward/Double.hs +1/−5
• src/Numeric/AD/Mode/Kahn.hs +1/−5
• src/Numeric/AD/Mode/Kahn/Double.hs +240/−0
• src/Numeric/AD/Mode/Reverse.hs +95/−33
• src/Numeric/AD/Mode/Reverse/Double.hs +91/−33
• 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 +1/−6
• 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 +105/−30
• src/Numeric/AD/Rank1/Forward/Double.hs +103/−24
• src/Numeric/AD/Rank1/Halley.hs +4/−2
• src/Numeric/AD/Rank1/Halley/Double.hs +122/−0
• src/Numeric/AD/Rank1/Kahn.hs +13/−215
• 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/doctests.hs +0/−25
• travis/cabal-apt-install +0/−27
• travis/config +0/−16

@@ -31,3 +31,4 @@ cabal.project.local~ .HTF/ .ghc.environment.*+Makefile

@@ -3,8 +3,13 @@ - 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}

@@ -1,164 +0,0 @@-# This Travis job script has been generated by a script via-#-#   haskell-ci '--output=.travis.yml' '--config=cabal.haskell-ci' 'cabal.project'-#-# To regenerate the script (for example after adjusting tested-with) run-#-#   haskell-ci regenerate-#-# For more information, see https://github.com/haskell-CI/haskell-ci-#-# version: 0.10-#-version: ~> 1.0-language: c-os: linux-dist: xenial-git:-  # whether to recursively clone submodules-  submodules: false-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-    - $HOME/.hlint-before_cache:-  - rm -fv $CABALHOME/packages/hackage.haskell.org/build-reports.log-  # remove files that are regenerated by 'cabal update'-  - rm -fv $CABALHOME/packages/hackage.haskell.org/00-index.*-  - rm -fv $CABALHOME/packages/hackage.haskell.org/*.json-  - rm -fv $CABALHOME/packages/hackage.haskell.org/01-index.cache-  - rm -fv $CABALHOME/packages/hackage.haskell.org/01-index.tar-  - rm -fv $CABALHOME/packages/hackage.haskell.org/01-index.tar.idx-  - rm -rfv $CABALHOME/packages/head.hackage-jobs:-  include:-    - compiler: ghc-8.10.1-      addons: {"apt":{"sources":[{"sourceline":"deb http://ppa.launchpad.net/hvr/ghc/ubuntu xenial main","key_url":"https://keyserver.ubuntu.com/pks/lookup?op=get&search=0x063dab2bdc0b3f9fcebc378bff3aeacef6f88286"}],"packages":["ghc-8.10.1","cabal-install-3.2"]}}-      os: linux-    - compiler: ghc-8.8.3-      addons: {"apt":{"sources":[{"sourceline":"deb http://ppa.launchpad.net/hvr/ghc/ubuntu xenial main","key_url":"https://keyserver.ubuntu.com/pks/lookup?op=get&search=0x063dab2bdc0b3f9fcebc378bff3aeacef6f88286"}],"packages":["ghc-8.8.3","cabal-install-3.2"]}}-      os: linux-    - compiler: ghc-8.6.5-      addons: {"apt":{"sources":[{"sourceline":"deb http://ppa.launchpad.net/hvr/ghc/ubuntu xenial main","key_url":"https://keyserver.ubuntu.com/pks/lookup?op=get&search=0x063dab2bdc0b3f9fcebc378bff3aeacef6f88286"}],"packages":["ghc-8.6.5","cabal-install-3.2"]}}-      os: linux-    - compiler: ghc-8.4.4-      addons: {"apt":{"sources":[{"sourceline":"deb http://ppa.launchpad.net/hvr/ghc/ubuntu xenial main","key_url":"https://keyserver.ubuntu.com/pks/lookup?op=get&search=0x063dab2bdc0b3f9fcebc378bff3aeacef6f88286"}],"packages":["ghc-8.4.4","cabal-install-3.2"]}}-      os: linux-    - compiler: ghc-8.2.2-      addons: {"apt":{"sources":[{"sourceline":"deb http://ppa.launchpad.net/hvr/ghc/ubuntu xenial main","key_url":"https://keyserver.ubuntu.com/pks/lookup?op=get&search=0x063dab2bdc0b3f9fcebc378bff3aeacef6f88286"}],"packages":["ghc-8.2.2","cabal-install-3.2"]}}-      os: linux-    - compiler: ghc-8.0.2-      addons: {"apt":{"sources":[{"sourceline":"deb http://ppa.launchpad.net/hvr/ghc/ubuntu xenial main","key_url":"https://keyserver.ubuntu.com/pks/lookup?op=get&search=0x063dab2bdc0b3f9fcebc378bff3aeacef6f88286"}],"packages":["ghc-8.0.2","cabal-install-3.2"]}}-      os: linux-    - compiler: ghc-7.10.3-      addons: {"apt":{"sources":[{"sourceline":"deb http://ppa.launchpad.net/hvr/ghc/ubuntu xenial main","key_url":"https://keyserver.ubuntu.com/pks/lookup?op=get&search=0x063dab2bdc0b3f9fcebc378bff3aeacef6f88286"}],"packages":["ghc-7.10.3","cabal-install-3.2"]}}-      os: linux-    - compiler: ghc-7.8.4-      addons: {"apt":{"sources":[{"sourceline":"deb http://ppa.launchpad.net/hvr/ghc/ubuntu xenial main","key_url":"https://keyserver.ubuntu.com/pks/lookup?op=get&search=0x063dab2bdc0b3f9fcebc378bff3aeacef6f88286"}],"packages":["ghc-7.8.4","cabal-install-3.2"]}}-      os: linux-    - compiler: ghc-7.6.3-      addons: {"apt":{"sources":[{"sourceline":"deb http://ppa.launchpad.net/hvr/ghc/ubuntu xenial main","key_url":"https://keyserver.ubuntu.com/pks/lookup?op=get&search=0x063dab2bdc0b3f9fcebc378bff3aeacef6f88286"}],"packages":["ghc-7.6.3","cabal-install-3.2"]}}-      os: linux-    - compiler: ghc-7.4.2-      addons: {"apt":{"sources":[{"sourceline":"deb http://ppa.launchpad.net/hvr/ghc/ubuntu xenial main","key_url":"https://keyserver.ubuntu.com/pks/lookup?op=get&search=0x063dab2bdc0b3f9fcebc378bff3aeacef6f88286"}],"packages":["ghc-7.4.2","cabal-install-3.2"]}}-      os: linux-before_install:-  - HC=$(echo "/opt/$CC/bin/ghc" | sed 's/-/\//')-  - WITHCOMPILER="-w $HC"-  - HADDOCK=$(echo "/opt/$CC/bin/haddock" | sed 's/-/\//')-  - HCPKG="$HC-pkg"-  - unset CC-  - CABAL=/opt/ghc/bin/cabal-  - CABALHOME=$HOME/.cabal-  - export PATH="$CABALHOME/bin:$PATH"-  - TOP=$(pwd)-  - "HCNUMVER=$(${HC} --numeric-version|perl -ne '/^(\\d+)\\.(\\d+)\\.(\\d+)(\\.(\\d+))?$/; print(10000 * $1 + 100 * $2 + ($3 == 0 ? $5 != 1 : $3))')"-  - echo $HCNUMVER-  - CABAL="$CABAL -vnormal+nowrap"-  - set -o pipefail-  - TEST=--enable-tests-  - BENCH=--enable-benchmarks-  - HEADHACKAGE=false-  - rm -f $CABALHOME/config-  - |-    echo "verbose: normal +nowrap +markoutput"          >> $CABALHOME/config-    echo "remote-build-reporting: anonymous"            >> $CABALHOME/config-    echo "write-ghc-environment-files: always"          >> $CABALHOME/config-    echo "remote-repo-cache: $CABALHOME/packages"       >> $CABALHOME/config-    echo "logs-dir:          $CABALHOME/logs"           >> $CABALHOME/config-    echo "world-file:        $CABALHOME/world"          >> $CABALHOME/config-    echo "extra-prog-path:   $CABALHOME/bin"            >> $CABALHOME/config-    echo "symlink-bindir:    $CABALHOME/bin"            >> $CABALHOME/config-    echo "installdir:        $CABALHOME/bin"            >> $CABALHOME/config-    echo "build-summary:     $CABALHOME/logs/build.log" >> $CABALHOME/config-    echo "store-dir:         $CABALHOME/store"          >> $CABALHOME/config-    echo "install-dirs user"                            >> $CABALHOME/config-    echo "  prefix: $CABALHOME"                         >> $CABALHOME/config-    echo "repository hackage.haskell.org"               >> $CABALHOME/config-    echo "  url: http://hackage.haskell.org/"           >> $CABALHOME/config-install:-  - ${CABAL} --version-  - echo "$(${HC} --version) [$(${HC} --print-project-git-commit-id 2> /dev/null || echo '?')]"-  - |-    echo "program-default-options"                >> $CABALHOME/config-    echo "  ghc-options: $GHCJOBS +RTS -M6G -RTS" >> $CABALHOME/config-  - cat $CABALHOME/config-  - rm -fv cabal.project cabal.project.local cabal.project.freeze-  - travis_retry ${CABAL} v2-update -v-  # Generate cabal.project-  - rm -rf cabal.project cabal.project.local cabal.project.freeze-  - touch cabal.project-  - |-    echo "packages: ." >> cabal.project-  - if [ $HCNUMVER -ge 80200 ] ; then echo 'package ad' >> cabal.project ; fi-  - "if [ $HCNUMVER -ge 80200 ] ; then echo '  ghc-options: -Werror=missing-methods' >> cabal.project ; fi"-  - |-  - "for pkg in $($HCPKG list --simple-output); do echo $pkg | sed 's/-[^-]*$//' | (grep -vE -- '^(ad)$' || true) | sed 's/^/constraints: /' | sed 's/$/ installed/' >> cabal.project.local; done"-  - cat cabal.project || true-  - cat cabal.project.local || true-  - if [ -f "./configure.ac" ]; then (cd "." && autoreconf -i); fi-  - ${CABAL} v2-freeze $WITHCOMPILER ${TEST} ${BENCH}-  - "cat cabal.project.freeze | sed -E 's/^(constraints: *| *)//' | sed 's/any.//'"-  - rm  cabal.project.freeze-  - travis_wait 40 ${CABAL} v2-build $WITHCOMPILER ${TEST} ${BENCH} --dep -j2 all-script:-  - DISTDIR=$(mktemp -d /tmp/dist-test.XXXX)-  # Packaging...-  - ${CABAL} v2-sdist all-  # Unpacking...-  - mv dist-newstyle/sdist/*.tar.gz ${DISTDIR}/-  - cd ${DISTDIR} || false-  - find . -maxdepth 1 -type f -name '*.tar.gz' -exec tar -xvf '{}' \;-  - find . -maxdepth 1 -type f -name '*.tar.gz' -exec rm       '{}' \;-  - PKGDIR_ad="$(find . -maxdepth 1 -type d -regex '.*/ad-[0-9.]*')"-  # Generate cabal.project-  - rm -rf cabal.project cabal.project.local cabal.project.freeze-  - touch cabal.project-  - |-    echo "packages: ${PKGDIR_ad}" >> cabal.project-  - if [ $HCNUMVER -ge 80200 ] ; then echo 'package ad' >> cabal.project ; fi-  - "if [ $HCNUMVER -ge 80200 ] ; then echo '  ghc-options: -Werror=missing-methods' >> cabal.project ; fi"-  - |-  - "for pkg in $($HCPKG list --simple-output); do echo $pkg | sed 's/-[^-]*$//' | (grep -vE -- '^(ad)$' || true) | sed 's/^/constraints: /' | sed 's/$/ installed/' >> cabal.project.local; done"-  - cat cabal.project || true-  - cat cabal.project.local || true-  # Building with tests and benchmarks...-  # build & run tests, build benchmarks-  - ${CABAL} v2-build $WITHCOMPILER ${TEST} ${BENCH} all-  # Testing...-  - ${CABAL} v2-test $WITHCOMPILER ${TEST} ${BENCH} all-  # cabal check...-  - (cd ${PKGDIR_ad} && ${CABAL} -vnormal check)-  # haddock...-  - ${CABAL} v2-haddock $WITHCOMPILER --with-haddock $HADDOCK ${TEST} ${BENCH} all--# REGENDATA ("0.10",["--output=.travis.yml","--config=cabal.haskell-ci","cabal.project"])-# EOF

@@ -1,3 +1,36 @@+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.

@@ -1,7 +1,7 @@ ad == -[![Hackage](https://img.shields.io/hackage/v/ad.svg)](https://hackage.haskell.org/package/ad) [![Build Status](https://secure.travis-ci.org/ekmett/ad.png?branch=master)](http://travis-ci.org/ekmett/ad)+[![Hackage](https://img.shields.io/hackage/v/ad.svg)](https://hackage.haskell.org/package/ad) [![Build Status](https://github.com/ekmett/ad/workflows/Haskell-CI/badge.svg)](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. 

@@ -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

@@ -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

@@ -1,8 +1,8 @@ name:          ad-version:       4.4.1+version:       4.5 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,30 +10,26 @@ 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.5-             , GHC == 8.8.3-             , GHC == 8.10.1+             , GHC == 8.8.4+             , GHC == 8.10.7+             , GHC == 9.0.1+             , GHC == 9.2.1 synopsis:      Automatic Differentiation extra-source-files:   .gitignore   .hlint.yaml-  .travis.yml   .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.     .@@ -91,17 +87,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@@ -111,6 +104,7 @@     GeneralizedNewtypeDeriving     MultiParamTypeClasses     PatternGuards+    PatternSynonyms     Rank2Types     ScopedTypeVariables     TypeFamilies@@ -118,17 +112,18 @@     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.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.5.2.0 && < 0.6    if impl(ghc < 7.8)     build-depends: tagged >= 0.7 && < 1@@ -140,42 +135,66 @@    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@@ -183,20 +202,7 @@     ghc-options: -Wno-star-is-type    ghc-options: -fspec-constr -fdicts-cheap -O2---- Verify the results of the examples-test-suite doctests-  default-language: Haskell2010-  type:    exitcode-stdio-1.0-  main-is: doctests.hs-  build-depends:-    ad,-    base,-    directory,-    doctest >= 0.9.0.1 && < 0.18,-    filepath-  ghc-options: -Wall -threaded-  hs-source-dirs: tests+  x-docspec-extra-packages: distributive  benchmark blackscholes   default-language: Haskell2010

@@ -0,0 +1,125 @@+#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];+      if (v == 0.0) continue;++      int i = pTape->lnk[idx*2];+      double x = pTape->val[idx*2];+      if (x != 0.0)+      {+        buffer[i] += v*x;+      }++      int j = pTape->lnk[idx*2 + 1]; +      double y = pTape->val[idx*2 + 1];+      if (y != 0.0)+      {+        buffer[j] += v*y;+      }+    }+    idx += 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);+  }+}

@@ -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,7 @@   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+instance BODY2(Num a, Enum a) Enum HEAD where   succ             = lift1 succ (const 1)   pred             = lift1 pred (const 1)   toEnum           = auto . toEnum@@ -66,10 +69,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 +88,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 +97,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)

@@ -0,0 +1,326 @@++{-# 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 Data.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 GHC.Exts as Exts+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) = 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 IsList List where+  type Item List = SCALAR_TYPE+  fromList (x:xs) = Cons x (fromList xs)+  fromList [] = Nil+  toList Nil = []+  toList (Cons x xs) = x : 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

@@ -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!

@@ -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@@ -138,10 +137,8 @@   ) 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@@ -189,7 +186,14 @@ -- -- 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.@@ -197,21 +201,42 @@ -- > 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

@@ -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!

@@ -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 #-}

@@ -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 #-}

@@ -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

@@ -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,16 +41,8 @@   ) 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.AD.Internal.Combinators@@ -97,6 +89,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 +116,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 +181,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"

@@ -0,0 +1,177 @@+{-# 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.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"

@@ -2,7 +2,7 @@  ----------------------------------------------------------------------------- -- |--- Copyright   :  (c) Edward Kmett 2010-2019+-- Copyright   :  (c) Edward Kmett 2010-2021 -- License     :  BSD3 -- Maintainer  :  ekmett@gmail.com -- Stability   :  experimental

@@ -10,7 +10,7 @@  ----------------------------------------------------------------------------- -- |--- Copyright   :  (c) Edward Kmett 2010-2015+-- Copyright   :  (c) Edward Kmett 2010-2021 -- License     :  BSD3 -- Maintainer  :  ekmett@gmail.com -- Stability   :  experimental@@ -37,14 +37,8 @@   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.AD.Internal.Combinators@@ -52,10 +46,6 @@ 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@@ -97,6 +87,10 @@   isKnownZero Zero = True   isKnownZero _    = False +  asKnownConstant Zero = Just 0+  asKnownConstant (Lift a) = Just a+  asKnownConstant _ = Nothing+   isKnownConstant Forward{} = False   isKnownConstant _ = True @@ -120,12 +114,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 @@ -180,7 +168,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

@@ -10,7 +10,7 @@  ----------------------------------------------------------------------------- ---- |----- Copyright   :  (c) Edward Kmett 2010-2015+---- Copyright   :  (c) Edward Kmett 2010-2021 ---- License     :  BSD3 ---- Maintainer  :  ekmett@gmail.com ---- Stability   :  experimental@@ -33,16 +33,9 @@   , 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.AD.Internal.Combinators import Numeric.AD.Internal.Identity@@ -74,6 +67,9 @@   isKnownZero (ForwardDouble 0 0) = True   isKnownZero _ = False +  asKnownConstant (ForwardDouble x 0) = Just x+  asKnownConstant _ = Nothing+    isKnownConstant (ForwardDouble _ 0) = True   isKnownConstant _ = False @@ -108,98 +104,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 -> (recip $ sqrt (2 * pi)) * exp (- x * x / 2)--instance InvErf ForwardDouble 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)+#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@@ -229,3 +138,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))

@@ -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

@@ -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,11 +48,6 @@   , 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@@ -99,12 +94,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 +112,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 +156,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

@@ -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>

@@ -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>

@@ -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 ------------------------------------------------------------------------------@@ -44,5 +40,8 @@ instance (Mode t, Mode (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)

@@ -3,14 +3,12 @@ {-# LANGUAGE TypeFamilies #-} {-# 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 +26,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 +51,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 +69,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 +89,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 +177,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)

@@ -13,7 +13,7 @@  ----------------------------------------------------------------------------- -- |--- Copyright   :  (c) Edward Kmett 2012-2015+-- Copyright   :  (c) Edward Kmett 2012-2021 -- License     :  BSD3 -- Maintainer  :  ekmett@gmail.com -- Stability   :  experimental@@ -36,6 +36,7 @@   , Head(..)   , Cells(..)   , reifyTape+  , reifyTypeableTape   , partials   , partialArrayOf   , partialMapOf@@ -63,11 +64,7 @@ import Data.Number.Erf import Data.Proxy import Data.Reflection-#if __GLASGOW_HASKELL__ < 710-import Data.Traversable (Traversable, mapM)-#else import Data.Traversable (mapM)-#endif import Data.Typeable import Numeric.AD.Internal.Combinators import Numeric.AD.Internal.Identity@@ -77,17 +74,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@@ -127,13 +118,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@@ -141,6 +130,10 @@   isKnownZero Zero = True   isKnownZero _    = False +  asKnownConstant Zero = Just 0+  asKnownConstant (Lift n) = Just n+  asKnownConstant _ = Nothing+   isKnownConstant Reverse{} = False   isKnownConstant _ = True @@ -153,12 +146,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@@ -167,7 +154,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 @@ -201,9 +188,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.@@ -263,6 +250,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

@@ -1,4 +1,5 @@ {-# LANGUAGE CPP #-}+{-# LANGUAGE BangPatterns #-} {-# LANGUAGE DeriveDataTypeable #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-}@@ -11,12 +12,15 @@ {-# OPTIONS_GHC -fno-full-laziness #-} {-# OPTIONS_HADDOCK not-home #-} +#ifdef AD_FFI+{-# LANGUAGE ForeignFunctionInterface #-}+#endif+ module Numeric.AD.Internal.Reverse.Double   ( ReverseDouble(..)   , Tape(..)-  , Head(..)-  , Cells(..)   , reifyTape+  , reifyTypeableTape   , partials   , partialArrayOf   , partialMapOf@@ -32,23 +36,29 @@   , primal   ) where -import Data.Functor-import Control.Monad hiding (mapM)+#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 Control.Monad.Trans.State import Data.Array.ST-import Data.Array 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-#if __GLASGOW_HASKELL__ < 710-import Data.Traversable (Traversable, mapM)-#else import Data.Traversable (mapM)-#endif import Data.Typeable import Numeric.AD.Internal.Combinators import Numeric.AD.Internal.Identity@@ -56,19 +66,70 @@ import Numeric.AD.Mode import Prelude hiding (mapM) import System.IO.Unsafe (unsafePerformIO)-import Unsafe.Coerce -#ifdef HLINT-{-# ANN module "HLint: ignore Reduce duplication" #-}-#endif+#ifdef AD_FFI --- evil untyped tape-#ifndef HLINT+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 0 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 -> Double -> Cells -> Cells-  Binary :: {-# UNPACK #-} !Int -> {-# UNPACK #-} !Int -> Double -> Double -> Cells -> Cells-#endif+  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@@ -76,10 +137,54 @@ dropCells n (Unary _ _ xs)      = (dropCells $! n - 1) xs dropCells n (Binary _ _ _ _ xs) = (dropCells $! n - 1) xs -data Head = Head {-# UNPACK #-} !Int Cells+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@@ -108,18 +213,19 @@ binarily f di dj i b j c = ReverseDouble (unsafePerformIO (modifyTape (Proxy :: Proxy s) (bin i j di dj))) $! f b c {-# INLINE binarily #-} -#ifndef HLINT+#endif+ data ReverseDouble s where   Zero :: ReverseDouble s-  Lift :: Double -> ReverseDouble s-  ReverseDouble :: {-# UNPACK #-} !Int -> Double -> ReverseDouble s+  Lift :: {-# UNPACK #-} !Double -> ReverseDouble s+  ReverseDouble :: {-# UNPACK #-} !Int -> {-# UNPACK #-} !Double -> ReverseDouble s   deriving (Show, Typeable)-#endif -instance (Reifies s Tape) => Mode (ReverseDouble s) where+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@@ -131,24 +237,18 @@   a ^* b = lift1 (* b) (\_ -> auto b) a   a ^/ b = lift1 (/ b) (\_ -> auto (recip b)) a -(<+>) :: (Reifies s Tape) => ReverseDouble s -> ReverseDouble s -> ReverseDouble s+(<+>) :: Reifies s Tape => ReverseDouble s -> ReverseDouble s -> ReverseDouble s (<+>)  = binary (+) 1 1 -(<**>) :: (Reifies s Tape) => ReverseDouble s -> ReverseDouble s -> ReverseDouble s-Zero <**> y      = auto (0 ** primal y)-_    <**> Zero   = auto 1-x    <**> Lift y = lift1 (**y) (\z -> y *^ z ** Id (y - 1)) x-x    <**> y      = lift2_ (**) (\z xi yi -> (yi * xi ** (yi - 1), z * log xi)) x y- primal :: ReverseDouble s -> Double primal Zero = 0 primal (Lift a) = a primal (ReverseDouble _ a) = a -instance (Reifies s Tape) => Jacobian (ReverseDouble s) where+instance Reifies s Tape => Jacobian (ReverseDouble s) where   type D (ReverseDouble s) = Id Double -  unary f _         (Zero)   = Lift (f 0)+  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 @@ -179,158 +279,36 @@     a = f pb pc     (dadb, dadc) = df (Id a) (Id pb) (Id pc) --instance (Reifies s Tape) => Eq (ReverseDouble s) where-  a == b = primal a == primal b--instance (Reifies s Tape) => Ord (ReverseDouble s) where-  compare a b = compare (primal a) (primal b)--instance (Reifies s Tape) => Num (ReverseDouble s) where-  fromInteger 0  = zero-  fromInteger n = auto (fromInteger n)-  (+)          = (<+>) -- binary (+) 1 1-  (-)          = binary (-) (auto 1) (auto (-1)) -- TODO: <-> ? as it is, this might be pretty bad for Tower-  (*)          = lift2 (*) (\x y -> (y, x))-  negate       = lift1 negate (const (auto (-1)))-  abs          = lift1 abs signum-  signum a     = lift1 signum (const zero) a--instance (Reifies s Tape) => Fractional (ReverseDouble s) where-  fromRational 0 = zero-  fromRational r = auto (fromRational r)-  x / y        = x * recip y-  recip        = lift1_ recip (const . negate . join (*))--instance (Reifies s Tape) => Floating (ReverseDouble s) where-  pi       = auto pi-  exp      = lift1_ exp const-  log      = lift1 log recip-  logBase x y = log y / log x-  sqrt     = lift1_ sqrt (\z _ -> recip (auto 2 * z))-  (**)     = (<**>)-  --x ** y-  --   | isKnownZero y     = 1-  --   | isKnownConstant y, y' <- primal y = lift1 (** y') ((y'*) . (**(y'-1))) x-  --   | otherwise         = lift2_ (**) (\z xi yi -> (yi * z / xi, z * log1 xi)) x y-  sin      = lift1 sin cos-  cos      = lift1 cos $ negate . sin-  tan      = lift1 tan $ recip . join (*) . cos-  asin     = lift1 asin $ \x -> recip (sqrt (auto 1 - join (*) x))-  acos     = lift1 acos $ \x -> negate (recip (sqrt (1 - join (*) x)))-  atan     = lift1 atan $ \x -> recip (1 + join (*) x)-  sinh     = lift1 sinh cosh-  cosh     = lift1 cosh sinh-  tanh     = lift1 tanh $ recip . join (*) . cosh-  asinh    = lift1 asinh $ \x -> recip (sqrt (1 + join (*) x))-  acosh    = lift1 acosh $ \x -> recip (sqrt (join (*) x - 1))-  atanh    = lift1 atanh $ \x -> recip (1 - join (*) x)--instance (Reifies s Tape) => Enum (ReverseDouble s) where-  succ             = lift1 succ (const 1)-  pred             = lift1 pred (const 1)-  toEnum           = auto . toEnum-  fromEnum a       = fromEnum (primal a)-  enumFrom a       = withPrimal a <$> enumFrom (primal a)-  enumFromTo a b   = withPrimal a <$> enumFromTo (primal a) (primal b)-  enumFromThen a b = zipWith (fromBy a delta) [0..] $ enumFromThen (primal a) (primal b) where delta = b - a-  enumFromThenTo a b c = zipWith (fromBy a delta) [0..] $ enumFromThenTo (primal a) (primal b) (primal c) where delta = b - a--instance (Reifies s Tape) => Real (ReverseDouble s) where-  toRational      = toRational . primal--instance (Reifies s Tape) => RealFloat (ReverseDouble s) where-  floatRadix      = floatRadix . primal-  floatDigits     = floatDigits . primal-  floatRange      = floatRange . primal-  decodeFloat     = decodeFloat . primal-  encodeFloat m e = auto (encodeFloat m e)-  isNaN           = isNaN . primal-  isInfinite      = isInfinite . primal-  isDenormalized  = isDenormalized . primal-  isNegativeZero  = isNegativeZero . primal-  isIEEE          = isIEEE . primal-  exponent = exponent-  scaleFloat n = unary (scaleFloat n) (scaleFloat n 1)-  significand x =  unary significand (scaleFloat (- floatDigits x) 1) x-  atan2 = lift2 atan2 $ \vx vy -> let r = recip (join (*) vx + join (*) vy) in (vy * r, negate vx * r)--instance (Reifies s Tape) => RealFrac (ReverseDouble s) where-  properFraction a = (w, a `withPrimal` pb) where-    pa = primal a-    (w, pb) = properFraction pa-  truncate = truncate . primal-  round    = round . primal-  ceiling  = ceiling . primal-  floor    = floor . primal--instance (Reifies s Tape) => Erf (ReverseDouble s) where-  erf = lift1 erf $ \x -> (2 / sqrt pi) * exp (negate x * x)-  erfc = lift1 erfc $ \x -> ((-2) / sqrt pi) * exp (negate x * x)-  normcdf = lift1 normcdf $ \x -> (recip $ sqrt (2 * pi)) * exp (- x * x / 2)+mul :: Reifies s Tape => ReverseDouble s -> ReverseDouble s -> ReverseDouble s+mul = lift2 (*) (\x y -> (y, x)) -instance (Reifies s Tape) => InvErf (ReverseDouble s) where-  inverf = lift1_ inverf $ \x _ -> sqrt pi / 2 * exp (x * x)-  inverfc = lift1_ inverfc $ \x _ -> negate (sqrt pi / 2) * exp (x * x)-  invnormcdf = lift1_ invnormcdf $ \x _ -> sqrt (2 * pi) * exp (x * x / 2)+#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 :: 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' :: Reifies s Tape => Proxy s -> ReverseDouble s -> (Double, Double) derivativeOf' p r = (primal r, derivativeOf p r) {-# INLINE derivativeOf' #-} --- | Used internally to push sensitivities down the chain.-backPropagate :: Int -> Cells -> STArray s Int Double -> ST s Int-backPropagate k Nil _ = return k-backPropagate k (Unary i g xs) ss = do-  da <- readArray ss k-  db <- readArray ss i-  writeArray ss i $! db + unsafeCoerce g*da-  (backPropagate $! k - 1) xs ss-backPropagate k (Binary i j g h xs) ss = do-  da <- readArray ss k-  db <- readArray ss i-  writeArray ss i $! db + unsafeCoerce g*da-  dc <- readArray ss j-  writeArray ss j $! dc + unsafeCoerce h*da-  (backPropagate $! k - 1) xs ss --- | Extract the partials from the current chain for a given AD variable.-partials :: forall s. (Reifies s Tape) => ReverseDouble s -> [Double]-partials Zero        = []-partials (Lift _)    = []-partials (ReverseDouble k _) = map (sensitivities !) [0..vs] where-  Head n t = unsafePerformIO $ readIORef (getTape (reflect (Proxy :: Proxy s)))-  tk = dropCells (n - k) t-  (vs,sensitivities) = runST $ do-    ss <- newArray (0, k) 0-    writeArray ss k 1-    v <- backPropagate k tk ss-    as <- Unsafe.unsafeFreeze ss-    return (v, as)- -- | Return an 'Array' of 'partials' given bounds for the variable IDs.-partialArrayOf :: (Reifies s Tape) => Proxy s -> (Int, Int) -> ReverseDouble s -> Array Int Double+partialArrayOf :: 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 :: Reifies s Tape => Proxy s -> ReverseDouble s-> IntMap Double partialMapOf _ = fromDistinctAscList . zip [0..] . partials {-# INLINE partialMapOf #-} --- | Construct a tape that starts with @n@ variables.-reifyTape :: Int -> (forall s. Reifies s Tape => Proxy s -> r) -> r-reifyTape vs k = unsafePerformIO $ do-  h <- newIORef (Head vs Nil)-  return (reify (Tape h) k)-{-# NOINLINE reifyTape #-}- var :: Double -> Int -> ReverseDouble s var a v = ReverseDouble v a @@ -346,11 +324,11 @@ 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 :: 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 :: 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 :: 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

@@ -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,19 @@   ) 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.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 +76,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 +88,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 +119,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 +155,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 +178,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 +228,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

@@ -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"

@@ -0,0 +1,257 @@+{-# 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.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)

@@ -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,9 +35,6 @@   ) 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)@@ -60,7 +57,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 +118,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 +162,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 +179,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>

@@ -0,0 +1,270 @@+{-# 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.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>

@@ -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

@@ -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

@@ -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,19 +16,10 @@   , 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 @@ -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

@@ -1,16 +1,18 @@-{-# 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 UndecidableInstances #-} ----------------------------------------------------------------------------- -- |--- Copyright   :  (c) Edward Kmett 2010-2015+-- Copyright   :  (c) Edward Kmett 2010-2021 -- License     :  BSD3 -- Maintainer  :  ekmett@gmail.com -- Stability   :  experimental@@ -22,6 +24,8 @@   (   -- * AD modes     Mode(..)+  , pattern KnownZero+  , pattern Auto   ) where  import Numeric.Natural@@ -36,16 +40,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 +75,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   (^/) = (/)

@@ -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' #-}

@@ -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' #-}

@@ -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

@@ -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

@@ -1,4 +1,3 @@-{-# LANGUAGE CPP #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE FunctionalDependencies #-} {-# LANGUAGE MultiParamTypeClasses #-}@@ -7,7 +6,7 @@ {-# LANGUAGE UndecidableInstances #-} ----------------------------------------------------------------------------- -- |--- Copyright   :  (c) Edward Kmett 2010-2015+-- Copyright   :  (c) Edward Kmett 2010-2021 -- License     :  BSD3 -- Maintainer  :  ekmett@gmail.com -- Stability   :  experimental@@ -45,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(..))

@@ -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)

@@ -1,4 +1,3 @@-{-# LANGUAGE CPP #-} {-# LANGUAGE Rank2Types #-} {-# LANGUAGE BangPatterns #-} {-# LANGUAGE MultiParamTypeClasses #-}@@ -8,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@@ -44,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@@ -66,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 #-} @@ -75,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' #-}@@ -88,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 #-} @@ -99,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' #-}@@ -109,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 #-} @@ -120,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@@ -135,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 #-} @@ -147,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@@ -158,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.@@ -169,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.@@ -196,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 #-} @@ -206,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 #-}

@@ -1,4 +1,3 @@-{-# LANGUAGE CPP #-} {-# LANGUAGE Rank2Types #-} {-# LANGUAGE BangPatterns #-} {-# LANGUAGE MultiParamTypeClasses #-}@@ -8,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@@ -18,7 +17,7 @@ -- Data.Reflection -- -- This version is specialized to `Double` enabling the entire--- structure+-- structure to be unboxed. -- ----------------------------------------------------------------------------- @@ -47,10 +46,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@@ -71,8 +67,12 @@ -- -- >>> grad (\[x,y] -> x**y) [0,2] -- [0.0,NaN]-grad :: (Traversable f) => (forall s. Reifies s Tape => f (ReverseDouble s) -> ReverseDouble s) -> f Double -> f Double-grad f as = reifyTape (snd bds) $ \p -> unbind vs $! partialArrayOf p bds $! f vs where+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 #-} @@ -80,8 +80,12 @@ -- -- >>> grad' (\[x,y,z] -> x*y+z) [1,2,3] -- (5.0,[2.0,1.0,1.0])-grad' :: (Traversable f) => (forall s. Reifies s Tape => f (ReverseDouble s) -> ReverseDouble s) -> f Double -> (Double, f Double)-grad' f as = reifyTape (snd bds) $ \p -> case f vs of+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' #-}@@ -93,8 +97,13 @@ -- 'grad' == 'gradWith' (\_ dx -> dx) -- 'id' == 'gradWith' 'const' -- @-gradWith :: (Traversable f) => (Double -> Double -> b) -> (forall s. Reifies s Tape => f (ReverseDouble s) -> ReverseDouble s) -> f Double -> f b-gradWith g f as = reifyTape (snd bds) $ \p -> unbindWith g vs $! partialArrayOf p bds $! f vs+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 #-} @@ -104,8 +113,13 @@ -- @ -- 'grad'' == 'gradWith'' (\_ dx -> dx) -- @-gradWith' :: (Traversable f) => (Double -> Double -> b) -> (forall s. Reifies s Tape => f (ReverseDouble s) -> ReverseDouble s) -> f Double -> (Double, f b)-gradWith' g f as = reifyTape (snd bds) $ \p -> case f vs of+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' #-}@@ -114,8 +128,12 @@ -- -- >>> jacobian (\[x,y] -> [y,x,x*y]) [2,1] -- [[0.0,1.0],[1.0,0.0],[1.0,2.0]]-jacobian :: (Traversable f, Functor g) => (forall s. Reifies s Tape => f (ReverseDouble s) -> g (ReverseDouble s)) -> f Double -> g (f Double)-jacobian f as = reifyTape (snd bds) $ \p -> unbind vs . partialArrayOf p bds <$> f vs where+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 #-} @@ -125,8 +143,12 @@ -- -- >>> jacobian' (\[x,y] -> [y,x,x*y]) [2,1] -- [(1.0,[0.0,1.0]),(2.0,[1.0,0.0]),(2.0,[1.0,2.0])]-jacobian' :: (Traversable f, Functor g) => (forall s. Reifies s Tape => f (ReverseDouble s) -> g (ReverseDouble s)) -> f Double -> g (Double, f Double)-jacobian' f as = reifyTape (snd bds) $ \p ->+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@@ -140,8 +162,13 @@ -- 'jacobian' == 'jacobianWith' (\_ dx -> dx) -- 'jacobianWith' 'const' == (\f x -> 'const' x '<$>' f x) -- @-jacobianWith :: (Traversable f, Functor g) => (Double -> Double -> b) -> (forall s. Reifies s Tape => f (ReverseDouble s) -> g (ReverseDouble s)) -> f Double -> g (f b)-jacobianWith g f as = reifyTape (snd bds) $ \p -> unbindWith g vs . partialArrayOf p bds <$> f vs where+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 #-} @@ -152,8 +179,13 @@ -- -- @'jacobian'' == 'jacobianWith'' (\_ dx -> dx)@ ---jacobianWith' :: (Traversable f, Functor g) => (Double -> Double -> b) -> (forall s. Reifies s Tape => f (ReverseDouble s) -> g (ReverseDouble s)) -> f Double -> g (Double, f b)-jacobianWith' g f as = reifyTape (snd bds) $ \p ->+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@@ -163,8 +195,11 @@ -- -- >>> diff sin 0 -- 1.0-diff :: (forall s. Reifies s Tape => ReverseDouble s -> ReverseDouble s) -> Double -> Double-diff f a = reifyTape 1 $ \p -> derivativeOf p $! f (var a 0)+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.@@ -174,8 +209,11 @@ -- -- >>> diff' exp 0 -- (1.0,1.0)-diff' :: (forall s. Reifies s Tape => ReverseDouble s -> ReverseDouble s) -> Double -> (Double, Double)-diff' f a = reifyTape 1 $ \p -> derivativeOf' p $! f (var a 0)+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.@@ -183,16 +221,24 @@ -- >>> diffF (\a -> [sin a, cos a]) 0 -- [1.0,0.0] ---diffF :: (Functor f) => (forall s. Reifies s Tape => ReverseDouble s -> f (ReverseDouble s)) -> Double -> f Double-diffF f a = reifyTape 1 $ \p -> derivativeOf p <$> f (var a 0)+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 => ReverseDouble s -> f (ReverseDouble s)) -> Double -> f (Double, Double)-diffF' f a = reifyTape 1 $ \p -> derivativeOf' p <$> f (var a 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.@@ -201,7 +247,13 @@ -- -- >>> hessian (\[x,y] -> x*y) [1,2] -- [[0.0,1.0],[1.0,0.0]]-hessian :: (Traversable f) => (forall s s'. (Reifies s R.Tape, Reifies s' Tape) => f (On (R.Reverse s (ReverseDouble s'))) -> (On (R.Reverse s (ReverseDouble s')))) -> f Double -> f (f Double)+hessian+  :: 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 #-} @@ -211,6 +263,12 @@ -- -- >>> hessianF (\[x,y] -> [x*y,x+y,exp x*cos y]) [1,2 :: Double] -- [[[0.0,1.0],[1.0,0.0]],[[0.0,0.0],[0.0,0.0]],[[-1.1312043837568135,-2.4717266720048188],[-2.4717266720048188,1.1312043837568135]]]-hessianF :: (Traversable f, Functor g) => (forall s s'. (Reifies s R.Tape, Reifies s' Tape) => f (On (R.Reverse s (ReverseDouble s'))) -> g (On (R.Reverse s (ReverseDouble s')))) -> f Double -> g (f (f Double))+hessianF+  :: (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 #-}

@@ -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

@@ -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' #-}

@@ -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)

@@ -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 #-}

@@ -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

@@ -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

@@ -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' #-}

@@ -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' #-}

@@ -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')@@ -161,14 +211,22 @@ -- | Compute the gradient of a function using forward mode AD. -- -- 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+  :: (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.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'+  :: (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' #-}@@ -176,7 +234,12 @@ -- | 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.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+  :: (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 #-} @@ -187,17 +250,29 @@ -- -- >>> 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' #-}

@@ -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')@@ -140,14 +201,22 @@ -- | Compute the gradient of a function using forward mode AD. -- -- 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+  :: 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.Reverse.grad'' for @n@ inputs, in exchange for better space utilization.-grad' :: Traversable f => (f ForwardDouble -> ForwardDouble) -> f Double -> (Double, f Double)+grad'+  :: Traversable f+  => (f ForwardDouble -> ForwardDouble)+  -> f Double+  -> (Double, f Double) grad' f as = (primal b, tangent <$> bs)     where         (b, bs) = bind' f as@@ -156,7 +225,12 @@ -- | 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.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+  :: Traversable f+  => (Double -> Double -> b)+  -> (f ForwardDouble -> ForwardDouble)+  -> f Double+  -> f b gradWith g f = bindWith g (tangent . f) {-# INLINE gradWith #-} @@ -167,6 +241,11 @@ -- -- >>> 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' #-}

@@ -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

@@ -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 #-}

@@ -1,222 +1,20 @@ {-# LANGUAGE CPP #-}-{-# LANGUAGE FlexibleInstances #-}-{-# LANGUAGE FunctionalDependencies #-}-{-# LANGUAGE MultiParamTypeClasses #-}-{-# LANGUAGE Rank2Types #-}-{-# LANGUAGE ScopedTypeVariables #-}-{-# 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"

@@ -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"

@@ -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"

@@ -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')

@@ -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'

@@ -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

@@ -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!