ad 4.3 → 4.3.1

raw patch · 3 files changed

+50/−73 lines, 3 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• CHANGELOG.markdown +6/−0
• ad.cabal +2/−2
• src/Numeric/AD/Internal/Sparse.hs +42/−71

@@ -1,3 +1,9 @@+4.3.1+-----+* Further improvements have been made in the performance of `Sparse` mode, at least asymptotically, when used on functions with many variables.+  Since this is the target use-case for `Sparse` in the first place, this seems like a good trade-off. Note: this results in an API change, but+  only in the API of an `Internal` module, so this is treated as a minor version bump.+ 4.3 --- * Made drastic improvements in the performance of `Tower` and `Sparse` modes thanks to the help of Björn von Sydow.

@@ -1,5 +1,5 @@ name:          ad-version:       4.3+version:       4.3.1 license:       BSD3 license-File:  LICENSE copyright:     (c) Edward Kmett 2010-2015,@@ -12,7 +12,7 @@ bug-reports:   http://github.com/ekmett/ad/issues build-type:    Custom cabal-version: >= 1.10-tested-with:   GHC==7.0.1, GHC == 7.0.4, GHC == 7.2.2, GHC == 7.4.2, GHC == 7.6.3, GHC == 7.8.4, GHC == 7.10.1+tested-with:   GHC==7.0.1, GHC == 7.0.4, GHC == 7.2.2, GHC == 7.4.2, GHC == 7.6.3, GHC == 7.8.4, GHC == 7.10.1, GHC == 7.10.2 synopsis:      Automatic Differentiation extra-source-files:   .ghci

@@ -23,9 +23,9 @@ -- Handle with care. ----------------------------------------------------------------------------- module Numeric.AD.Internal.Sparse-  ( Index(..)-  , emptyIndex-  , addToIndex+  ( Monomial(..)+  , emptyMonomial+  , addToMonomial   , indices   , Sparse(..)   , apply@@ -40,7 +40,6 @@   , Grad(..)   , Grads(..)   , terms-  , deriv   , primal   ) where @@ -60,18 +59,18 @@ import Numeric.AD.Jacobian import Numeric.AD.Mode -newtype Index = Index (IntMap Int)+newtype Monomial = Monomial (IntMap Int) -emptyIndex :: Index-emptyIndex = Index IntMap.empty-{-# INLINE emptyIndex #-}+emptyMonomial :: Monomial+emptyMonomial = Monomial IntMap.empty+{-# INLINE emptyMonomial #-} -addToIndex :: Int -> Index -> Index-addToIndex k (Index m) = Index (insertWith (+) k 1 m)-{-# INLINE addToIndex #-}+addToMonomial :: Int -> Monomial -> Monomial+addToMonomial k (Monomial m) = Monomial (insertWith (+) k 1 m)+{-# INLINE addToMonomial #-} -indices :: Index -> [Int]-indices (Index as) = uncurry (flip replicate) `concatMap` toAscList as+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@@ -84,24 +83,6 @@   | Zero   deriving (Show, Data, Typeable) -{---These functions are now unused.--dropMap :: Int -> IntMap a -> IntMap a-dropMap n = snd . IntMap.split (n - 1)-{-# INLINE dropMap #-}--times :: Num a => Sparse a -> Int -> Sparse a -> Sparse a-times Zero _ _ = Zero-times _ _ Zero = Zero-times a@(Sparse pa da) n b@(Sparse pb db) = Sparse (pa * pb) $-  unionWith (+)-    (fmap (* b) (dropMap n da))-    (fmap (a *) (dropMap n db))-{-# INLINE times #-}--}- vars :: (Traversable f, Num a) => f a -> f (Sparse a) vars = snd . mapAccumL var 0 where   var !n a = (n + 1, Sparse a $ singleton n $ auto 1)@@ -127,13 +108,19 @@  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 emptyIndex <$> fns) where+ds fs (as@(Sparse a _)) = a :< (go emptyMonomial <$> fns) where   fns = skeleton fs-  -- go :: Index -> Int -> Cofree f a+  -- go :: Monomial -> Int -> Cofree f a   go ix i = partial (indices ix') as :< (go ix' <$> fns) where-    ix' = addToIndex i ix+    ix' = addToMonomial i ix {-# INLINE ds #-} +partialS :: Num a => [Int] -> Sparse a -> Sparse a+partialS []     a             = a+partialS (n:ns) (Sparse _ da) = partialS ns $ findWithDefault Zero n da+partialS _      Zero          = Zero+{-# INLINE partialS #-}+ partial :: Num a => [Int] -> Sparse a -> a partial []     (Sparse a _)  = a partial (n:ns) (Sparse _ da) = partial ns $ findWithDefault (auto 0) n da@@ -265,56 +252,40 @@ isZero _ = False  -- |--- A monomial is used to indicate order of differentiation.--- For a k-ary function, it represented as a list of k non-negative Ints.--- MI [n_0,n_1...n_{k-1}] denotes differentiation n_0 times with respect--- to variable 0, n_1 times to variable 1, etc.--- Trailing zeros omitted for efficiency.------ Add 1 to variable k (i.e.differentiate once more wrt variable k).-incMonomial :: Int -> [Int] -> [Int]-incMonomial k [] = replicate k 0 ++ [1]-incMonomial 0 (a:as) = a+1:as-incMonomial k (a:as) = a:incMonomial (k-1) as---- deriv f mi is the derivative of f of order mi (including higher derivatives).-deriv :: Sparse a -> [Int] -> Sparse a-deriv f mi = indx 0 mi f where-  indx _ [] f = f-  indx _ _ Zero = Zero-  indx v (0:as) f = indx (v+1) as f-  indx v (a:as) (Sparse _ df) = maybe Zero (indx v (a-1 : as)) (lookup v df)- -- The value of the derivative of (f*g) of order mi is---       sum [a*primal (deriv f b)*primal (deriv g c) | (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 Zeros 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+--+-- @+-- '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 :: [Int]-> [(Integer,[Int],[Int])]-terms [] = [(1,[],[])]-terms (a:as) = concatMap (f ps) (zip (bins!!a) [0..a]) where-  ps = terms as+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"-  f ps (b,k) = map (\(w,ks,is) -> (w*b,(k:ks),(a-k:is))) ps  mul :: Num a => Sparse a -> Sparse a -> Sparse a mul Zero _ = Zero mul _ Zero = Zero-mul f@(Sparse _ am) g@(Sparse _ bm) = Sparse (primal f * primal g) (derivs 0 []) where+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' = incMonomial w mi+        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 = deriv f mif-          gder = deriv g mig-  kMax = max (maximum (-1:IntMap.keys am)) (maximum (-1:IntMap.keys bm))+          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)