ad 1.3.1 → 1.4

raw patch · 9 files changed

+247/−80 lines, 9 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

• Numeric/AD/Internal/Classes.hs +23/−8
• Numeric/AD/Internal/Composition.hs +2/−1
• Numeric/AD/Internal/Dense.hs +65/−35
• Numeric/AD/Internal/Forward.hs +80/−11
• Numeric/AD/Internal/Identity.hs +2/−1
• Numeric/AD/Internal/Reverse.hs +15/−1
• Numeric/AD/Internal/Sparse.hs +52/−19
• Numeric/AD/Internal/Tower.hs +4/−2
• ad.cabal +4/−2

@@ -1,4 +1,5 @@-{-# LANGUAGE Rank2Types, TypeFamilies, FlexibleInstances, MultiParamTypeClasses, FlexibleContexts, FunctionalDependencies, UndecidableInstances, GeneralizedNewtypeDeriving, TemplateHaskell, CPP #-}+{-# LANGUAGE Rank2Types, TypeFamilies, FlexibleInstances, MultiParamTypeClasses, PatternGuards, CPP #-}+{-# LANGUAGE FlexibleContexts, FunctionalDependencies, UndecidableInstances, GeneralizedNewtypeDeriving, TemplateHaskell #-} -- {-# OPTIONS_HADDOCK hide #-} ----------------------------------------------------------------------------- -- |@@ -25,12 +26,12 @@     , Iso(..)     ) where -import Control.Applicative+import Control.Applicative hiding ((<**>)) import Data.Char import Language.Haskell.TH import Numeric.AD.Internal.Combinators (on) -infixl 8 **!+infixr 8 **!, <**> infixl 7 *!, /!, ^*, *^, ^/ infixl 6 +!, -!, <+> infix 4 ==!@@ -82,7 +83,14 @@     maxBound1       :: (Num a, Bounded a) => t a  class Lifted t => Mode t where+    -- | allowed to return False for items with a zero derivative, but we'll give more NaNs than strictly necessary+    isKnownConstant :: t a -> Bool+    isKnownConstant _ = False +    -- | allowed to return False for zero, but we give more NaN's than strictly necessary then+    isKnownZero :: Num a => t a -> Bool+    isKnownZero _ = False+     -- | Embed a constant     lift  :: Num a => a -> t a @@ -98,6 +106,10 @@     -- | Scalar division     (^/) :: Fractional a => t a -> a -> t a +    -- | Exponentiation, this should be overloaded if you can figure out anything about what is constant!+    (<**>) :: Floating a => t a -> t a -> t a+--  x <**> y = lift2_ (**) (\z xi yi -> (yi *! z /! xi, z *! log1 xi)) x y+     -- | > 'zero' = 'lift' 0     zero :: Num a => t a @@ -183,14 +195,14 @@ deriveLifted f _t = do         [InstanceD cxt0 type0 dec0] <- lifted         return <$> instanceD (cxt (f (return <$> cxt0))) (return type0) (return <$> dec0)-    where -      lifted = [d| +    where+      lifted = [d|        instance Lifted $_t where         (==!)         = (==) `on` primal         compare1      = compare `on` primal         maxBound1     = lift maxBound         minBound1     = lift minBound-        showsPrec1 d  = showsPrec d . primal +        showsPrec1 d  = showsPrec d . primal         fromInteger1  = lift . fromInteger         (+!)          = (<+>) -- binary (+) one one         (-!)          = binary (-) one negOne -- TODO: <-> ? as it is, this might be pretty bad for Tower@@ -201,13 +213,16 @@         fromRational1 = lift . fromRational         x /! y        = x *! recip1 y         recip1        = lift1_ recip (const . negate1 . square1)-         pi1       = lift pi         exp1      = lift1_ exp const         log1      = lift1 log recip1         logBase1 x y = log1 y /! log1 x         sqrt1     = lift1_ sqrt (\z _ -> recip1 (lift 2 *! z))-        (**!)     = lift2_ (**) (\z x y -> (y *! z /! x, z *! log1 x)) -- error at 0 ** n+        (**!)     = (<**>)+        --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         sin1      = lift1 sin cos1         cos1      = lift1 cos $ negate1 . sin1         tan1 x    = sin1 x /! cos1 x

@@ -18,7 +18,7 @@     , decomposeMode     ) where -import Control.Applicative+import Control.Applicative hiding ((<**>)) import Data.Data (Data(..), mkDataType, DataType, mkConstr, Constr, constrIndex, Fixity(..)) import Data.Typeable (Typeable1(..), Typeable(..), TyCon, mkTyCon3, mkTyConApp, typeOfDefault, gcast1) import Data.Foldable (Foldable(foldMap))@@ -84,6 +84,7 @@     a *^ ComposeMode b = ComposeMode (lift a *^ b)     ComposeMode a ^* b = ComposeMode (a ^* lift b)     ComposeMode a ^/ b = ComposeMode (a ^/ lift b)+    ComposeMode a <**> ComposeMode b = ComposeMode (a <**> b)  instance (Mode f, Mode g) => Lifted (ComposeMode f g) where     showsPrec1 n (ComposeMode a) = showsPrec1 n a

