Skip to content

Commit

Permalink
whitespace, xform u (b x y) to ((u.).b) x y to ((.).(.)) u b x y
Browse files Browse the repository at this point in the history
  • Loading branch information
barak committed Apr 10, 2009
1 parent d810e06 commit e4fd1c4
Show file tree
Hide file tree
Showing 2 changed files with 67 additions and 31 deletions.
83 changes: 57 additions & 26 deletions Numeric/FAD.hs
Original file line number Original file line Diff line number Diff line change
Expand Up @@ -282,9 +282,10 @@ liftA2_ f df x y = z
(dfdx, dfdy) = df z x y (dfdx, dfdy) = df z x y


-- | The 'liftA1disc' function lifts a scalar function with numeric -- | The 'liftA1disc' function lifts a scalar function with numeric
-- input and discrete output from into the derivative tower domain. -- input and discrete output from the primal domain into the
-- derivative tower domain.
liftA1disc :: Num a => (a -> c) -> Tower tag a -> c liftA1disc :: Num a => (a -> c) -> Tower tag a -> c
liftA1disc f x = f (primal x) liftA1disc = (. primal)


-- | The 'liftA2disc' function lifts a binary function with numeric -- | The 'liftA2disc' function lifts a binary function with numeric
-- inputs and discrete output from into the derivative tower domain. -- inputs and discrete output from into the derivative tower domain.
Expand Down Expand Up @@ -588,74 +589,99 @@ diffsUU f = fromTower . apply f
-- | The 'diffsUF' function calculates an infinite list of derivatives -- | The 'diffsUF' function calculates an infinite list of derivatives
-- of a scalar-to-nonscalar function. The 0-th element of the list is -- of a scalar-to-nonscalar function. The 0-th element of the list is
-- the primal value, the 1-st element is the first derivative, etc. -- the primal value, the 1-st element is the first derivative, etc.
diffsUF :: (Num a, Num b, Functor f, Foldable f) => (forall tag. Tower tag a -> f (Tower tag b)) -> a -> [f b] diffsUF :: (Num a, Num b, Functor f, Foldable f) =>
(forall tag. Tower tag a -> f (Tower tag b))
-> a -> [f b]
diffsUF f = transposePadF . fmap fromTower . apply f diffsUF f = transposePadF . fmap fromTower . apply f


-- | The 'diffsMU' function calculates an infinite list of derivatives -- | The 'diffsMU' function calculates an infinite list of derivatives
-- of a nonscalar-to-scalar function. The 0-th element of the list is -- of a nonscalar-to-scalar function. The 0-th element of the list is
-- the primal value, the 1-st element is the first derivative, etc. -- the primal value, the 1-st element is the first derivative, etc.
-- The input is a (possibly truncated) list of the primal, first -- The input is a (possibly truncated) list of the primal, first
-- derivative, etc, of the input. -- derivative, etc, of the input.
diffsMU :: (Num a, Num b) => (forall tag. [Tower tag a] -> Tower tag b) -> [[a]] -> [b] diffsMU :: (Num a, Num b) =>
(forall tag. [Tower tag a] -> Tower tag b)
-> [[a]] -> [b]
diffsMU f = fromTower . f . map toTower . transposePad diffsMU f = fromTower . f . map toTower . transposePad


-- | The 'diffsMF' function calculates an infinite list of derivatives -- | The 'diffsMF' function calculates an infinite list of derivatives
-- of a nonscalar-to-nonscalar function. The 0-th element of the list -- of a nonscalar-to-nonscalar function. The 0-th element of the list
-- is the primal value, the 1-st element is the first derivative, etc. -- is the primal value, the 1-st element is the first derivative, etc.
-- The input is a (possibly truncated) list of the primal, first -- The input is a (possibly truncated) list of the primal, first
-- derivative, etc, of the input. -- derivative, etc, of the input.
diffsMF :: (Num a, Num b, Functor f, Foldable f) => (forall tag. [Tower tag a] -> f (Tower tag b)) -> [[a]] -> [f b] diffsMF :: (Num a, Num b, Functor f, Foldable f) =>
(forall tag. [Tower tag a] -> f (Tower tag b))
-> [[a]] -> [f b]
diffsMF f = transposePadF . fmap fromTower . f . map toTower . transposePad diffsMF f = transposePadF . fmap fromTower . f . map toTower . transposePad


-- Variants of diffsXX names diffs0XX, which zero-pad the output list -- Variants of diffsXX names diffs0XX, which zero-pad the output list


