ad 4.5.4 → 4.5.5

raw patch · 5 files changed

+124/−20 lines, 5 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• CHANGELOG.markdown +8/−0
• ad.cabal +5/−3
• cbits/tape.c +10/−7
• src/Numeric/AD/Internal/Reverse/Double.hs +1/−1
• tests/Regression.hs +100/−9

@@ -1,3 +1,11 @@+4.5.5 [2024.01.28]+------------------+* `Numeric.AD.Mode.Reverse.Double` now handles IEEE floating-point special+  values (e.g., `NaN` and `Inf`) correctly when `ad` is compiled with `+ffi`.+  Note that this increase in floating-point accuracy may come at a slight+  performance penalty in certain applications. If this negatively impacts your+  application, please mention this at https://github.com/ekmett/ad/issues/106.+ 4.5.4 [2023.02.19] ------------------ * Add a `Num (Scalar (Scalar t))` constraint to `On`'s `Mode` instance, which is

@@ -1,5 +1,5 @@ name:          ad-version:       4.5.4+version:       4.5.5 license:       BSD3 license-File:  LICENSE copyright:     (c) Edward Kmett 2010-2021,@@ -19,7 +19,9 @@              , GHC == 8.8.4              , GHC == 8.10.7              , GHC == 9.0.2-             , GHC == 9.2.2+             , GHC == 9.2.7+             , GHC == 9.4.5+             , GHC == 9.6.2 synopsis:      Automatic Differentiation extra-source-files:   .gitignore@@ -116,7 +118,7 @@     array            >= 0.4     && < 0.6,     base             >= 4.9     && < 5,     comonad          >= 4       && < 6,-    containers       >= 0.5     && < 0.7,+    containers       >= 0.5     && < 0.8,     data-reify       >= 0.6     && < 0.7,     erf              >= 2.0     && < 2.1,     free             >= 4.6.1   && < 6,

@@ -88,20 +88,23 @@     while (--idx >= 0)     {       double v = buffer[idx + pTape->offset];-      if (v == 0.0) continue; +      // TODO: if we do not care about handling IEEE floating point special values (NaN, Inf) correctly+      //       then we can skip the rest of the loop body in case v == 0+      //       see also https://github.com/ekmett/ad/issues/106+       int i = pTape->lnk[idx*2];-      double x = pTape->val[idx*2];-      if (x != 0.0)+      if (i >= 0)       {-        buffer[i] += v*x;+        double x = v * pTape->val[idx*2];+        if (x != 0) buffer[i] += x;       }        int j = pTape->lnk[idx*2 + 1]; -      double y = pTape->val[idx*2 + 1];-      if (y != 0.0)+      if (j >= 0)       {-        buffer[j] += v*y;+        double y = v * pTape->val[idx*2 + 1];+        if (y != 0) buffer[j] += y;       }     }     idx += 1 + pTape->offset;

@@ -115,7 +115,7 @@ -- | 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+unarily f di i b = ReverseDouble (unsafePerformIO (pushTape (Proxy :: Proxy s) i (-1) di 0.0)) $! f b {-# INLINE unarily #-}  -- | This is used to create a new entry on the chain given a binary function, its derivatives with respect to its inputs,

@@ -1,25 +1,116 @@+{-# LANGUAGE NoMonomorphismRestriction #-}+{-# LANGUAGE RankNTypes #-}+ module Main (main) where  import qualified Numeric.AD.Mode.Reverse as R import qualified Numeric.AD.Mode.Reverse.Double as RD +import Text.Printf import Test.Tasty import Test.Tasty.HUnit +type Diff = (forall a. Floating a => a -> a) -> Double -> Double+type Grad = (forall a. Floating a => [a] -> a) -> [Double] -> [Double]+type Jacobian = (forall a. Floating a => [a] -> [a]) -> [Double] -> [[Double]]+type Hessian = (forall a. Floating a => [a] -> a) -> [Double] -> [[Double]]+ main :: IO () main = defaultMain tests  tests :: TestTree-tests = testGroup "Regression tests"-  [ testCase "#97" $-      assertBool "Reverse.diff and Reverse.Double.diff should behave identically" $-      nearZero $ R.diff f (0 :: Double) - RD.diff f (0 :: Double)-  ]+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 :: String -> Diff -> Grad -> Jacobian -> Hessian -> TestTree+mode name diff grad jacobian hessian = testGroup name [basic diff grad jacobian hessian, issue97 diff, issue104 diff grad]++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]+  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]+  tjacobian = testCase "jacobian" $ do+    assertNearMatrix [[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+ -- 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-f :: Num a => a -> a-f = sum . replicate 5000+issue97 :: Diff -> TestTree+issue97 diff = testCase "issue-97" $ assertNear 5000 $ diff f 0 where f = sum . replicate 5000 -nearZero :: (Fractional a, Ord a) => a -> Bool-nearZero a = abs a <= 1e-12+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]+    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]+    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]+    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]+    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)++assertNear :: Double -> Double -> Assertion+assertNear a b = near a b @? expect a b++assertNearList :: [Double] -> [Double] -> Assertion+assertNearList a b = nearList a b @? expect a b++assertNearMatrix :: [[Double]] -> [[Double]] -> Assertion+assertNearMatrix a b = nearMatrix a b @? expect a b++expect :: Show a => a -> a -> String+expect a b = printf "expected %s but got %s" (show a) (show b)++nan :: Double+nan = 0 / 0++inf :: Double+inf = 1 / 0