ad 4.5.5 → 4.5.6
raw patch · 16 files changed
+134/−70 lines, 16 filesdep ~containersPVP ok
version bump matches the API change (PVP)
Dependency ranges changed: containers
API changes (from Hackage documentation)
Files
- CHANGELOG.markdown +5/−0
- ad.cabal +1/−1
- include/instances.h +5/−0
- include/internal_kahn.h +3/−2
- src/Numeric/AD/Internal/Dense.hs +1/−0
- src/Numeric/AD/Internal/Dense/Representable.hs +1/−0
- src/Numeric/AD/Internal/Forward.hs +1/−0
- src/Numeric/AD/Internal/Forward/Double.hs +1/−0
- src/Numeric/AD/Internal/Kahn.hs +3/−2
- src/Numeric/AD/Internal/Reverse.hs +1/−0
- src/Numeric/AD/Internal/Reverse/Double.hs +1/−0
- src/Numeric/AD/Internal/Sparse.hs +1/−0
- src/Numeric/AD/Internal/Sparse/Double.hs +1/−0
- src/Numeric/AD/Internal/Tower.hs +1/−0
- src/Numeric/AD/Internal/Tower/Double.hs +1/−0
- tests/Regression.hs +107/−65
CHANGELOG.markdown view
@@ -1,3 +1,8 @@+4.5.6 [2024.05.01]+------------------+* Add specialized implementations of `log1p`, `expm1`, `log1pexp`, and+ `log1mexp` in `Floating` instances.+ 4.5.5 [2024.01.28] ------------------ * `Numeric.AD.Mode.Reverse.Double` now handles IEEE floating-point special
ad.cabal view
@@ -1,5 +1,5 @@ name: ad-version: 4.5.5+version: 4.5.6 license: BSD3 license-File: LICENSE copyright: (c) Edward Kmett 2010-2021,
include/instances.h view
@@ -59,6 +59,11 @@ acosh = lift1 acosh $ \x -> recip (sqrt (join (*) x - 1)) atanh = lift1 atanh $ \x -> recip (1 - join (*) x) + log1p = lift1 log1p $ recip . (+) 1+ expm1 = lift1 expm1 exp+ log1pexp = lift1 log1pexp $ recip . (+) 1 . exp . negate+ log1mexp = lift1 log1mexp $ recip . negate . expm1 . negate+ instance BODY2(Num a, Enum a) Enum HEAD where succ = lift1 succ (const 1) pred = lift1 pred (const 1)
include/internal_kahn.h view
@@ -53,7 +53,7 @@ import Control.Monad.ST import Control.Monad hiding (mapM) import Control.Monad.Trans.State-import Data.List (foldl')+import qualified Data.List as List (foldl') import Data.Array.ST import Data.Array.IArray import qualified Data.Array as A@@ -67,6 +67,7 @@ import Data.Data (Data) import Data.Typeable (Typeable) import qualified GHC.Exts as Exts+import Numeric import Numeric.AD.Internal.Combinators import Numeric.AD.Internal.Identity import Numeric.AD.Jacobian@@ -234,7 +235,7 @@ backPropagate vmap ss return ss - sbounds ((a,_):as) = foldl' (\(lo,hi) (b,_) -> let lo' = min lo b; hi' = max hi b in lo' `seq` hi' `seq` (lo', hi')) (a,a) as+ sbounds ((a,_):as) = List.foldl' (\(lo,hi) (b,_) -> let lo' = min lo b; hi' = max hi b in lo' `seq` hi' `seq` (lo', hi')) (a,a) as sbounds _ = undefined -- the graph can't be empty, it contains the output node! successors :: Tape Int -> [Int]
src/Numeric/AD/Internal/Dense.hs view
@@ -45,6 +45,7 @@ import Data.Traversable (mapAccumL) import Data.Data () import Data.Number.Erf+import Numeric import Numeric.AD.Internal.Combinators import Numeric.AD.Internal.Identity import Numeric.AD.Jacobian
src/Numeric/AD/Internal/Dense/Representable.hs view
@@ -37,6 +37,7 @@ import Data.Typeable () import Data.Data () import Data.Number.Erf+import Numeric import Numeric.AD.Internal.Combinators import Numeric.AD.Internal.Identity import Numeric.AD.Jacobian
src/Numeric/AD/Internal/Forward.hs view
@@ -41,6 +41,7 @@ import Data.Traversable (mapAccumL) import Data.Data import Data.Number.Erf+import Numeric import Numeric.AD.Internal.Combinators import Numeric.AD.Internal.Identity import Numeric.AD.Jacobian
src/Numeric/AD/Internal/Forward/Double.hs view
@@ -37,6 +37,7 @@ import Data.Traversable (mapAccumL) import Control.Monad (join) import Data.Number.Erf+import Numeric import Numeric.AD.Internal.Combinators import Numeric.AD.Internal.Identity import Numeric.AD.Jacobian
src/Numeric/AD/Internal/Kahn.hs view
@@ -51,7 +51,7 @@ import Control.Monad.ST import Control.Monad hiding (mapM) import Control.Monad.Trans.State-import Data.List (foldl')+import qualified Data.List as List (foldl') import Data.Array.ST import Data.Array import Data.IntMap (IntMap, fromListWith, findWithDefault)@@ -62,6 +62,7 @@ import System.IO.Unsafe (unsafePerformIO) import Data.Data (Data) import Data.Typeable (Typeable)+import Numeric import Numeric.AD.Internal.Combinators import Numeric.AD.Internal.Identity import Numeric.AD.Jacobian@@ -222,7 +223,7 @@ backPropagate vmap ss return ss - sbounds ((a,_):as) = foldl' (\(lo,hi) (b,_) -> let lo' = min lo b; hi' = max hi b in lo' `seq` hi' `seq` (lo', hi')) (a,a) as+ sbounds ((a,_):as) = List.foldl' (\(lo,hi) (b,_) -> let lo' = min lo b; hi' = max hi b in lo' `seq` hi' `seq` (lo', hi')) (a,a) as sbounds _ = undefined -- the graph can't be empty, it contains the output node! successors :: Tape a Int -> [Int]
src/Numeric/AD/Internal/Reverse.hs view
@@ -66,6 +66,7 @@ import Data.Reflection import Data.Traversable (mapM) import Data.Typeable+import Numeric import Numeric.AD.Internal.Combinators import Numeric.AD.Internal.Identity import Numeric.AD.Jacobian
src/Numeric/AD/Internal/Reverse/Double.hs view
@@ -60,6 +60,7 @@ import Data.Reflection import Data.Traversable (mapM) import Data.Typeable+import Numeric import Numeric.AD.Internal.Combinators import Numeric.AD.Internal.Identity import Numeric.AD.Jacobian
src/Numeric/AD/Internal/Sparse.hs view
@@ -52,6 +52,7 @@ import Data.Number.Erf import Data.Traversable import Data.Typeable ()+import Numeric import Numeric.AD.Internal.Combinators import Numeric.AD.Internal.Sparse.Common import Numeric.AD.Jacobian
src/Numeric/AD/Internal/Sparse/Double.hs view
@@ -52,6 +52,7 @@ import Data.Number.Erf import Data.Traversable import Data.Typeable ()+import Numeric import Numeric.AD.Internal.Combinators import Numeric.AD.Internal.Sparse.Common import Numeric.AD.Jacobian
src/Numeric/AD/Internal/Tower.hs view
@@ -40,6 +40,7 @@ import Data.Data (Data) import Data.Number.Erf import Data.Typeable (Typeable)+import Numeric import Numeric.AD.Internal.Combinators import Numeric.AD.Jacobian import Numeric.AD.Mode
src/Numeric/AD/Internal/Tower/Double.hs view
@@ -46,6 +46,7 @@ import Data.Data (Data) import Data.Number.Erf import Data.Typeable (Typeable)+import Numeric import Numeric.AD.Internal.Combinators import Numeric.AD.Jacobian import Numeric.AD.Mode
tests/Regression.hs view
@@ -1,8 +1,12 @@+{-# LANGUAGE CPP #-} {-# LANGUAGE NoMonomorphismRestriction #-} {-# LANGUAGE RankNTypes #-} module Main (main) where +import Numeric+import qualified Numeric.AD.Mode.Forward as F+import qualified Numeric.AD.Mode.Forward.Double as FD import qualified Numeric.AD.Mode.Reverse as R import qualified Numeric.AD.Mode.Reverse.Double as RD @@ -10,7 +14,7 @@ import Test.Tasty import Test.Tasty.HUnit -type Diff = (forall a. Floating a => a -> a) -> Double -> Double+type Diff' = (forall a. Floating a => a -> a) -> Double -> (Double, Double) type Grad = (forall a. Floating a => [a] -> a) -> [Double] -> [Double] type Jacobian = (forall a. Floating a => [a] -> [a]) -> [Double] -> [[Double]] type Hessian = (forall a. Floating a => [a] -> a) -> [Double] -> [[Double]]@@ -18,96 +22,134 @@ main :: IO () main = defaultMain tests +-- TODO: the forward-double tests are currently failing due to discrepancies between the modes+-- see also https://github.com/ekmett/ad/issues/109 and https://github.com/ekmett/ad/pull/110 tests :: TestTree tests = testGroup "tests" [- mode "reverse" (\ f -> R.diff f) (\ f -> R.grad f) (\ f -> R.jacobian f) (\ f -> R.hessian f),- mode "reverse-double" (\ f -> RD.diff f) (\ f -> RD.grad f) (\ f -> RD.jacobian f) (\ f -> RD.hessian f)]+ mode "forward" (\ f -> F.diff' f) (\ f -> F.grad f) (\ f -> F.jacobian f) (\ f -> F.jacobian $ F.grad f),+ --mode "forward-double" (\ f -> FD.diff' f) (\ f -> FD.grad f) (\ f -> FD.jacobian f) (\ f -> FD.jacobian $ F.grad f),+ mode "reverse" (\ f -> R.diff' f) (\ f -> R.grad f) (\ f -> R.jacobian f) (\ f -> R.hessian f),+ mode "reverse-double" (\ f -> RD.diff' f) (\ f -> RD.grad f) (\ f -> RD.jacobian f) (\ f -> RD.hessian f)] -mode :: String -> Diff -> Grad -> Jacobian -> Hessian -> TestTree-mode name diff grad jacobian hessian = testGroup name [basic diff grad jacobian hessian, issue97 diff, issue104 diff grad]+mode :: String -> Diff' -> Grad -> Jacobian -> Hessian -> TestTree+mode name diff grad jacobian hessian = testGroup name [+ basic diff grad jacobian hessian,+ issue97 diff,+ issue104 diff grad,+ issue108 diff] -basic :: Diff -> Grad -> Jacobian -> Hessian -> TestTree+basic :: Diff' -> Grad -> Jacobian -> Hessian -> TestTree basic diff grad jacobian hessian = testGroup "basic" [tdiff, tgrad, tjacobian, thessian] where tdiff = testCase "diff" $ do- assertNearList [11, 5.5, 3, 3.5, 7, 13.5, 23, 35.5, 51] $ diff p <$> [-2, -1.5, -1, -0.5, 0, 0.5, 1, 1.5, 2]- assertNearList [nan, inf, 1, 0.5, 0.25] $ diff sqrt <$> [-1, 0, 0.25, 1, 4]- assertNearList [1, 0, 1] $ [diff sin, diff cos, diff tan] <*> [0]- assertNearList [-1, 0, 1] $ diff abs <$> [-1, 0, 1]- assertNearList [1, exp 1, inf, 1] $ [diff exp, diff log] <*> [0, 1]+ expect (list eq) [11, 5.5, 3, 3.5, 7, 13.5, 23, 35.5, 51] $ snd . diff p <$> [-2, -1.5, -1, -0.5, 0, 0.5, 1, 1.5, 2]+ expect (list eq) [nan, inf, 1, 0.5, 0.25] $ snd . diff sqrt <$> [-1, 0, 0.25, 1, 4]+ expect (list eq) [1, 0, 1] $ [snd . diff sin, snd . diff cos, snd . diff tan] <*> [0]+ expect (list eq) [-1, 0, 1] $ snd . diff abs <$> [-1, 0, 1]+ expect (list eq) [1, exp 1, inf, 1] $ [snd . diff exp, snd . diff log] <*> [0, 1] tgrad = testCase "grad" $ do- assertNearList [2, 1, 1] $ grad f [1, 2, 3]- assertNearList [1, 0.25] $ grad h [2, 8]- assertNearList [0, nan] $ grad power [0, 2]+ expect (list eq) [2, 1, 1] $ grad f [1, 2, 3]+ expect (list eq) [1, 0.25] $ grad h [2, 8]+ expect (list eq) [0, nan] $ grad power [0, 2] tjacobian = testCase "jacobian" $ do- assertNearMatrix [[0, 1], [1, 0], [1, 2]] $ jacobian g [2, 1]+ expect (list $ list eq) [[0, 1], [1, 0], [1, 2]] $ jacobian g [2, 1] thessian = testCase "hessian" $ do- assertNearMatrix [[0, 1, 0], [1, 0, 0], [0, 0, 0]] $ hessian f [1, 2, 3]- assertNearMatrix [[0, 0], [0, 0]] $ hessian sum [1, 2]- assertNearMatrix [[0, 1], [1, 0]] $ hessian product [1, 2]- assertNearMatrix [[2, 1], [1, 0]] $ hessian power [1, 2]- sum = \ [x, y] -> x + y- product = \ [x, y] -> x * y- power = \ [x, y] -> x ** y- f = \ [x, y, z] -> x * y + z- g = \ [x, y] -> [y, x, x * y]- h = \ [x, y] -> sqrt $ x * y- p = \ x -> 12 + 7 * x + 5 * x ^ 2 + 2 * x ^ 3+ expect (list $ list eq) [[0, 1, 0], [1, 0, 0], [0, 0, 0]] $ hessian f [1, 2, 3]+ expect (list $ list eq) [[0, 0], [0, 0]] $ hessian sum [1, 2]+ expect (list $ list eq) [[0, 1], [1, 0]] $ hessian product [1, 2]+ expect (list $ list eq) [[2, 1], [1, 0]] $ hessian power [1, 2]+ sum [x, y] = x + y+ product [x, y] = x * y+ power [x, y] = x ** y+ f [x, y, z] = x * y + z+ g [x, y] = [y, x, x * y]+ h [x, y] = sqrt $ x * y+ p x = 12 + 7 * x + 5 * x ^ 2 + 2 * x ^ 3 -- Reverse.Double +ffi initializes the tape with a block of size 4096 -- The large term in this function forces the allocation of an additional block-issue97 :: Diff -> TestTree-issue97 diff = testCase "issue-97" $ assertNear 5000 $ diff f 0 where f = sum . replicate 5000+issue97 :: Diff' -> TestTree+issue97 diff = testCase "issue-97" $ expect eq 5000 $ snd $ diff f 0 where f = sum . replicate 5000 -issue104 :: Diff -> Grad -> TestTree+issue104 :: Diff' -> Grad -> TestTree issue104 diff grad = testGroup "issue-104" [inside, outside] where inside = testGroup "inside" [tdiff, tgrad] where tdiff = testCase "diff" $ do- assertNearList [nan, nan] $ diff (0 `f`) <$> [0, 1]- assertNearList [inf, 0.5] $ diff (1 `f`) <$> [0, 1]- assertNearList [nan, nan] $ diff (`f` 0) <$> [0, 1]- assertNearList [inf, 0.5] $ diff (`f` 1) <$> [0, 1]+ expect (list eq) [nan, nan] $ snd . diff (0 `f`) <$> [0, 1]+ expect (list eq) [inf, 0.5] $ snd . diff (1 `f`) <$> [0, 1]+ expect (list eq) [nan, nan] $ snd . diff (`f` 0) <$> [0, 1]+ expect (list eq) [inf, 0.5] $ snd . diff (`f` 1) <$> [0, 1] tgrad = testCase "grad" $ do- assertNearList [nan, nan] $ grad (binary f) [0, 0]- assertNearList [nan, inf] $ grad (binary f) [1, 0]- assertNearList [inf, nan] $ grad (binary f) [0, 1]- assertNearList [0.5, 0.5] $ grad (binary f) [1, 1]+ expect (list eq) [nan, nan] $ grad (binary f) [0, 0]+ expect (list eq) [nan, inf] $ grad (binary f) [1, 0]+ expect (list eq) [inf, nan] $ grad (binary f) [0, 1]+ expect (list eq) [0.5, 0.5] $ grad (binary f) [1, 1] f x y = sqrt $ x * y -- grad f [x, y] = [y / (2 * f x y), x / (2 * f x y)] outside = testGroup "outside" [tdiff, tgrad] where tdiff = testCase "diff" $ do- assertNearList [nan, 0.0] $ diff (0 `f`) <$> [0, 1]- assertNearList [inf, 0.5] $ diff (1 `f`) <$> [0, 1]- assertNearList [nan, 0.0] $ diff (`f` 0) <$> [0, 1]- assertNearList [inf, 0.5] $ diff (`f` 1) <$> [0, 1]+ expect (list eq) [nan, 0.0] $ snd . diff (0 `f`) <$> [0, 1]+ expect (list eq) [inf, 0.5] $ snd . diff (1 `f`) <$> [0, 1]+ expect (list eq) [nan, 0.0] $ snd . diff (`f` 0) <$> [0, 1]+ expect (list eq) [inf, 0.5] $ snd . diff (`f` 1) <$> [0, 1] tgrad = testCase "grad" $ do- assertNearList [nan, nan] $ grad (binary f) [0, 0]- assertNearList [0.0, inf] $ grad (binary f) [1, 0]- assertNearList [inf, 0.0] $ grad (binary f) [0, 1]- assertNearList [0.5, 0.5] $ grad (binary f) [1, 1]+ expect (list eq) [nan, nan] $ grad (binary f) [0, 0]+ expect (list eq) [0.0, inf] $ grad (binary f) [1, 0]+ expect (list eq) [inf, 0.0] $ grad (binary f) [0, 1]+ expect (list eq) [0.5, 0.5] $ grad (binary f) [1, 1] f x y = sqrt x * sqrt y -- grad f [x, y] = [sqrt y / 2 sqrt x, sqrt x / 2 sqrt y]- binary f = \ [x, y] -> f x y--near :: Double -> Double -> Bool-near a b = bothNaN || bothInfinite || abs (a - b) <= 1e-12 where- bothNaN = isNaN a && isNaN b- bothInfinite = signum a == signum b && isInfinite a && isInfinite b--nearList :: [Double] -> [Double] -> Bool-nearList as bs = length as == length bs && and (zipWith near as bs)--nearMatrix :: [[Double]] -> [[Double]] -> Bool-nearMatrix as bs = length as == length bs && and (zipWith nearList as bs)+ binary f [x, y] = f x y -assertNear :: Double -> Double -> Assertion-assertNear a b = near a b @? expect a b+issue108 :: Diff' -> TestTree+issue108 diff = testGroup "issue-108" [tlog1p, texpm1, tlog1pexp, tlog1mexp] where+ tlog1p = testCase "log1p" $ do+ equal (-inf, inf) $ diff log1p (-1)+ equal (-1.0000000000000007e-15, 1.000000000000001) $ diff log1p (-1e-15)+ equal (-1e-20, 1) $ diff log1p (-1e-20)+ equal (0, 1) $ diff log1p 0+ equal (1e-20, 1) $ diff log1p 1e-20+ equal (9.999999999999995e-16, 0.9999999999999989) $ diff log1p 1e-15+ equal (0.6931471805599453, 0.5) $ diff log1p 1+ texpm1 = testCase "expm1" $ do+ equal (-0.6321205588285577, 0.36787944117144233) $ diff expm1 (-1)+ equal (-9.999999999999995e-16, 0.999999999999999) $ diff expm1 (-1e-15)+ equal (-1e-20, 1) $ diff expm1 (-1e-20)+ equal (0, 1) $ diff expm1 0+ equal (1e-20, 1) $ diff expm1 1e-20+ equal (1.0000000000000007e-15, 1.000000000000001) $ diff expm1 1e-15+ equal (1.718281828459045, 2.718281828459045) $ diff expm1 1+ tlog1pexp = testCase "log1pexp" $ do+ equal (0, 0) $ diff log1pexp (-1000)+ equal (3.720075976020836e-44, 3.7200759760208356e-44) $ diff log1pexp (-100)+ equal (0.31326168751822286, 0.2689414213699951) $ diff log1pexp (-1)+ equal (0.6931471805599453, 0.5) $ diff log1pexp 0+ equal (1.3132616875182228, 0.7310585786300049) $ diff log1pexp 1+ equal (100, 1) $ diff log1pexp 100+ equal (1000, 1) $ diff log1pexp 1000+ tlog1mexp = testCase "log1mexp" $ do+ equal (-0, -0) $ diff log1mexp (-1000)+-- old versions of base have a faulty implementation of log1mexp, causing this case to fail+-- see also https://gitlab.haskell.org/ghc/ghc/-/issues/17125+#if MIN_VERSION_base(4, 13, 0)+ equal (-3.720075976020836e-44, -3.7200759760208356e-44) $ diff log1mexp (-100)+#endif+ equal (-0.45867514538708193, -0.5819767068693265) $ diff log1mexp (-1)+ equal (-0.9327521295671886, -1.5414940825367982) $ diff log1mexp (-0.5)+ equal (-2.3521684610440907, -9.50833194477505) $ diff log1mexp (-0.1)+ equal (-34.538776394910684, -9.999999999999994e14) $ diff log1mexp (-1e-15)+ equal (-46.051701859880914, -1e20) $ diff log1mexp (-1e-20)+ equal (-inf, -inf) $ diff log1mexp (-0)+ equal = expect $ \ (a, b) (c, d) -> eq a c && eq b d -assertNearList :: [Double] -> [Double] -> Assertion-assertNearList a b = nearList a b @? expect a b+-- TODO: ideally, we would consider `0` and `-0` to be different+-- however, zero signedness is currently not reliably propagated through some modes+-- see also https://github.com/ekmett/ad/issues/109 and https://github.com/ekmett/ad/pull/110+eq :: Double -> Double -> Bool+eq a b = isNaN a && isNaN b || a == b -assertNearMatrix :: [[Double]] -> [[Double]] -> Assertion-assertNearMatrix a b = nearMatrix a b @? expect a b+list :: (a -> a -> Bool) -> [a] -> [a] -> Bool+list eq as bs = length as == length bs && and (zipWith eq as bs) -expect :: Show a => a -> a -> String-expect a b = printf "expected %s but got %s" (show a) (show b)+expect :: HasCallStack => Show a => (a -> a -> Bool) -> a -> a -> Assertion+expect eq a b = eq a b @? printf "expected %s but got %s" (show a) (show b) nan :: Double nan = 0 / 0