-- | The 'diffs0UU' function is like 'diffsUU' except the output is zero padded. -- | The 'diffs0UU' function is like 'diffsUU' except the output is zero padded.
diffs0UU :: (Num a, Num b) => (forall tag. Tower tag a -> Tower tag b) -> a -> [b] diffs0UU :: (Num a, Num b) =>
(forall tag. Tower tag a -> Tower tag b)
-> a -> [b]
diffs0UU f = zeroPad . diffsUU f diffs0UU f = zeroPad . diffsUU f


-- | The 'diffs0UF' function is like 'diffsUF' except the output is zero padded. -- | The 'diffs0UF' function is like 'diffsUF' except the output is zero padded.
diffs0UF :: (Num a, Num b, Functor f, Foldable f) => (forall tag. Tower tag a -> f (Tower tag b)) -> a -> [f b] diffs0UF :: (Num a, Num b, Functor f, Foldable f) =>
(forall tag. Tower tag a -> f (Tower tag b))
-> a -> [f b]
diffs0UF f = zeroPadF . diffsUF f diffs0UF f = zeroPadF . diffsUF f


-- | The 'diffs0MU' function is like 'diffsMU' except the output is zero padded. -- | The 'diffs0MU' function is like 'diffsMU' except the output is zero padded.
diffs0MU :: (Num a, Num b) => (forall tag. [Tower tag a] -> Tower tag b) -> [[a]] -> [b] diffs0MU :: (Num a, Num b) =>
(forall tag. [Tower tag a] -> Tower tag b)
-> [[a]] -> [b]
diffs0MU f = zeroPad . diffsMU f diffs0MU f = zeroPad . diffsMU f


-- | The 'diffs0MF' function is like 'diffsMF' except the output is zero padded. -- | The 'diffs0MF' function is like 'diffsMF' except the output is zero padded.
diffs0MF :: (Num a, Num b, Functor f, Foldable f) => (forall tag. [Tower tag a] -> f (Tower tag b)) -> [[a]] -> [f b] diffs0MF :: (Num a, Num b, Functor f, Foldable f) =>
(forall tag. [Tower tag a] -> f (Tower tag b))
-> [[a]] -> [f b]
diffs0MF f = zeroPadF . diffsMF f diffs0MF f = zeroPadF . diffsMF f


-- Common access patterns -- Common access patterns


-- | The 'diff' function is a synonym for 'diffUU'. -- | The 'diff' function is a synonym for 'diffUU'.
diff :: (Num a, Num b) => (forall tag. Tower tag a -> Tower tag b) -> a -> b diff :: (Num a, Num b) =>
(forall tag. Tower tag a -> Tower tag b)
-> a -> b
diff = diffUU diff = diffUU


-- | The 'diff2' function is a synonym for 'diff2UU'. -- | The 'diff2' function is a synonym for 'diff2UU'.
diff2 :: (Num a, Num b) => (forall tag. Tower tag a -> Tower tag b) -> a -> (b, b) diff2 :: (Num a, Num b) =>
(forall tag. Tower tag a -> Tower tag b)
-> a -> (b, b)
diff2 = diff2UU diff2 = diff2UU


-- | The 'diffs' function is a synonym for 'diffsUU'. -- | The 'diffs' function is a synonym for 'diffsUU'.
diffs :: (Num a, Num b) => (forall tag. Tower tag a -> Tower tag b) -> a -> [b] diffs :: (Num a, Num b) =>
(forall tag. Tower tag a -> Tower tag b)
-> a -> [b]
diffs = diffsUU diffs = diffsUU


-- | The 'diffs0' function is a synonym for 'diffs0UU'. -- | The 'diffs0' function is a synonym for 'diffs0UU'.
diffs0 :: (Num a, Num b) => (forall tag. Tower tag a -> Tower tag b) -> a -> [b] diffs0 :: (Num a, Num b) =>
(forall tag. Tower tag a -> Tower tag b)
-> a -> [b]
diffs0 = diffs0UU diffs0 = diffs0UU


-- | The 'grad' function calculates the gradient of a -- | The 'grad' function calculates the gradient of a
-- nonscalar-to-scalar function, using n invocations of forward AD, -- nonscalar-to-scalar function, using n invocations of forward AD,
-- where n is the input dimmensionality. NOTE: this is O(n) -- where n is the input dimmensionality. NOTE: this is O(n)
-- inefficient as compared to reverse AD. -- inefficient as compared to reverse AD.
grad :: (Num a, Num b) => (forall tag. [Tower tag a] -> Tower tag b) -> [a] -> [b] grad :: (Num a, Num b) =>
(forall tag. [Tower tag a] -> Tower tag b)
-> [a] -> [b]
-- grad f = head . jacobian ((:[]) . f) -- Robot face, robot claw! -- grad f = head . jacobian ((:[]) . f) -- Robot face, robot claw!
grad f xs = map (diffMU f xs) (identity xs) grad f xs = map (diffMU f xs) (identity xs)