@@ -44,19 +44,22 @@ data Dense f a     = Lift !a     | Dense !a (f a)+    | Zero  instance Show a => Show (Dense f a) where-    showsPrec n (Lift a) = showsPrec n a-    showsPrec n (Dense a _) = showsPrec n a+    showsPrec d (Lift a)    = showsPrec d a+    showsPrec d (Dense a _) = showsPrec d a+    showsPrec _ Zero        = showString "0"  ds :: f a -> AD (Dense f) a -> f a-ds _    (AD (Dense _ da)) = da-ds z _                = z+ds _ (AD (Dense _ da)) = da+ds z _ = z {-# INLINE ds #-} -ds' :: f a -> AD (Dense f) a -> (a, f a)-ds' _    (AD (Dense a da)) = (a, da)-ds' z (AD (Lift a))    = (a, z)+ds' :: Num a => f a -> AD (Dense f) a -> (a, f a)+ds' _ (AD (Dense a da)) = (a, da)+ds' z (AD (Lift a)) = (a, z)+ds' z (AD Zero) = (0, z) {-# INLINE ds' #-}  -- Bind variables and count inputs@@ -72,75 +75,102 @@ {-# INLINE apply #-}  instance Primal (Dense f) where+    primal Zero = 0     primal (Lift a) = a     primal (Dense a _) = a  instance (Traversable f, Lifted (Dense f)) => Mode (Dense f) where     lift = Lift-    Lift a <+> Lift b = Lift (a + b)-    Lift a <+> Dense b db = Dense (a + b) db-    Dense a da <+> Lift b = Dense (a + b) da+    zero = Zero++    Zero <+> a = a+    a <+> Zero = a+    Lift a     <+> Lift b     = Lift (a + b)+    Lift a     <+> Dense b db = Dense (a + b) db+    Dense a da <+> Lift b     = Dense (a + b) da     Dense a da <+> Dense b db = Dense (a + b) $ zipWithT (+) da db-    a *^ Lift b    = Lift (a * b)++    _ <**> Zero   = lift 1+    x <**> Lift y = lift1 (**y) (\z -> (y *^ z ** Id (y-1))) x+    x <**> y      = lift2_ (**) (\z xi yi -> (yi *! z /! xi, z *! log1 xi)) x y++    _ *^ Zero       = Zero+    a *^ Lift b     = Lift (a * b)     a *^ Dense b db = Dense (a * b) $ fmap (a*) db-    Lift a    ^* b = Lift (a * b)+    Zero       ^* _ = Zero+    Lift a     ^* b = Lift (a * b)     Dense a da ^* b = Dense (a * b) $ fmap (*b) da-    Lift a    ^/ b = Lift (a / b)+    Zero       ^/ _ = Zero+    Lift a     ^/ b = Lift (a / b)     Dense a da ^/ b = Dense (a / b) $ fmap (/b) da  instance (Traversable f, Lifted (Dense f)) => Jacobian (Dense f) where     type D (Dense f) = Id+    unary f _         Zero        = Lift (f 0)     unary f _         (Lift b)    = Lift (f b)     unary f (Id dadb) (Dense b db) = Dense (f b) (fmap (dadb *) db) +    lift1 f _  Zero        = Lift (f 0)     lift1 f _  (Lift b)    = Lift (f b)     lift1 f df (Dense b db) = Dense (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 (Dense b db) = Dense a (fmap (dadb *) db)         where             a = f b             Id dadb = df (Id a) (Id b) -    binary f _ _ (Lift b) (Lift c)-        = Lift (f b c)-    binary f _ (Id dadc) (Lift b) (Dense c dc)-        = Dense (f b c) $ fmap (* dadc) dc-    binary f (Id dadb) _ (Dense b db) (Lift c)-        = Dense (f b c) $ fmap (dadb *) db-    binary f (Id dadb) (Id dadc) (Dense b db) (Dense c dc)-        = Dense (f b c) $ zipWithT productRule db dc-        where-            productRule dbi dci = dadb * dbi + dci * dadc+    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         (Dense c dc) = Dense (f 0 c) $ fmap (* dadc) dc+    binary f _         (Id dadc) (Lift b)     (Dense c dc) = Dense (f b c) $ fmap (* dadc) dc+    binary f (Id dadb) _         (Dense b db) Zero         = Dense (f b 0) $ fmap (dadb *) db+    binary f (Id dadb) _         (Dense b db) (Lift c)     = Dense (f b c) $ fmap (dadb *) db+    binary f (Id dadb) (Id dadc) (Dense b db) (Dense c dc) = Dense (f b c) $ zipWithT productRule db dc+        where productRule dbi dci = dadb * dbi + dci * dadc -    lift2 f _  (Lift b) (Lift c)-        = Lift (f b c)-    lift2 f df (Lift b) (Dense c dc)-        = Dense (f b c) $ fmap (*dadc) dc-        where-            (_, Id dadc) = df (Id b) (Id c)-    lift2 f df (Dense b db) (Lift c)-        = Dense (f b c) $ fmap (dadb*) db-        where-            (Id dadb, _) = df (Id b) (Id c)+    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         (Dense c dc) = Dense (f 0 c) $ fmap (*dadc) dc where dadc = runId (snd (df (Id 0) (Id c)))+    lift2 f df (Lift b)     (Dense c dc) = Dense (f b c) $ fmap (*dadc) dc where dadc = runId (snd (df (Id b) (Id c)))+    lift2 f df (Dense b db) Zero         = Dense (f b 0) $ fmap (dadb*) db where dadb = runId (fst (df (Id b) (Id 0)))+    lift2 f df (Dense b db) (Lift c)     = Dense (f b c) $ fmap (dadb*) db where dadb = runId (fst (df (Id b) (Id c)))     lift2 f df (Dense b db) (Dense c dc) = Dense (f b c) da         where             (Id dadb, Id dadc) = df (Id b) (Id c)             da = zipWithT 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     (Dense c dc)+        = Dense a $ fmap (*dadc) dc+        where+            a = f 0 c+            (_, Id dadc) = df (Id a) (Id 0) (Id c)     lift2_ f df (Lift b) (Dense c dc)         = Dense a $ fmap (*dadc) dc         where-            a = (f b c)+            a = f b c             (_, Id dadc) = df (Id a) (Id b) (Id c)+    lift2_ f df (Dense b db) Zero+        = Dense a $ fmap (dadb*) db+        where+            a = f b 0+            (Id dadb, _) = df (Id a) (Id b) (Id 0)     lift2_ f df (Dense b db) (Lift c)         = Dense a $ fmap (dadb*) db         where-            a = (f b c)+            a = f b c             (Id dadb, _) = df (Id a) (Id b) (Id c)     lift2_ f df (Dense b db) (Dense c dc)         = Dense a $ zipWithT productRule db dc

@@ -37,14 +37,21 @@ import Numeric.AD.Internal.Classes import Numeric.AD.Internal.Identity -data Forward a = Forward !a a deriving (Show, Data, Typeable)+data Forward a+  = Forward !a a+  | Lift !a+  | Zero+  deriving (Show, Data, Typeable) -tangent :: AD Forward a -> a+tangent :: Num a => AD Forward a -> a tangent (AD (Forward _ da)) = da+tangent _ = 0 {-# INLINE tangent #-} -unbundle :: AD Forward a -> (a, a)+unbundle :: Num a => AD Forward a -> (a, a) unbundle (AD (Forward a da)) = (a, da)+unbundle (AD Zero) = (0,0)+unbundle (AD (Lift a)) = (a, 0) {-# INLINE unbundle #-}  bundle :: a -> a -> AD Forward a@@ -57,33 +64,95 @@  instance Primal Forward where     primal (Forward a _) = a+    primal (Lift a) = a+    primal Zero = 0  instance Lifted Forward => Mode Forward where-    lift a = Forward a 0+    lift = Lift+    zero = Zero++    isKnownZero Zero = True+    isKnownZero _    = False++    isKnownConstant Forward{} = False+    isKnownConstant _ = True++    Zero <+> a = a+    a <+> Zero = a     Forward a da <+> Forward b db = Forward (a + b) (da + db)+    Forward a da <+> Lift b = Forward (a + b) da+    Lift a <+> Forward b db = Forward (a + b) db+    Lift a <+> Lift b = Lift (a + b)++    _ <**> Zero = lift 1+    x <**> Lift y = lift1 (**y) (\z -> (y *^ z ** Id (y-1))) x+    x <**> y = lift2_ (**) (\z xi yi -> (yi *! z /! xi, z *! log1 xi)) x y+     a *^ Forward b db = Forward (a * b) (a * db)+    a *^ Lift b = Lift (a * b)+    _ *^ Zero = Zero+     Forward a da ^* b = Forward (a * b) (da * b)+    Lift a ^* b = Lift (a * b)+    Zero ^* _ = Zero+     Forward a da ^/ b = Forward (a / b) (da / b)+    Lift a ^/ b = Lift (a / b)+    Zero ^/ _ = Zero  instance Lifted Forward => Jacobian Forward where     type D Forward = Id++     unary f (Id dadb) (Forward b db) = Forward (f b) (dadb * db)+    unary f _         (Lift b)       = Lift (f b)+    unary f _         Zero           = Lift (f 0)++    lift1 f _ Zero            = Lift (f 0)+    lift1 f _  (Lift b)       = Lift (f b)     lift1 f df (Forward b db) = Forward (f b) (dadb * db)         where             Id dadb = df (Id b)++    lift1_ f _  Zero           = Lift (f 0)+    lift1_ f _  (Lift b)       = Lift (f b)     lift1_ f df (Forward b db) = Forward a da         where             a = f b             Id da = df (Id a) (Id b) ^* db -    binary f (Id dadb) (Id dadc) (Forward b db) (Forward c dc) = Forward (f b c) da-        where-            da = dadb * db + dc * dadc+    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           (Forward c dc) = Forward (f 0 c) $ dc * dadc+    binary f _         (Id dadc) (Lift b)       (Forward c dc) = Forward (f b c) $ dc * dadc+    binary f (Id dadb) _         (Forward b db) Zero           = Forward (f b 0) $ dadb * db+    binary f (Id dadb) _         (Forward b db) (Lift c)       = Forward (f b c) $ dadb * db+    binary f (Id dadb) (Id dadc) (Forward b db) (Forward c dc) = Forward (f b c) $ dadb * db + dc * 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           (Forward c dc) = Forward (f 0 c) $ dc * runId (snd (df (Id 0) (Id c)))+    lift2 f df (Lift b)       (Forward c dc) = Forward (f b c) $ dc * runId (snd (df (Id b) (Id c)))+    lift2 f df (Forward b db) Zero           = Forward (f b 0) $ runId (fst (df (Id b) (Id 0))) * db+    lift2 f df (Forward b db) (Lift c)       = Forward (f b c) $ runId (fst (df (Id b) (Id c))) * db     lift2 f df (Forward b db) (Forward c dc) = Forward a da         where             a = f b c             (Id dadb, Id dadc) = df (Id b) (Id c)             da = dadb * db + dc * 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           (Forward c dc) = Forward a $ dc * runId (snd (df (Id a) (Id 0) (Id c))) where a = f 0 c+    lift2_ f df (Lift b)       (Forward c dc) = Forward a $ dc * runId (snd (df (Id a) (Id b) (Id c))) where a = f b c+    lift2_ f df (Forward b db) Zero           = Forward a $ runId (fst (df (Id a) (Id b) (Id 0))) * db where a = f b 0+    lift2_ f df (Forward b db) (Lift c)       = Forward a $ runId (fst (df (Id a) (Id b) (Id c))) * db where a = f b c     lift2_ f df (Forward b db) (Forward c dc) = Forward a da         where             a = f b c@@ -96,13 +165,13 @@ bind f as = snd $ mapAccumL outer (0 :: Int) as     where         outer !i _ = (i + 1, f $ snd $ mapAccumL (inner i) 0 as)-        inner !i !j a = (j + 1, bundle a $ if i == j then 1 else 0)+        inner !i !j a = (j + 1, if i == j then bundle a 1 else AD Zero)  bind' :: (Traversable f, Num a) => (f (AD Forward a) -> b) -> f a -> (b, f b) bind' f as = dropIx $ mapAccumL outer (0 :: Int, b0) as     where         outer (!i, _) _ = let b = f $ snd $ mapAccumL (inner i) (0 :: Int) as in ((i + 1, b), b)-        inner !i !j a = (j + 1, bundle a $ if i == j then 1 else 0)+        inner !i !j a = (j + 1, if i == j then bundle a 1 else AD Zero)         b0 = f (lift <$> as)         dropIx ((_,b),bs) = (b,bs) @@ -110,13 +179,13 @@ bindWith g f as = snd $ mapAccumL outer (0 :: Int) as     where         outer !i a = (i + 1, g a $ f $ snd $ mapAccumL (inner i) 0 as)-        inner !i !j a = (j + 1, bundle a $ if i == j then 1 else 0)+        inner !i !j a = (j + 1, if i == j then bundle a 1 else AD Zero)  bindWith' :: (Traversable f, Num a) => (a -> b -> c) -> (f (AD Forward a) -> b) -> f a -> (b, f c) bindWith' g f as = dropIx $ mapAccumL outer (0 :: Int, b0) as     where         outer (!i, _) a = let b = f $ snd $ mapAccumL (inner i) (0 :: Int) as in ((i + 1, b), g a b)-        inner !i !j a = (j + 1, bundle a $ if i == j then 1 else 0)+        inner !i !j a = (j + 1, if i == j then bundle a 1 else AD Zero)         b0 = f (lift <$> as)         dropIx ((_,b),bs) = (b,bs) 

@@ -27,7 +27,7 @@ import Data.Traversable (Traversable, traverse) import Data.Foldable (Foldable, foldMap) -newtype Id a = Id a deriving+newtype Id a = Id { runId :: a } deriving     (Iso a, Eq, Ord, Show, Enum, Bounded, Num, Real, Fractional, Floating, RealFrac, RealFloat, Monoid, Data, Typeable)  probe :: a -> AD Id a@@ -133,6 +133,7 @@     Id a ^* b = Id (a * b)     a *^ Id b = Id (a * b)     Id a <+> Id b = Id (a + b)+    Id a <**> Id b = Id (a ** b)  instance Primal Id where     primal (Id a) = a

@@ -60,7 +60,8 @@  -- | A @Tape@ records the information needed back propagate from the output to each input during 'Reverse' 'Mode' AD. data Tape a t-    = Lift !a+    = Zero+    | Lift !a     | Var !a {-# UNPACK #-} !Int     | Binary !a a a t t     | Unary !a a t@@ -74,6 +75,7 @@ instance MuRef (Reverse a) where     type DeRef (Reverse a) = Tape a +    mapDeRef _ (Reverse Zero) = pure Zero     mapDeRef _ (Reverse (Lift a)) = pure (Lift a)     mapDeRef _ (Reverse (Var a v)) = pure (Var a v)     mapDeRef f (Reverse (Binary a dadb dadc b c)) = Binary a dadb dadc <$> f b <*> f c@@ -81,12 +83,18 @@  instance Lifted Reverse => Mode Reverse where     lift a = Reverse (Lift a)+    zero   = Reverse Zero     (<+>)  = binary (+) one one     a *^ b = lift1 (a *) (\_ -> lift a) b     a ^* b = lift1 (* b) (\_ -> lift b) a     a ^/ b = lift1 (/ b) (\_ -> lift (recip b)) a +    _ <**> Reverse Zero     = lift 1+    x <**> Reverse (Lift y) = lift1 (**y) (\z -> (y *^ z ** Id (y-1))) x+    x <**> y                = lift2_ (**) (\z xi yi -> (yi *! z /! xi, z *! log1 xi)) x y+ instance Primal Reverse where+    primal (Reverse Zero) = 0     primal (Reverse (Lift a)) = a     primal (Reverse (Var a _)) = a     primal (Reverse (Binary a _ _ _ _)) = a@@ -95,6 +103,7 @@ instance Lifted Reverse => Jacobian Reverse where     type D Reverse = Id +    unary f _         (Reverse Zero)     = Reverse (Lift (f 0))     unary f _         (Reverse (Lift a)) = Reverse (Lift (f a))     unary f (Id dadb) b                  = Reverse (Unary (f (primal b)) dadb b) @@ -105,8 +114,13 @@         where pb = primal b               a = f pb +    binary f _         _         (Reverse Zero)     (Reverse Zero)     = Reverse (Lift (f 0 0))+    binary f _         _         (Reverse Zero)     (Reverse (Lift c)) = Reverse (Lift (f 0 c))+    binary f _         _         (Reverse (Lift b)) (Reverse Zero)     = Reverse (Lift (f b 0))     binary f _         _         (Reverse (Lift b)) (Reverse (Lift c)) = Reverse (Lift (f b c))+    binary f _         (Id dadc) (Reverse Zero)     c                  = Reverse (Unary (f 0 (primal c)) dadc c)     binary f _         (Id dadc) (Reverse (Lift b)) c                  = Reverse (Unary (f b (primal c)) dadc c)+    binary f (Id dadb) _         b                  (Reverse Zero)     = Reverse (Unary (f (primal b) 0) dadb b)     binary f (Id dadb) _         b                  (Reverse (Lift c)) = Reverse (Unary (f (primal b) c) dadb b)     binary f (Id dadb) (Id dadc) b                  c                  = Reverse (Binary (f (primal b) (primal c)) dadb dadc b c) 

@@ -1,6 +1,6 @@ {-# LANGUAGE BangPatterns, TemplateHaskell, TypeFamilies, TypeOperators, FlexibleContexts, UndecidableInstances, DeriveDataTypeable, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances #-} {-# OPTIONS_GHC -fno-warn-name-shadowing #-}-module Numeric.AD.Internal.Sparse +module Numeric.AD.Internal.Sparse     ( Index(..)     , emptyIndex     , addToIndex@@ -20,13 +20,13 @@     ) where  import Prelude hiding (lookup)-import Control.Applicative+import Control.Applicative hiding ((<**>)) import Numeric.AD.Internal.Classes import Control.Comonad.Cofree import Numeric.AD.Internal.Types import Data.Data import Data.Typeable ()-import qualified Data.IntMap as IntMap +import qualified Data.IntMap as IntMap import Data.IntMap (IntMap, mapWithKey, unionWith, findWithDefault, toAscList, singleton, insertWith, lookup) import Data.Traversable import Language.Haskell.TH@@ -48,30 +48,35 @@ -- | 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 +-- there are only (n + k - 1) choose k distinct nth partial derivatives of a -- function with k inputs.-data Sparse a = Sparse !a (IntMap (Sparse a)) deriving (Show, Data, Typeable)+data Sparse a+  = Sparse !a (IntMap (Sparse a))+  | Zero+  deriving (Show, Data, Typeable)  -- | drop keys below a given value dropMap :: Int -> IntMap a -> IntMap a-dropMap n = snd . IntMap.split (n - 1) +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 (Sparse a as) n (Sparse b bs) = Sparse (a * b) $-    unionWith (<+>) +    unionWith (<+>)         (fmap (^* b) (dropMap n as))         (fmap (a *^) (dropMap n bs)) {-# INLINE times #-}  vars :: (Traversable f, Num a) => f a -> f (AD Sparse a)-vars = snd . mapAccumL var 0 +vars = snd . mapAccumL var 0     where         var !n a = (n + 1, AD $ Sparse a $ singleton n $ lift 1) {-# INLINE vars #-}  apply :: (Traversable f, Num a) => (f (AD Sparse a) -> b) -> f a -> b-apply f = f . vars +apply f = f . vars {-# INLINE apply #-}  skeleton :: Traversable f => f a -> f Int@@ -79,14 +84,17 @@ {-# INLINE skeleton #-}  d :: (Traversable f, Num a) => f b -> AD Sparse a -> f a+d fs (AD Zero) = 0 <$ fs d fs (AD (Sparse _ da)) = snd $ mapAccumL (\ !n _ -> (n + 1, maybe 0 primal $ lookup n da)) 0 fs {-# INLINE d #-}  d' :: (Traversable f, Num a) => f a -> AD Sparse a -> (a, f a)-d' fs (AD (Sparse a da)) = (a , snd $ mapAccumL (\ !n _ -> (n + 1, maybe 0 primal $ lookup n da)) 0 fs)+d' fs (AD Zero) = (0, 0 <$ fs)+d' fs (AD (Sparse a da)) = (a, snd $ mapAccumL (\ !n _ -> (n + 1, maybe 0 primal $ lookup n da)) 0 fs) {-# INLINE d' #-}  ds :: (Traversable f, Num a) => f b -> AD Sparse a -> Cofree f a+ds fs (AD Zero) = r where r = 0 :< (r <$ fs) ds fs (AD as@(Sparse a _)) = a :< (go emptyIndex <$> fns)     where         fns = skeleton fs@@ -101,7 +109,7 @@ {-# INLINE vvars #-}  vapply :: Num a => (Vector (AD Sparse a) -> b) -> Vector a -> b-vapply f = f . vvars +vapply f = f . vvars {-# INLINE vapply #-}  @@ -124,6 +132,7 @@ partial :: Num a => [Int] -> Sparse a -> a partial []     (Sparse a _)  = a partial (n:ns) (Sparse _ da) = partial ns $ findWithDefault (lift 0) n da+partial _      Zero          = 0 {-# INLINE partial #-}  spartial :: Num a => [Int] -> Sparse a -> Maybe a@@ -131,42 +140,66 @@ spartial (n:ns) (Sparse _ da) = do     a' <- lookup n da     spartial ns a'+spartial _  Zero         = Nothing {-# INLINE spartial #-} -- instance Primal Sparse where     primal (Sparse a _) = a+    primal Zero = 0  instance Lifted Sparse => Mode Sparse where-    lift a = Sparse a (IntMap.empty)+    lift a = Sparse a IntMap.empty+    zero = Zero+    _ <**> Zero = lift 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 *! z /! xi, z *! log1 xi)) x y+    Zero <+> a = a+    a <+> Zero = a     Sparse a as <+> Sparse b bs = Sparse (a + b) $ unionWith (<+>) as bs+    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  instance Lifted Sparse => Jacobian Sparse where     type D Sparse = Sparse+    unary f _ Zero = lift (f 0)     unary f dadb (Sparse pb bs) = Sparse (f pb) $ mapWithKey (times dadb) bs++    lift1 f _ Zero = lift (f 0)     lift1 f df b@(Sparse pb bs) = Sparse (f pb) $ mapWithKey (times (df b)) bs++    lift1_ f _  Zero = lift (f 0)     lift1_ f df b@(Sparse pb bs) = a where         a = Sparse (f pb) $ mapWithKey (times (df a b)) bs -    binary f dadb dadc (Sparse pb db) (Sparse pc dc) = Sparse (f pb pc) $ -        unionWith (<+>) +    binary f _    _    Zero           Zero           = lift (f 0 0)+    binary f _    dadc Zero           (Sparse pc dc) = Sparse (f 0  pc) $ mapWithKey (times dadc) dc+    binary f dadb _    (Sparse pb db) Zero           = Sparse (f pb 0 ) $ mapWithKey (times dadb) db+    binary f dadb dadc (Sparse pb db) (Sparse pc dc) = Sparse (f pb pc) $+        unionWith (<+>)             (mapWithKey (times dadb) db)             (mapWithKey (times dadc) dc) +    lift2 f _  Zero             Zero = lift (f 0 0)+    lift2 f df Zero c@(Sparse pc dc) = Sparse (f 0 pc) $ mapWithKey (times dadc) dc where dadc = snd (df zero c)+    lift2 f df b@(Sparse pb db) Zero = Sparse (f pb 0) $ mapWithKey (times 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 (<+>) +        da = unionWith (<+>)             (mapWithKey (times dadb) db)             (mapWithKey (times dadc) dc)-        ++    lift2_ f _  Zero             Zero = lift (f 0 0)+    lift2_ f df b@(Sparse pb db) Zero = a where a = Sparse (f pb 0) (mapWithKey (times (fst (df a b zero))) db)+    lift2_ f df Zero c@(Sparse pc dc) = a where a = Sparse (f 0 pc) (mapWithKey (times (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 (<+>) +        da = unionWith (<+>)             (mapWithKey (times dadb) db)             (mapWithKey (times dadc) dc) 

@@ -28,7 +28,7 @@     ) where  import Prelude hiding (all)-import Control.Applicative+import Control.Applicative hiding ((<**>)) import Data.Foldable import Data.Data (Data) import Data.Typeable (Typeable)@@ -105,6 +105,9 @@ instance Lifted Tower => Mode Tower where     lift a = Tower [a]     zero = Tower []+    _ <**> Tower []  = lift 1+    x <**> Tower [y] = lift1 (**y) (\z -> (y *^ z <**> Tower [y-1])) x+    x <**> y         = lift2_ (**) (\z xi yi -> (yi *! z /! xi, z *! log1 xi)) x y      Tower [] <+> bs = bs     as <+> Tower [] = as@@ -115,7 +118,6 @@      a *^ Tower bs = Tower (map (a*) bs)     Tower as ^* b = Tower (map (*b) as)-     Tower as ^/ b = Tower (map (/b) as)  instance Lifted Tower => Jacobian Tower where

@@ -1,5 +1,5 @@ name:         ad-version:      1.3.1+version:      1.4 license:      BSD3 license-File: LICENSE copyright:    (c) Edward Kmett 2010-2011,@@ -62,7 +62,9 @@     .     Changes since 1.3     .-    * Dependency bump to be compatible with ghc 7.4.1+    * Dependency bump to be compatible with ghc 7.4.1 and mtl 2.1+    .+    * Work on diff (**2) 0     .     Changes since 1.2     .