rad 0.1.4 → 0.1.6
raw patch · 2 files changed
+60/−46 lines, 2 filesdep ~data-reify
Dependency ranges changed: data-reify
Files
- Numeric/RAD.hs +58/−44
- rad.cabal +2/−2
Numeric/RAD.hs view
@@ -66,10 +66,10 @@ newtype RAD s a = RAD (Tape a (RAD s a)) data Tape a t- = C a - | V a Int- | B a a a t t- | U a a t + = Literal a + | Var a Int+ | Binary a a a t t+ | Unary a a t instance Show a => Show (RAD s a) where showsPrec d = disc1 (showsPrec d)@@ -77,27 +77,27 @@ -- | The 'lift' function injects a primal number into the RAD data type with a 0 derivative. -- If reverse-mode AD numbers formed a monad, then 'lift' would be 'return'. lift :: a -> RAD s a -lift = RAD . C +lift = RAD . Literal {-# INLINE lift #-} primal :: RAD s a -> a-primal (RAD (C y)) = y-primal (RAD (V y _)) = y-primal (RAD (B y _ _ _ _)) = y-primal (RAD (U y _ _)) = y+primal (RAD (Literal y)) = y+primal (RAD (Var y _)) = y+primal (RAD (Binary y _ _ _ _)) = y+primal (RAD (Unary y _ _)) = y {-# INLINE primal #-} var :: a -> Int -> RAD s a -var a v = RAD (V a v)+var a v = RAD (Var a v) -- TODO: A higher-order data-reify -- mapDeRef :: (Applicative f) => (forall a . Num a => RAD s a -> f (u a)) -> a -> f (Tape a (u a)) instance MuRef (RAD s a) where type DeRef (RAD s a) = Tape a- mapDeRef f (RAD (C a)) = pure (C a)- mapDeRef f (RAD (V a v)) = pure (V a v)- mapDeRef f (RAD (B a jb jc x1 x2)) = B a jb jc <$> f x1 <*> f x2- mapDeRef f (RAD (U a j x)) = U a j <$> f x+ mapDeRef f (RAD (Literal a)) = pure (Literal a)+ mapDeRef f (RAD (Var a v)) = pure (Var a v)+ mapDeRef f (RAD (Binary a jb jc x1 x2)) = Binary a jb jc <$> f x1 <*> f x2+ mapDeRef f (RAD (Unary a j x)) = Unary a j <$> f x on :: (a -> a -> c) -> (b -> a) -> b -> b -> c on f g a b = f (g a) (g b)@@ -113,25 +113,25 @@ minBound = lift minBound unary_ :: (a -> a) -> a -> RAD s a -> RAD s a-unary_ f _ (RAD (C b)) = RAD (C (f b))-unary_ f g b = RAD (U (disc1 f b) g b)+unary_ f _ (RAD (Literal b)) = RAD (Literal (f b))+unary_ f g b = RAD (Unary (disc1 f b) g b) {-# INLINE unary_ #-} unary :: (a -> a) -> (a -> a) -> RAD s a -> RAD s a-unary f _ (RAD (C b)) = RAD (C (f b))-unary f g b = RAD (U (disc1 f b) (disc1 g b) b)+unary f _ (RAD (Literal b)) = RAD (Literal (f b))+unary f g b = RAD (Unary (disc1 f b) (disc1 g b) b) {-# INLINE unary #-} binary_ :: (a -> a -> a) -> a -> a -> RAD s a -> RAD s a -> RAD s a-binary_ f _ _ (RAD (C b)) (RAD (C c)) = RAD (C (f b c))-binary_ f gb gc b c = RAD (B (f vb vc) gb gc b c)+binary_ f _ _ (RAD (Literal b)) (RAD (Literal c)) = RAD (Literal (f b c))+binary_ f gb gc b c = RAD (Binary (f vb vc) gb gc b c) where vb = primal b; vc = primal c {-# INLINE binary_ #-} -- binary_ with partials binary :: (a -> a -> a) -> (a -> a) -> (a -> a) -> RAD s a -> RAD s a -> RAD s a-binary f _ _ (RAD (C b)) (RAD (C c)) = RAD (C (f b c))-binary f gb gc b c = RAD (B (f vb vc) (gb vc) (gc vb) b c)+binary f _ _ (RAD (Literal b)) (RAD (Literal c)) = RAD (Literal (f b c))+binary f gb gc b c = RAD (Binary (f vb vc) (gb vc) (gc vb) b c) where vb = primal b; vc = primal c {-# INLINE binary #-} @@ -148,12 +148,12 @@ {-# INLINE disc3 #-} from :: Num a => RAD s a -> a -> RAD s a-from (RAD (C a)) x = RAD (C x)-from a x = RAD (U x 1 a)+from (RAD (Literal a)) x = RAD (Literal x)+from a x = RAD (Unary x 1 a) fromBy :: Num a => RAD s a -> RAD s a -> Int -> a -> RAD s a -fromBy (RAD (C a)) _ _ x = RAD (C x)-fromBy a delta n x = RAD (B x 1 (fromIntegral n) a delta)+fromBy (RAD (Literal a)) _ _ x = RAD (Literal x)+fromBy a delta n x = RAD (Binary x 1 (fromIntegral n) a delta) instance (Num a, Enum a) => Enum (RAD s a) where succ = unary_ succ 1@@ -206,16 +206,16 @@ significand x = unary_ significand (scaleFloat (- floatDigits x) 1) x - atan2 (RAD (C x)) (RAD (C y)) = RAD (C (atan2 x y))- atan2 x y = RAD (B (atan2 vx vy) (vy*r) (-vx*r) x y)+ atan2 (RAD (Literal x)) (RAD (Literal y)) = RAD (Literal (atan2 x y))+ atan2 x y = RAD (Binary (atan2 vx vy) (vy*r) (-vx*r) x y) where vx = primal x vy = primal y r = recip (vx^2 + vy^2) instance RealFrac a => RealFrac (RAD s a) where- properFraction (RAD (C a)) = (w, RAD (C p))+ properFraction (RAD (Literal a)) = (w, RAD (Literal p)) where (w, p) = properFraction a- properFraction a = (w, RAD (U p 1 a))+ properFraction a = (w, RAD (Unary p 1 a)) where (w, p) = properFraction (primal a) truncate = disc1 truncate round = disc1 truncate@@ -232,8 +232,8 @@ exp = unary exp exp log = unary log recip sqrt = unary sqrt (recip . (2*) . sqrt)- RAD (C x) ** RAD (C y) = lift (x ** y)- x ** y = RAD (B vz (vy*vz/vx) (vz*log vx) x y)+ RAD (Literal x) ** RAD (Literal y) = lift (x ** y)+ x ** y = RAD (Binary vz (vy*vz/vx) (vz*log vx) x y) where vx = primal x vy = primal y vz = vx ** vy@@ -252,11 +252,11 @@ backprop :: (Ix t, Ord t, Num a) => (Vertex -> (Tape a t, t, [t])) -> STArray s t a -> Vertex -> ST s () backprop vmap ss v = do case node of- U _ g b -> do+ Unary _ g b -> do da <- readArray ss i db <- readArray ss b writeArray ss b (db + g*da)- B _ gb gc b c -> do+ Binary _ gb gc b c -> do da <- readArray ss i db <- readArray ss b writeArray ss b (db + gb*da)@@ -266,9 +266,8 @@ where (node, i, _) = vmap v - runTape :: Num a => (Int, Int) -> RAD s a -> Array Int a -runTape vbounds tape = accumArray (+) 0 vbounds [ (id, sensitivities ! ix) | (ix, V _ id) <- xs ]+runTape vbounds tape = accumArray (+) 0 vbounds [ (id, sensitivities ! ix) | (ix, Var _ id) <- xs ] where Reified.Graph xs start = unsafePerformIO $ reifyGraph tape (g, vmap) = graphFromEdges' (edgeSet <$> filter nonConst xs)@@ -280,12 +279,27 @@ return ss sbounds ((a,_):as) = foldl' (\(lo,hi) (b,_) -> (min lo b, max hi b)) (a,a) as edgeSet (i, t) = (t, i, successors t)- nonConst (_, C{}) = False+ nonConst (_, Literal{}) = False nonConst _ = True- successors (U _ _ b) = [b]- successors (B _ _ _ b c) = [b,c]- successors _ = []+ successors (Unary _ _ b) = [b]+ successors (Binary _ _ _ b c) = [b,c]+ successors _ = [] + -- this isn't _quite_ right, as it should allow negative zeros to multiply through+ -- but then we have to know what an isNegativeZero looks like, and that rather limits+ -- our underlying data types we can permit.+ -- this approach however, allows for the occasional cycles to be resolved in the + -- dependency graph by breaking the cycle on 0 edges.++ -- test x = y where y = y * 0 + x++ -- successors (Unary _ db b) = edge db b []+ -- successors (Binary _ db dc b c) = edge db b (edge dc c [])+ -- successors _ = [] ++ -- edge 0 x xs = xs+ -- edge _ x xs = x : xs+ d :: Num a => RAD s a -> a d r = runTape (0,0) r ! 0 @@ -318,10 +332,10 @@ bind xs = (r,(0,s)) where (r,s) = runS (mapM freshVar xs) 0- freshVar a = S (\s -> let s' = s + 1 in s' `seq` (RAD (V a s), s'))+ freshVar a = S (\s -> let s' = s + 1 in s' `seq` (RAD (Var a s), s')) unbind :: Functor f => f (RAD s b) -> Array Int a -> f a -unbind xs ys = fmap (\(RAD (V _ i)) -> ys ! i) xs+unbind xs ys = fmap (\(RAD (Var _ i)) -> ys ! i) xs -- | The 'diff2UU' function calculates the value and derivative, as a -- pair, of a scalar-to-scalar function.@@ -358,7 +372,7 @@ -- significantly greater than the input dimensionality you should use 'Numeric.FAD.jacobian' instead. jacobian :: (Traversable f, Functor g, Num a) => (forall s. f (RAD s a) -> g (RAD s a)) -> f a -> g (f a) jacobian f as = unbind s . runTape bounds <$> f s- where (s,bounds) = bind as+ where (s, bounds) = bind as -- | The 'jacobian2' function calcualtes both the result and the Jacobian of a -- nonscalar-to-nonscalar function, using m invocations of reverse AD,@@ -366,7 +380,7 @@ -- 'fmap snd' on the result will recover the result of 'jacobian' jacobian2 :: (Traversable f, Functor g, Num a) => (forall s. f (RAD s a) -> g (RAD s a)) -> f a -> g (a, f a) jacobian2 f as = row <$> f s- where (s,bounds) = bind as+ where (s, bounds) = bind as row a = (primal a, unbind s (runTape bounds a)) -- | The 'zeroNewton' function finds a zero of a scalar function using
rad.cabal view
@@ -1,5 +1,5 @@ Name: rad-Version: 0.1.4+Version: 0.1.6 License: BSD3 License-File: LICENSE Copyright: Edward Kmett 2010@@ -14,7 +14,7 @@ Category: Math Build-Type: Simple Build-Depends: base >= 4 && < 5,- data-reify >= 0.5 && < 0.6, + data-reify >= 0.5 && < 0.7, containers >= 0.2 && < 0.4, array >= 0.2 && < 0.4 Exposed-Modules: Numeric.RAD