-- | The 'jacobian' function calcualtes the Jacobian of a -- | The 'jacobian' function calcualtes the Jacobian of a
-- nonscalar-to-nonscalar function, using n invocations of forward AD, -- nonscalar-to-nonscalar function, using n invocations of forward AD,
-- where n is the input dimmensionality. -- where n is the input dimmensionality.
jacobian :: (Num a, Num b) => jacobian :: (Num a, Num b) =>
(forall tag. [Tower tag a] -> [Tower tag b]) -> [a] -> [[b]] (forall tag. [Tower tag a] -> [Tower tag b])
-> [a] -> [[b]]
jacobian f xs = transpose $ map (diffMF f xs) (identity xs) jacobian f xs = transpose $ map (diffMF f xs) (identity xs)


-- | The 'dualToPair' function converts a tower of derivatives to a -- | The 'dualToPair' function converts a tower of derivatives to a
Expand All @@ -671,12 +697,10 @@ fdualsToPair fxs = (fmap primal fxs, fmap tangent fxs)


-- | The 'zipWithBundle' function zip two lists of numbers into a list -- | The 'zipWithBundle' function zip two lists of numbers into a list
-- of derivative towers with the given primal values andd first -- of derivative towers with the given primal values andd first
-- derivatives. Like @zipWith Bundle@ except that the two lists -- derivatives. The two lists should have the same length.
-- should be the same length.
zipWithBundle :: Num a => [a] -> [a] -> [Tower tag a] zipWithBundle :: Num a => [a] -> [a] -> [Tower tag a]
zipWithBundle [] [] = [] zipWithBundle = zipWithDefaults ((flip bundle) . lift) e e
zipWithBundle (x:xs) (y:ys) = (bundle x (lift y)):(zipWithBundle xs ys) where e = error "zipWithBundle arguments, lengths differ"
zipWithBundle _ _ = error "zipWithBundle arguments, lengths differ"


-- | The 'primalUU' function lowers a function over dual numbers to a -- | The 'primalUU' function lowers a function over dual numbers to a
-- function in the primal domain, where the function is -- function in the primal domain, where the function is
Expand All @@ -688,13 +712,15 @@ primalUU f = primal . f . lift
-- | The 'primalUF' function lowers a function over dual numbers to a -- | The 'primalUF' function lowers a function over dual numbers to a
-- function over primals, where the function is scalar-to-nonscalar. -- function over primals, where the function is scalar-to-nonscalar.
primalUF :: (Num a, Num b, Functor fb) => primalUF :: (Num a, Num b, Functor fb) =>
(forall tag. Tower tag a -> fb (Tower tag b)) -> a -> (fb b) (forall tag. Tower tag a -> fb (Tower tag b))
-> a -> (fb b)
primalUF f = fmap primal . f . lift primalUF f = fmap primal . f . lift


-- | The 'primalFU' function lowers a function over dual numbers to a -- | The 'primalFU' function lowers a function over dual numbers to a
-- function over primals where the function is nonscalar-to-scalar. -- function over primals where the function is nonscalar-to-scalar.
primalFU :: (Num a, Num b, Functor fa) => primalFU :: (Num a, Num b, Functor fa) =>
(forall tag. fa (Tower tag a) -> Tower tag b) -> (fa a) -> b (forall tag. fa (Tower tag a) -> Tower tag b)
-> (fa a) -> b
primalFU f = primal . f . fmap lift primalFU f = primal . f . fmap lift


-- | The 'primalFF' function lowers a function over dual numbers to a -- | The 'primalFF' function lowers a function over dual numbers to a
Expand Down Expand Up @@ -724,7 +750,9 @@ show2d = ("["++) . (++"]\n") . (foldl1 $ (++) . (++"\n ")) . map show
-- list of increasingly higher-order approximations. -- list of increasingly higher-order approximations.
-- --
-- EXAMPLE: @taylor exp 0 1@ -- EXAMPLE: @taylor exp 0 1@
taylor :: Fractional a => (forall tag. Tower tag a -> Tower tag a) -> a -> a -> [a] taylor :: Fractional a =>
(forall tag. Tower tag a -> Tower tag a)
-> a -> a -> [a]


taylor f x dx = scanl1 (+) taylor f x dx = scanl1 (+)
$ zipWith3 (\x y z -> x*y*z) $ zipWith3 (\x y z -> x*y*z)
Expand Down Expand Up @@ -770,7 +798,8 @@ taylor2 f x y dx dy =
-- @take 10 $ zeroNewton ((+1).(^2)) (1 :+ 1) -- converge to (0 :+ 1)@ -- @take 10 $ zeroNewton ((+1).(^2)) (1 :+ 1) -- converge to (0 :+ 1)@
-- --
zeroNewton :: Fractional a => zeroNewton :: Fractional a =>
(forall tag. Tower tag a -> Tower tag a) -> a -> [a] (forall tag. Tower tag a -> Tower tag a)
-> a -> [a]
zeroNewton f x0 = iterate (\x -> let (y,y') = diff2UU f x in x - y/y') x0 zeroNewton f x0 = iterate (\x -> let (y,y') = diff2UU f x in x - y/y') x0


-- | The 'inverseNewton' function inverts a scalar function using -- | The 'inverseNewton' function inverts a scalar function using
Expand All @@ -789,15 +818,17 @@ inverseNewton f x0 y = zeroNewton (\x -> (f x) - (lift y)) x0
-- function using Newton's method; its output is a stream of -- function using Newton's method; its output is a stream of
-- increasingly accurate results. (Modulo the usual caveats.) -- increasingly accurate results. (Modulo the usual caveats.)
fixedPointNewton :: Fractional a => fixedPointNewton :: Fractional a =>
(forall tag. Tower tag a -> Tower tag a) -> a -> [a] (forall tag. Tower tag a -> Tower tag a)
-> a -> [a]
fixedPointNewton f x0 = zeroNewton (\x -> (f x) - x) x0 fixedPointNewton f x0 = zeroNewton (\x -> (f x) - x) x0


-- | The 'extremumNewton' function finds an extremum of a scalar -- | The 'extremumNewton' function finds an extremum of a scalar
-- function using Newton's method; produces a stream of increasingly -- function using Newton's method; produces a stream of increasingly
-- accurate results. (Modulo the usual caveats.) -- accurate results. (Modulo the usual caveats.)
extremumNewton :: Fractional a => extremumNewton :: Fractional a =>
(forall tag. forall tag1. (forall tag. forall tag1.
Tower tag1 (Tower tag a) -> Tower tag1 (Tower tag a)) Tower tag1 (Tower tag a)
-> Tower tag1 (Tower tag a))
-> a -> [a] -> a -> [a]
extremumNewton f x0 = zeroNewton (diffUU f) x0 extremumNewton f x0 = zeroNewton (diffUU f) x0


Expand Down
15 changes: 10 additions & 5 deletions Test.hs
Original file line number Original file line Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ nearbyHybrid accuracy x1 x2 = abs (x1 - x2) < accuracy * maximum (map abs [x1,
infix 4 ~= infix 4 ~=


(~~=) :: (Fractional t, Ord t) => [t] -> [t] -> Bool (~~=) :: (Fractional t, Ord t) => [t] -> [t] -> Bool
(~~=) xs ys = and $ zipWithDefaults (~=) notNumber notNumber xs ys (~~=) = ((.).(.)) and $ zipWithDefaults (~=) notNumber notNumber
where notNumber = 0/0 where notNumber = 0/0
infix 4 ~~= infix 4 ~~=


Expand Down Expand Up @@ -99,10 +99,12 @@ prop_diffs_4 =


-- General routines for testing Taylor series accuracy -- General routines for testing Taylor series accuracy


taylor_accurate :: (Ord a, Fractional a) => (forall tag. Tower tag a -> Tower tag a) -> Int -> a -> a -> Bool taylor_accurate :: (Ord a, Fractional a) =>
(forall tag. Tower tag a -> Tower tag a)
-> Int -> a -> a -> Bool


taylor_accurate f n x dx = s !! 0 ~= f0 x && taylor_accurate f n x dx = s !!~ 0 ~= f0 x &&
s !!~ n ~= f0 (x+dx) s !!~ n ~= f0 (x+dx)
where s = taylor f x dx where s = taylor f x dx
f0 = primalUU f f0 = primalUU f


Expand All @@ -113,7 +115,10 @@ taylor_accurate_p f n dLo dHi x d =
dLo <= d && d <= dHi ==> taylor_accurate f n x d dLo <= d && d <= dHi ==> taylor_accurate f n x d


taylor2_accurate :: (Ord a, Fractional a) => taylor2_accurate :: (Ord a, Fractional a) =>
(forall tag0 tag. Tower tag0 (Tower tag a) -> Tower tag0 (Tower tag a) -> Tower tag0 (Tower tag a)) (forall tag0 tag.
Tower tag0 (Tower tag a)
-> Tower tag0 (Tower tag a)
-> Tower tag0 (Tower tag a))
-> Int -> Int -> a -> a -> a -> a -> Bool -> Int -> Int -> a -> a -> a -> a -> Bool


taylor2_accurate f nx ny x y dx dy = taylor2_accurate f nx ny x y dx dy =
Expand Down

0 comments on commit e4fd1c4

Please sign in to comment.