diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright (c) 2010, Edward Kmett
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above
+      copyright notice, this list of conditions and the following
+      disclaimer in the documentation and/or other materials provided
+      with the distribution.
+
+    * Neither the name of Edward Kmett nor the names of other
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/Numeric/AD.hs b/Numeric/AD.hs
new file mode 100644
--- /dev/null
+++ b/Numeric/AD.hs
@@ -0,0 +1,85 @@
+{-# LANGUAGE Rank2Types, TypeFamilies #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.AD
+-- Copyright   :  (c) Edward Kmett 2010
+-- License     :  BSD3
+-- Maintainer  :  ekmett@gmail.com
+-- Stability   :  experimental
+-- Portability :  GHC only 
+--
+-- Mixed-Mode Automatic Differentiation.
+-- 
+-- Each combinator exported from this module chooses an appropriate AD mode.
+-----------------------------------------------------------------------------
+
+module Numeric.AD 
+    ( 
+    -- * Gradients
+      grad, grad2
+
+    -- * Jacobians
+    , jacobian, jacobian2
+
+    -- * Synonyms
+    , diff
+    , diff2
+    , diffs
+    , diffs0
+
+    -- * Derivatives (Forward)
+    , diffUU
+    , diffUF
+
+    , diff2UU
+    , diff2UF
+
+    -- * Derivatives (Reverse)
+    , diffFU
+    , diff2FU
+
+    -- * Derivatives (Tower)
+    , diffsUU
+    , diffsUF
+
+    , diffs0UU
+    , diffs0UF
+
+    -- * Taylor Series (Tower)
+    , taylor
+    , taylor0
+
+    -- * Exposed Types
+    , AD(..)
+    , Mode(..)
+    ) where
+
+import Data.Traversable (Traversable)
+import Data.Foldable (Foldable, foldr')
+import Control.Applicative
+import Numeric.AD.Classes  (Mode(..))
+import Numeric.AD.Internal (AD(..), probe, unprobe)
+import Numeric.AD.Forward  (diff, diffUU, diff2, diff2UU, diffUF, diff2UF)
+import Numeric.AD.Tower    (diffsUU, diffs0UU , diffsUF, diffs0UF , diffs, diffs0, taylor, taylor0) 
+import Numeric.AD.Reverse  (diffFU, diff2FU, grad, grad2)
+
+import qualified Numeric.AD.Forward as Forward
+import qualified Numeric.AD.Reverse as Reverse
+
+-- | Calculate the Jacobian of a non-scalar-to-non-scalar function, automatically choosing between forward and reverse mode AD based on the number of inputs and outputs
+jacobian :: (Traversable f, Traversable g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f a)
+jacobian f bs = snd <$> jacobian2 f bs
+{-# INLINE jacobian #-}
+
+-- | Calculate the answer and Jacobian of a non-scalar-to-non-scalar function, automatically choosing between forward- and reverse- mode AD based on the relative, number of inputs and outputs. If you need to support functions where the output is only a 'Functor', consider using 'jacobianT' from "Numeric.AD.Forward" or 'jacobian2' from "Numeric.AD.Reverse" directly.
+jacobian2 :: (Traversable f, Traversable g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f a)
+jacobian2 f bs | n == 0    = fmap (\x -> (unprobe x, bs)) as
+               | n > m     = Reverse.jacobian2 f bs
+               | otherwise = Forward.jacobian2 f bs
+    where
+        as = f (probe <$> bs)
+        n = size bs
+        m = size as
+        size :: Foldable f => f a -> Int
+        size = foldr' (\_ b -> 1 + b) 0 
+{-# INLINE jacobian2 #-}
diff --git a/Numeric/AD/Classes.hs b/Numeric/AD/Classes.hs
new file mode 100644
--- /dev/null
+++ b/Numeric/AD/Classes.hs
@@ -0,0 +1,321 @@
+{-# LANGUAGE Rank2Types, TypeFamilies, FlexibleInstances, MultiParamTypeClasses, FlexibleContexts, FunctionalDependencies, UndecidableInstances, GeneralizedNewtypeDeriving, TemplateHaskell #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.AD.Classes
+-- Copyright   :  (c) Edward Kmett 2010
+-- License     :  BSD3
+-- Maintainer  :  ekmett@gmail.com
+-- Stability   :  experimental
+-- Portability :  GHC only 
+--
+-----------------------------------------------------------------------------
+
+module Numeric.AD.Classes
+    ( 
+    -- * AD modes
+      Mode(..) 
+    , one
+    -- * Automatically Deriving AD
+    , Jacobian(..)
+    , Primal(..)
+    , deriveLifted
+    , deriveNumeric
+    , Lifted(..)
+    ) where
+
+import Control.Applicative
+import Data.Char
+import Language.Haskell.TH
+-- import Text.Show
+
+infixl 8 **!
+infixl 7 *!, /!, ^*, *^, ^/
+infixl 6 +!, -!, <+> 
+infix 4 ==!
+
+-- | Higher-order versions of the stock numerical methods.
+class Lifted t where
+-- class Show1 t where
+    showsPrec1 :: Show a => Int -> t a -> ShowS
+--    show1 :: Show a => t a -> String
+--    showList1 :: Show a => [t a] -> String -> String
+
+-- class Eq1 t where
+    (==!) :: (Num a, Eq a) => t a -> t a -> Bool
+    -- (/=!) :: (Num a, Eq a) => t a -> t a -> Bool
+
+-- class Eq1 => Ord1 t where
+    compare1 :: (Num a, Ord a) => t a -> t a -> Ordering
+    -- (<!) :: (Num a, Ord a) => t a -> t a -> Bool
+    -- (>=!) :: (Num a, Ord a) => t a -> t a -> Bool
+    -- (>!) :: (Num a, Ord a) => t a -> t a -> Bool
+    -- (<=!) :: (Num a, Ord a) => t a -> t a -> Bool
+    -- min1 :: (Num a, Ord a) => t a -> t a -> t a
+    -- max1 :: (Num a, Ord a) => t a -> t a -> t a
+
+-- class (Show1 t, Eq t) => Num1 t where
+    fromInteger1 :: Num a => Integer -> t a
+    (+!),(-!),(*!) :: Num a => t a -> t a -> t a
+    negate1, abs1, signum1 :: Num a => t a -> t a 
+
+-- class Num1 t => Fractional1 t where
+    (/!) :: Fractional a => t a -> t a -> t a
+    recip1 :: Fractional a => t a -> t a
+    fromRational1 :: Fractional a => Rational -> t a 
+
+-- class (Num1 t, Ord1 t) => Real1 t 
+    toRational1 :: Real a => t a -> Rational -- unsafe
+
+-- class Fractional1 t => Floating1 t 
+    pi1 :: Floating a => t a 
+    exp1, log1, sqrt1 :: Floating a => t a -> t a
+    (**!), logBase1 :: Floating a => t a -> t a -> t a 
+    sin1, cos1, tan1, asin1, acos1, atan1 :: Floating a => t a -> t a
+    sinh1, cosh1, tanh1, asinh1, acosh1, atanh1 :: Floating a => t a -> t a
+
+-- class (Real1 t, Fractional1 t) => RealFrac1 t where
+    properFraction1 :: (RealFrac a, Integral b) => t a -> (b, t a)
+
+    truncate1, round1, ceiling1, floor1 :: (RealFrac a, Integral b) => t a -> b
+
+-- class (RealFrac1 t, Floating1 t) => RealFloat1 t where
+    floatRadix1 :: RealFloat a => t a -> Integer
+    floatDigits1 :: RealFloat a => t a -> Int
+    floatRange1  :: RealFloat a => t a -> (Int, Int)
+    decodeFloat1 :: RealFloat a => t a -> (Integer, Int)
+    encodeFloat1  :: RealFloat a => Integer -> Int -> t a
+    exponent1     :: RealFloat a => t a -> Int
+    significand1  :: RealFloat a => t a -> t a
+    scaleFloat1   :: RealFloat a => Int -> t a -> t a 
+    isNaN1, isInfinite1, isDenormalized1, isNegativeZero1, isIEEE1 :: RealFloat a => t a -> Bool
+    atan21 :: RealFloat a => t a -> t a -> t a
+
+-- class Enum1 t where
+    succ1, pred1    :: (Num a, Enum a) => t a -> t a
+    toEnum1         :: (Num a, Enum a) => Int -> t a
+    fromEnum1       :: (Num a, Enum a) => t a -> Int
+    enumFrom1       :: (Num a, Enum a) => t a -> [t a]
+    enumFromThen1   :: (Num a, Enum a) => t a -> t a -> [t a]
+    enumFromTo1     :: (Num a, Enum a) => t a -> t a -> [t a]
+    enumFromThenTo1 :: (Num a, Enum a) => t a -> t a -> t a -> [t a]
+
+-- class Bounded1 t where
+    minBound1 :: (Num a, Bounded a) => t a 
+    maxBound1 :: (Num a, Bounded a) => t a
+
+class Lifted t => Mode t where
+
+    -- | Embed a constant
+    lift  :: Num a => a -> t a 
+
+    -- | Vector sum
+    (<+>) :: Num a => t a -> t a -> t a
+
+    -- | Scalar-vector multiplication
+    (*^) :: Num a => a -> t a -> t a
+
+    -- | Vector-scalar multiplication
+    (^*) :: Num a => t a -> a -> t a
+
+    -- | Scalar division
+    (^/) :: Fractional a => t a -> a -> t a 
+
+    -- | > 'zero' = 'lift' 0
+    zero :: Num a => t a
+
+    a *^ b = lift a *! b
+    a ^* b = a *! lift b
+
+    a ^/ b = a ^* recip b
+
+    zero = lift 0
+
+-- | > 'one' = 'lift' 1
+one :: (Mode t, Num a) => t a
+one = lift 1
+{-# INLINE one #-}
+
+negOne :: (Mode t, Num a) => t a
+negOne = lift (-1)
+{-# INLINE negOne #-}
+
+-- | 'Primal' is used by 'deriveMode' but is not exposed 
+-- via the 'Mode' class to prevent its abuse by end users
+-- via the AD data type. 
+--
+-- It provides direct access to the result, stripped of its derivative information,
+-- but this is unsafe in general as (lift . primal) would discard derivative
+-- information. The end user is protected from accidentally using this function
+-- by the universal quantification on the various combinators we expose.
+
+class Primal t where
+    primal :: Num a => t a -> a
+
+-- | 'Jacobian' is used by 'deriveMode' but is not exposed
+-- via 'Mode' to prevent its abuse by end users
+-- via the 'AD' data type. 
+class (Mode t, Mode (D t)) => Jacobian t where
+    type D t :: * -> *
+
+    unary  :: Num a => (a -> a) -> D t a -> t a -> t a
+    lift1  :: Num a => (a -> a) -> (D t a -> D t a) -> t a -> t a
+    lift1_ :: Num a => (a -> a) -> (D t a -> D t a -> D t a) -> t a -> t a
+
+    binary :: Num a => (a -> a -> a) -> D t a -> D t a -> t a -> t a -> t a
+    lift2  :: Num a => (a -> a -> a) -> (D t a -> D t a -> (D t a, D t a)) -> t a -> t a -> t a
+    lift2_ :: Num a => (a -> a -> a) -> (D t a -> D t a -> D t a -> (D t a, D t a)) -> t a -> t a -> t a
+
+withPrimal :: (Jacobian t, Num a) => t a -> a -> t a
+withPrimal t a = unary (const a) one t
+
+fromBy :: (Jacobian t, Num a) => t a -> t a -> Int -> a -> t a 
+fromBy a delta n x = binary (\_ _ -> x) one (fromIntegral1 n) a delta
+
+fromIntegral1 :: (Integral n, Lifted t, Num a) => n -> t a
+fromIntegral1 = fromInteger1 . fromIntegral
+{-# INLINE fromIntegral1 #-}
+
+square1 :: (Lifted t, Num a) => t a -> t a
+square1 x = x *! x 
+{-# INLINE square1 #-}
+
+on :: (a -> a -> c) -> (b -> a) -> b -> b -> c
+on f g a b = f (g a) (g b)
+
+discrete1 :: (Primal t, Num a) => (a -> c) -> t a -> c
+discrete1 f x = f (primal x)
+
+discrete2 :: (Primal t, Num a) => (a -> a -> c) -> t a -> t a -> c
+discrete2 f x y = f (primal x) (primal y)
+
+discrete3 :: (Primal t, Num a) => (a -> a -> a -> d) -> t a -> t a -> t a -> d
+discrete3 f x y z = f (primal x) (primal y) (primal z)
+
+-- | @'deriveLifted' t@ provides
+--
+-- > instance Lifted $t 
+--
+-- given supplied instances for
+--
+-- > instance Lifted $t => Primal $t where ... 
+-- > instance Lifted $t => Jacobian $t where ...
+-- 
+-- The seemingly redundant @'Lifted' $t@ constraints are caused by Template Haskell staging restrictions.
+deriveLifted :: Q Type -> Q [Dec]
+deriveLifted _t = [d| 
+    instance Lifted $_t where
+        (==!)         = (==) `on` primal
+        compare1      = compare `on` primal
+        maxBound1     = lift maxBound
+        minBound1     = lift minBound
+        showsPrec1    = showsPrec
+        fromInteger1  = lift . fromInteger
+        (+!)          = (<+>) -- binary (+) one one
+        (-!)          = binary (-) one negOne -- TODO: <-> ? as it is, this might be pretty bad for Tower
+        (*!)          = lift2 (*) (\x y -> (y, x))
+        negate1       = lift1 negate (const negOne)
+        abs1          = lift1 abs signum1
+        signum1       = lift1 signum (const zero)
+        fromRational1 = lift . fromRational
+        (/!)          = lift2 (/) $ \x y -> (recip1 y, x)
+        recip1        = lift1 recip (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 
+        sin1      = lift1 sin cos1
+        cos1      = lift1 cos $ \x -> negate1 (sin1 x)
+        tan1 x    = sin1 x /! cos1 x 
+        asin1     = lift1 asin $ \x -> recip1 (sqrt1 (one -! square1 x))
+        acos1     = lift1 acos $ \x -> negate1 (recip1 (sqrt1 (one -! square1 x)))
+        atan1     = lift1 atan $ \x -> recip1 (one +! square1 x)
+        sinh1     = lift1 sinh cosh1
+        cosh1     = lift1 cosh sinh1
+        tanh1 x   = sinh1 x /! cosh1 x
+        asinh1    = lift1 asinh $ \x -> recip1 (sqrt1 (one +! square1 x))
+        acosh1    = lift1 acosh $ \x -> recip1 (sqrt1 (square1 x -! one))
+        atanh1    = lift1 atanh $ \x -> recip1 (one -! square1 x)
+    
+        succ1                 = lift1 succ (const one)
+        pred1                 = lift1 pred (const one)
+        toEnum1               = lift . toEnum
+        fromEnum1             = discrete1 fromEnum 
+        enumFrom1 a           = withPrimal a <$> discrete1 enumFrom a
+        enumFromTo1 a b       = withPrimal a <$> discrete2 enumFromTo a b
+        enumFromThen1 a b     = zipWith (fromBy a delta) [0..] $ discrete2 enumFromThen a b where delta = b -! a
+        enumFromThenTo1 a b c = zipWith (fromBy a delta) [0..] $ discrete3 enumFromThenTo a b c where delta = b -! a
+    
+        toRational1      = discrete1 toRational
+        floatRadix1      = discrete1 floatRadix
+        floatDigits1     = discrete1 floatDigits
+        floatRange1      = discrete1 floatRange
+        decodeFloat1     = discrete1 decodeFloat
+        encodeFloat1 m e = lift (encodeFloat m e)
+        isNaN1           = discrete1 isNaN
+        isInfinite1      = discrete1 isInfinite
+        isDenormalized1  = discrete1 isDenormalized
+        isNegativeZero1  = discrete1 isNegativeZero
+        isIEEE1          = discrete1 isIEEE
+        exponent1 = exponent . primal
+        scaleFloat1 n = unary (scaleFloat n) (scaleFloat1 n one) 
+        significand1 x =  unary significand (scaleFloat1 (- floatDigits1 x) one) x
+        atan21 = lift2 atan2 $ \vx vy -> let r = recip1 (square1 vx +! square1 vy) in (vy *! r, negate1 vx *! r)
+        properFraction1 a = (w, a `withPrimal` pb) where 
+             pa = primal a 
+             (w, pb) = properFraction pa
+        truncate1 = discrete1 truncate
+        round1    = discrete1 round
+        ceiling1  = discrete1 ceiling
+        floor1    = discrete1 floor 
+    |]
+    
+varA :: Q Type
+varA = varT (mkName "a") 
+
+-- | Find all the members defined in the 'Lifted' data type
+liftedMembers :: Q [String]
+liftedMembers = do
+    ClassI (ClassD _ _ _ _ ds) <- reify ''Lifted
+    return [ nameBase n | SigD n _ <- ds]
+
+-- | @'deriveNumeric' f g@ provides the following instances:
+--
+-- > instance ('Lifted' $f, 'Num' a, 'Enum' a) => 'Enum' ($g a)
+-- > instance ('Lifted' $f, 'Num' a, 'Eq' a) => 'Eq' ($g a)
+-- > instance ('Lifted' $f, 'Num' a, 'Ord' a) => 'Ord' ($g a)
+-- > instance ('Lifted' $f, 'Num a, 'Bounded' a) => 'Bounded' ($g a)
+-- > instance ('Lifted' $f, 'Show' a) => 'Show' ($g a)
+-- > instance ('Lifted' $f, 'Num' a) => 'Num' ($g a)
+-- > instance ('Lifted' $f, 'Fractional' a) => 'Fractional' ($g a)
+-- > instance ('Lifted' $f, 'Floating' a) => 'Floating' ($g a)
+-- > instance ('Lifted' $f, 'RealFloat' a) => 'RealFloat' ($g a)
+-- > instance ('Lifted' $f, 'RealFrac' a) => 'RealFrac' ($g a)
+-- > instance ('Lifted' $f, 'Real' a) => 'Real' ($g a)
+deriveNumeric :: Q Type -> Q Type -> Q [Dec]
+deriveNumeric t' t = do
+    members <- liftedMembers
+    let keep n = nameBase n `elem` members
+    xs <- lowerInstance keep (classP ''Num [varA]:) t t' `mapM` [''Enum, ''Eq, ''Ord, ''Bounded]
+    ys <- lowerInstance keep id                     t t' `mapM` [''Show, ''Num, ''Fractional, ''Floating, ''RealFloat,''RealFrac, ''Real]
+    return (xs ++ ys)
+
+lowerInstance :: (Name -> Bool) -> ([Q Pred] -> [Q Pred]) -> Q Type -> Q Type -> Name -> Q Dec
+lowerInstance p f t t' n = do
+    ClassI (ClassD _ _ _ _ ds) <- reify n
+    instanceD (cxt (f [classP ''Lifted [t], classP n [varA]]))
+              (conT n `appT` (t' `appT` varA))
+              (concatMap lower1 ds)
+    where 
+        lower1 :: Dec -> [Q Dec]
+        lower1 (SigD n' _) | p n'' = [valD (varP n') (normalB (varE n'')) []] where n'' = primed n'
+        lower1 _          = []
+
+        primed n' = mkName $ base ++ [prime]
+            where 
+                base = nameBase n'
+                h = head base
+                prime | isSymbol h || h `elem` "/*-<>" = '!'
+                      | otherwise = '1'
diff --git a/Numeric/AD/Directed.hs b/Numeric/AD/Directed.hs
new file mode 100644
--- /dev/null
+++ b/Numeric/AD/Directed.hs
@@ -0,0 +1,103 @@
+{-# LANGUAGE Rank2Types #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.AD.Directed
+-- Copyright   :  (c) Edward Kmett 2010
+-- License     :  BSD3
+-- Maintainer  :  ekmett@gmail.com
+-- Stability   :  experimental
+-- Portability :  GHC only 
+--
+-- Allows the choice of AD 'Mode' to be specified at the term level for
+-- benchmarking or more complicated usage patterns.
+-----------------------------------------------------------------------------
+
+module Numeric.AD.Directed
+    ( 
+    -- * Derivatives
+      diffUU
+    , diff2UU
+    -- * Common access patterns
+    , diff
+    , diff2
+    -- * Jacobians
+    , jacobian
+    , jacobian2
+    -- * Gradients
+    , grad
+    , grad2
+    -- * Exposed Types
+    , Direction(..)
+    , Mode(..)
+    , AD(..)
+    ) where
+
+import Prelude hiding (reverse)
+import Numeric.AD.Classes
+import Numeric.AD.Internal
+import Data.Traversable (Traversable)
+import qualified Numeric.AD.Reverse as R
+import qualified Numeric.AD.Forward as F
+import qualified Numeric.AD.Tower as T
+import qualified Numeric.AD as M
+import Data.Ix
+
+-- TODO: use a data types a la carte approach, so we can expose more methods here
+-- rather than just the intersection of all of the functionality
+data Direction 
+    = Forward 
+    | Reverse 
+    | Tower 
+    | Mixed 
+    deriving (Show, Eq, Ord, Read, Bounded, Enum, Ix)
+
+diffUU :: Num a => Direction -> (forall s. Mode s => AD s a -> AD s a) -> a -> a
+diffUU Forward = F.diffUU
+diffUU Reverse = R.diffUU
+diffUU Tower = T.diffUU
+diffUU Mixed = F.diffUU
+{-# INLINE diffUU #-}
+
+diff2UU :: Num a => Direction -> (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a) 
+diff2UU Forward = F.diff2UU
+diff2UU Reverse = R.diff2UU
+diff2UU Tower = T.diff2UU
+diff2UU Mixed = F.diff2UU
+{-# INLINE diff2UU #-}
+
+diff :: Num a => Direction -> (forall s. Mode s => AD s a -> AD s a) -> a -> a
+diff = diffUU
+{-# INLINE diff #-}
+
+diff2 :: Num a => Direction -> (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a) 
+diff2 = diff2UU
+{-# INLINE diff2 #-}
+
+jacobian :: (Traversable f, Traversable g, Num a) => Direction -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f a)
+jacobian Forward = F.jacobian
+jacobian Reverse = R.jacobian
+jacobian Tower = error "jacobian Tower: unimplemented"
+jacobian Mixed = M.jacobian
+{-# INLINE jacobian #-}
+
+jacobian2 :: (Traversable f, Traversable g, Num a) => Direction -> (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f a)
+jacobian2 Forward = F.jacobian2
+jacobian2 Reverse = R.jacobian2
+jacobian2 Tower = error "jacobian2 Tower: unimplemented"
+jacobian2 Mixed = M.jacobian2
+{-# INLINE jacobian2 #-}
+
+grad :: (Traversable f, Num a) => Direction -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f a
+grad Forward = F.grad
+grad Reverse = R.grad
+grad Tower   = error "grad Tower: unimplemented"
+grad Mixed   = M.grad
+{-# INLINE grad #-}
+
+grad2 :: (Traversable f, Num a) => Direction -> (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f a)
+grad2 Forward = F.grad2
+grad2 Reverse = R.grad2
+grad2 Tower   = error "grad2 Tower: unimplemented"
+grad2 Mixed   = M.grad2
+{-# INLINE grad2 #-}
+
diff --git a/Numeric/AD/Forward.hs b/Numeric/AD/Forward.hs
new file mode 100644
--- /dev/null
+++ b/Numeric/AD/Forward.hs
@@ -0,0 +1,100 @@
+{-# LANGUAGE Rank2Types #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.AD.Forward
+-- Copyright   :  (c) Edward Kmett 2010
+-- License     :  BSD3
+-- Maintainer  :  ekmett@gmail.com
+-- Stability   :  experimental
+-- Portability :  GHC only 
+--
+-- Forward mode automatic differentiation
+--
+-----------------------------------------------------------------------------
+
+module Numeric.AD.Forward
+    ( 
+    -- * Gradient
+      grad
+    , grad2
+    -- * Jacobian
+    , jacobian
+    , jacobian2
+    , jacobianT
+    -- * Derivatives
+    , diffUU
+    , diff2UU
+    , diffUF
+    , diff2UF
+    -- * Synonyms
+    , diff
+    , diff2
+    -- * Exposed Types
+    , AD(..)
+    , Mode(..)
+    ) where
+
+import Data.Traversable (Traversable)
+import Control.Applicative
+import Numeric.AD.Classes
+import Numeric.AD.Internal
+import Numeric.AD.Internal.Forward
+
+-- | The 'diff2' function calculates the first derivative of scalar-to-scalar function by 'Forward' 'AD'
+diff :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> a
+diff = diffUU
+{-# INLINE diff #-}
+
+-- | The 'diff2' function calculates the result and first derivative of scalar-to-scalar function by 'Forward' 'AD'
+diff2 :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a) 
+diff2 = diff2UU
+{-# INLINE diff2 #-}
+
+-- | The 'diffUU' function calculates the first derivative of a scalar-to-scalar function by 'Forward' 'AD'
+diffUU :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> a
+diffUU f a = tangent $ apply f a
+{-# INLINE diffUU #-}
+
+-- | The 'diff2UU' function calculates the result and first derivative of scalar-to-scalar function by 'Forward' 'AD'
+diff2UU :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a) 
+diff2UU f a = unbundle $ apply f a
+{-# INLINE diff2UU #-}
+
+-- | The 'diffUF' function calculates the first derivative of scalar-to-nonscalar function by 'Forward' 'AD'
+diffUF :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f a
+diffUF f a = tangent <$> apply f a
+{-# INLINE diffUF #-}
+
+-- | The 'diff2UF' function calculates the result and first derivative of a scalar-to-non-scalar function by 'Forward' 'AD'
+diff2UF :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f (a, a) 
+diff2UF f a = unbundle <$> apply f a
+{-# INLINE diff2UF #-}
+
+-- A fast, simple transposed forward jacobian
+jacobianT :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> f (g a)
+jacobianT f = bind (fmap tangent . f)
+-- jacobianT f as = fmap tangent <$> bind f as
+{-# INLINE jacobianT #-}
+
+jacobian :: (Traversable f, Traversable g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f a)
+jacobian f as = transposeWith (const id) t p
+    where 
+        (p, t) = bind2 (fmap tangent . f) as 
+{-# INLINE jacobian #-}
+
+jacobian2 :: (Traversable f, Traversable g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f a)
+jacobian2 f as = transposeWith row t p
+    where 
+        (p, t) = bind2 f as 
+        row x as' = (primal x, tangent <$> as') 
+{-# INLINE jacobian2 #-}
+
+grad :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f a
+grad f = bind (tangent . f)
+{-# INLINE grad #-}
+
+grad2 :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f a)
+grad2 f as = (primal b, tangent <$> bs)
+    where 
+        (b, bs) = bind2 f as
+{-# INLINE grad2 #-}
diff --git a/Numeric/AD/Internal.hs b/Numeric/AD/Internal.hs
new file mode 100644
--- /dev/null
+++ b/Numeric/AD/Internal.hs
@@ -0,0 +1,128 @@
+{-# LANGUAGE GeneralizedNewtypeDeriving, TemplateHaskell, FlexibleContexts, FlexibleInstances #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.AD.Internal
+-- Copyright   :  (c) Edward Kmett 2010
+-- License     :  BSD3
+-- Maintainer  :  ekmett@gmail.com
+-- Stability   :  experimental
+-- Portability :  GHC only 
+--
+-----------------------------------------------------------------------------
+module Numeric.AD.Internal
+    ( zipWithT
+    , AD(..)
+    , Id(..)
+    , probe
+    , unprobe
+    ) where
+
+import Control.Applicative
+import Language.Haskell.TH
+import Numeric.AD.Classes
+import Data.Monoid
+import Data.Traversable (Traversable, mapAccumL)
+import Data.Foldable (Foldable, toList)
+
+zipWithT :: (Foldable f, Traversable g) => (a -> b -> c) -> f a -> g b -> g c 
+zipWithT f as = snd . mapAccumL (\(a:as') b -> (as', f a b)) (toList as)
+
+-- | 'AD' serves as a common wrapper for different 'Mode' instances, exposing a traditional 
+-- numerical tower. Universal quantification is used to limit the actions in user code to 
+-- machinery that will return the same answers under all AD modes, allowing us to use modes
+-- interchangeably as both the type level \"brand\" and dictionary, providing a common API.
+newtype AD f a = AD { runAD :: f a } deriving (Lifted, Mode, Primal)
+
+let f = varT (mkName "f") in deriveNumeric (conT ''AD `appT` f) f
+
+newtype Id a = Id a deriving
+    (Eq, Ord, Show, Enum, Bounded, Num, Real, Fractional, Floating, RealFrac, RealFloat, Monoid)
+
+probe :: a -> AD Id a
+probe a = AD (Id a)
+
+unprobe :: AD Id a -> a
+unprobe (AD (Id a)) = a
+
+instance Functor Id where
+    fmap f (Id a) = Id (f a)
+
+instance Applicative Id where
+    pure = Id
+    Id f <*> Id a = Id (f a)
+
+instance Monad Id where
+    return = Id
+    Id a >>= f = f a 
+
+instance Lifted Id where
+    (==!) = (==)
+    compare1 = compare
+    showsPrec1 = showsPrec
+    fromInteger1 = fromInteger
+    (+!) = (+)
+    (-!) = (-)
+    (*!) = (*) 
+    negate1 = negate
+    abs1 = abs
+    signum1 = signum
+    (/!) = (/)
+    recip1 = recip
+    fromRational1 = fromRational
+    toRational1 = toRational   
+    pi1 = pi
+    exp1 = exp
+    log1 = log
+    sqrt1 = sqrt
+    (**!) = (**)
+    logBase1 = logBase
+    sin1 = sin
+    cos1 = cos
+    tan1 = tan
+    asin1 = asin
+    acos1 = acos
+    atan1 = atan
+    sinh1 = sinh
+    cosh1 = cosh
+    tanh1 = tanh
+    asinh1 = asinh
+    acosh1 = acosh
+    atanh1 = atanh
+    properFraction1 = properFraction
+    truncate1 = truncate
+    round1 = round
+    ceiling1 = ceiling
+    floor1 = floor
+    floatRadix1 = floatRadix
+    floatDigits1 = floatDigits
+    floatRange1 = floatRange
+    decodeFloat1 = decodeFloat
+    encodeFloat1 = encodeFloat
+    exponent1 = exponent
+    significand1 = significand
+    scaleFloat1 = scaleFloat
+    isNaN1 = isNaN
+    isInfinite1 = isInfinite
+    isDenormalized1 = isDenormalized
+    isNegativeZero1 = isNegativeZero
+    isIEEE1 = isIEEE
+    atan21 = atan2
+    succ1 = succ
+    pred1 = pred
+    toEnum1 = toEnum
+    fromEnum1 = fromEnum
+    enumFrom1 = enumFrom
+    enumFromThen1 = enumFromThen
+    enumFromTo1 = enumFromTo
+    enumFromThenTo1 = enumFromThenTo
+    minBound1 = minBound
+    maxBound1 = maxBound
+
+instance Mode Id where
+    lift = Id
+    Id a ^* b = Id (a * b)
+    a *^ Id b = Id (a * b)
+    Id a <+> Id b = Id (a + b)
+
+instance Primal Id where
+    primal (Id a) = a
diff --git a/Numeric/AD/Internal/Forward.hs b/Numeric/AD/Internal/Forward.hs
new file mode 100644
--- /dev/null
+++ b/Numeric/AD/Internal/Forward.hs
@@ -0,0 +1,112 @@
+{-# LANGUAGE Rank2Types, TypeFamilies, DeriveDataTypeable, TemplateHaskell, UndecidableInstances, BangPatterns #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.AD.Internal.Forward
+-- Copyright   :  (c) Edward Kmett 2010
+-- License     :  BSD3
+-- Maintainer  :  ekmett@gmail.com
+-- Stability   :  experimental
+-- Portability :  GHC only 
+--
+-- Unsafe and often partial combinators intended for internal usage.
+--
+-- Handle with care.
+-----------------------------------------------------------------------------
+
+module Numeric.AD.Internal.Forward
+    ( Forward(..)
+    , tangent
+    , bundle
+    , unbundle
+    , apply
+    , bind
+    , bind2
+    , transposeWith
+    ) where
+
+import Language.Haskell.TH
+import Data.Typeable
+import Data.Traversable (Traversable, mapAccumL)
+import Data.Foldable (Foldable, toList)
+import Data.Data
+import Control.Applicative
+import Numeric.AD.Classes
+import Numeric.AD.Internal
+
+data Forward a = Forward a a deriving (Show, Data, Typeable)
+
+tangent :: AD Forward a -> a
+tangent (AD (Forward _ da)) = da
+{-# INLINE tangent #-}
+
+unbundle :: AD Forward a -> (a, a)
+unbundle (AD (Forward a da)) = (a, da)
+{-# INLINE unbundle #-}
+
+bundle :: a -> a -> AD Forward a
+bundle a da = AD (Forward a da)
+{-# INLINE bundle #-}
+
+apply :: Num a => (AD Forward a -> b) -> a -> b
+apply f a = f (bundle a 1)
+{-# INLINE apply #-}
+
+instance Primal Forward where
+    primal (Forward a _) = a
+
+instance Lifted Forward => Mode Forward where
+    lift a = Forward a 0
+    Forward a da <+> Forward b db = Forward (a + b) (da + db)
+    a *^ Forward b db = Forward (a * b) (a * db)
+    Forward a da ^* b = Forward (a * b) (da * b)
+    Forward a da ^/ b = Forward (a / b) (da / b)
+
+instance Lifted Forward => Jacobian Forward where
+    type D Forward = Id
+    unary f (Id dadb) (Forward b db) = Forward (f b) (dadb * db)
+    lift1 f df (Forward b db) = Forward (f b) (dadb * db)
+        where 
+            Id dadb = df (Id 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
+    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 df (Forward b db) (Forward c dc) = Forward a da
+        where 
+            a = f b c
+            (Id dadb, Id dadc) = df (Id a) (Id b) (Id c)
+            da = dadb * db + dc * dadc
+
+deriveLifted $ conT ''Forward
+
+bind :: (Traversable f, Num a) => (f (AD Forward a) -> b) -> f a -> f b
+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)
+
+bind2 :: (Traversable f, Num a) => (f (AD Forward a) -> b) -> f a -> (b, f b)
+bind2 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)
+        b0 = f (lift <$> as)
+        dropIx ((_,b),bs) = (b,bs)
+
+-- we can't transpose arbitrary traversables, since we can't construct one out of whole cloth, and the outer
+-- traversable could be empty. So instead we use one as a 'skeleton'
+transposeWith :: (Functor f, Foldable f, Traversable g) => (b -> f a -> c) -> f (g a) -> g b -> g c
+transposeWith f as = snd . mapAccumL go xss0
+    where 
+        go xss b = (tail <$> xss, f b (head <$> xss))
+        xss0 = toList <$> as
+
diff --git a/Numeric/AD/Internal/Reverse.hs b/Numeric/AD/Internal/Reverse.hs
new file mode 100644
--- /dev/null
+++ b/Numeric/AD/Internal/Reverse.hs
@@ -0,0 +1,202 @@
+{-# LANGUAGE Rank2Types, TypeFamilies, MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, FlexibleContexts, TemplateHaskell, UndecidableInstances #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.AD.Internal.Reverse
+-- Copyright   :  (c) Edward Kmett 2010
+-- License     :  BSD3
+-- Maintainer  :  ekmett@gmail.com
+-- Stability   :  experimental
+-- Portability :  GHC only
+--
+-- Reverse-Mode Automatic Differentiation implementation details
+--
+-- For reverse mode AD we use 'System.Mem.StableName.StableName' to recover sharing information from
+-- the tape to avoid combinatorial explosion, and thus run asymptotically faster
+-- than it could without such sharing information, but the use of side-effects
+-- contained herein is benign.
+--
+-----------------------------------------------------------------------------
+
+module Numeric.AD.Internal.Reverse
+    ( Reverse(..)
+    , Tape(..)
+    , partials
+    , partialArray
+    , partialMap
+    , derivative
+    , derivative2
+    , Var(..)
+    ) where
+
+import Prelude hiding (mapM)
+import Control.Applicative (Applicative(..),(<$>))
+import Control.Monad.ST
+import Control.Monad (forM_)
+import Data.List (foldl')
+import Data.Array.ST
+import Data.Array
+import Data.IntMap (IntMap, fromListWith, findWithDefault)
+import Data.Graph (graphFromEdges', topSort, Vertex)
+import Data.Reify (reifyGraph, MuRef(..))
+import qualified Data.Reify.Graph as Reified
+import Data.Traversable (Traversable, mapM)
+import System.IO.Unsafe (unsafePerformIO)
+import Language.Haskell.TH
+
+import Numeric.AD.Classes
+import Numeric.AD.Internal
+
+-- | A @Tape@ records the information needed back propagate from the output to each input during 'Reverse' 'Mode' AD.
+data Tape a t
+    = Lift a
+    | Var a {-# UNPACK #-} !Int
+    | Binary a a a t t
+    | Unary a a t
+    deriving (Show)
+
+-- | @Reverse@ is a 'Mode' using reverse-mode automatic differentiation that provides fast 'diffFU', 'diff2FU', 'grad', 'grad2' and a fast 'jacobian' when you have a significantly smaller number of outputs than inputs.
+newtype Reverse a = Reverse (Tape a (Reverse a)) deriving (Show)
+
+instance MuRef (Reverse a) where
+    type DeRef (Reverse a) = Tape a
+
+    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
+    mapDeRef f (Reverse (Unary a dadb b)) = Unary a dadb <$> f b
+
+instance Lifted Reverse => Mode Reverse where
+    lift a = Reverse (Lift a)
+    (<+>)  = 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
+
+instance Primal Reverse where
+    primal (Reverse (Lift a)) = a
+    primal (Reverse (Var a _)) = a
+    primal (Reverse (Binary a _ _ _ _)) = a
+    primal (Reverse (Unary a _ _)) = a
+
+instance Lifted Reverse => Jacobian Reverse where
+    type D Reverse = Id
+
+    unary f _         (Reverse (Lift a)) = Reverse (Lift (f a))
+    unary f (Id dadb) b                  = Reverse (Unary (f (primal b)) dadb b)
+
+    lift1 f df b = unary f (df (Id pb)) b
+        where pb = primal b
+
+    lift1_ f df b = unary (const a) (df (Id a) (Id pb)) b
+        where pb = primal b
+              a = f pb
+
+    binary f _         _         (Reverse (Lift b)) (Reverse (Lift c)) = Reverse (Lift (f b c))
+    binary f _         (Id dadc) (Reverse (Lift b)) c                  = Reverse (Unary (f b (primal c)) dadc c)
+    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)
+
+    lift2 f df b c = binary f dadb dadc b c
+        where (dadb, dadc) = df (Id (primal b)) (Id (primal c))
+
+    lift2_ f df b c = binary (\_ _ -> a) dadb dadc b c
+        where
+            pb = primal b
+            pc = primal c
+            a = f pb pc
+            (dadb, dadc) = df (Id a) (Id pb) (Id pc)
+
+deriveLifted (conT ''Reverse)
+-- deriveNumeric  ''Reverse
+
+derivative :: Num a => AD Reverse a -> a
+derivative = sum . map snd . partials
+{-# INLINE derivative #-}
+
+derivative2 :: Num a => AD Reverse a -> (a, a)
+derivative2 r = (primal r, derivative r)
+{-# INLINE derivative2 #-}
+
+-- | back propagate sensitivities along a tape.
+backPropagate :: Num a => (Vertex -> (Tape a Int, Int, [Int])) -> STArray s Int a -> Vertex -> ST s ()
+backPropagate vmap ss v = do
+        case node of
+            Unary _ g b -> do
+                da <- readArray ss i
+                db <- readArray ss b
+                writeArray ss b (db + g*da)
+            Binary _ gb gc b c -> do
+                da <- readArray ss i
+                db <- readArray ss b
+                writeArray ss b (db + gb*da)
+                dc <- readArray ss c
+                writeArray ss c (dc + gc*da)
+            _ -> return ()
+    where
+        (node, i, _) = vmap v
+
+        -- this isn't _quite_ right, as it should allow negative zeros to multiply through
+
+-- | This returns a list of contributions to the partials.
+-- The variable ids returned in the list are likely /not/ unique!
+partials :: Num a => AD Reverse a -> [(Int, a)]
+partials (AD tape) = [ (ident, sensitivities ! ix) | (ix, Var _ ident) <- xs ]
+    where
+        Reified.Graph xs start = unsafePerformIO $ reifyGraph tape
+        (g, vmap) = graphFromEdges' (edgeSet <$> filter nonConst xs)
+        sensitivities = runSTArray $ do
+            ss <- newArray (sbounds xs) 0
+            writeArray ss start 1
+            forM_ (topSort g) $
+                backPropagate vmap ss
+            return ss
+        sbounds ((a,_):as) = foldl' (\(lo,hi) (b,_) -> (min lo b, max hi b)) (a,a) as
+        sbounds _ = undefined -- the graph can't be empty, it contains the output node!
+        edgeSet (i, t) = (t, i, successors t)
+        nonConst (_, Lift{}) = False
+        nonConst _ = True
+        successors (Unary _ _ b) = [b]
+        successors (Binary _ _ _ b c) = [b,c]
+        successors _ = []
+
+-- | Return an 'Array' of 'partials' given bounds for the variable IDs.
+partialArray :: Num a => (Int, Int) -> AD Reverse a -> Array Int a
+partialArray vbounds tape = accumArray (+) 0 vbounds (partials tape)
+{-# INLINE partialArray #-}
+
+-- | Return an 'IntMap' of sparse partials
+partialMap :: Num a => AD Reverse a -> IntMap a
+partialMap = fromListWith (+) . partials
+{-# INLINE partialMap #-}
+
+-- A simple fresh variable supply monad
+newtype S a = S { runS :: Int -> (a,Int) }
+instance Monad S where
+    return a = S (\s -> (a,s))
+    S g >>= f = S (\s -> let (a,s') = g s in runS (f a) s')
+
+-- | Used to mark variables for inspection during the reverse pass
+class Var t a | t -> a where
+    var     :: a -> Int -> t
+    fromVar :: t -> Int
+
+    bind :: Traversable f => f a -> (f t, (Int,Int))
+    unbind :: Functor f => f t -> Array Int a -> f a
+    unbindMap :: (Functor f, Num a) => f t -> IntMap a -> f a
+
+    -- TODO: tweak bounds
+    bind xs = (r,(0,hi))
+        where
+        (r,hi) = runS (mapM freshVar xs) 0
+        freshVar a = S (\s -> let s' = s + 1 in s' `seq` (var a s, s'))
+    unbind xs ys = fmap (\v -> ys ! fromVar v) xs
+    unbindMap xs ys = fmap (\v -> findWithDefault 0 (fromVar v) ys) xs
+
+instance Var (Reverse a) a where
+    var a v = Reverse (Var a v)
+    fromVar (Reverse (Var _ v)) = v
+    fromVar _ = error "fromVar: not a Var"
+
+instance Var (AD Reverse a) a where
+    var a v = AD (var a v)
+    fromVar (AD v) = fromVar v
diff --git a/Numeric/AD/Internal/Tower.hs b/Numeric/AD/Internal/Tower.hs
new file mode 100644
--- /dev/null
+++ b/Numeric/AD/Internal/Tower.hs
@@ -0,0 +1,101 @@
+{-# LANGUAGE Rank2Types, TypeFamilies, FlexibleContexts, UndecidableInstances, TemplateHaskell #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      : Numeric.AD.Tower.Internal
+-- Copyright   : (c) Edward Kmett 2010
+-- License     : BSD3
+-- Maintainer  : ekmett@gmail.com
+-- Stability   : experimental
+-- Portability : GHC only 
+--
+-----------------------------------------------------------------------------
+
+module Numeric.AD.Internal.Tower
+    ( Tower(..)
+    , zeroPad
+    , d
+    , d2
+    , tangents
+    , bundle
+    , apply
+    , getADTower
+    ) where
+
+import Numeric.AD.Classes
+import Numeric.AD.Internal
+import Language.Haskell.TH
+
+-- | @Tower@ is an AD 'Mode' that calculates a tangent tower by forward AD, and provides fast 'diffsUU', 'diffsUF'
+newtype Tower a = Tower { getTower :: [a] } deriving (Show)
+
+-- Local combinators
+
+zeroPad :: Num a => [a] -> [a]
+zeroPad xs = xs ++ repeat 0 
+{-# INLINE zeroPad #-}
+
+d :: Num a => [a] -> a
+d (_:da:_) = da
+d _ = 0
+{-# INLINE d #-}
+
+d2 :: Num a => [a] -> (a, a) 
+d2 (a:da:_) = (a, da)
+d2 (a:_)    = (a, 0)
+d2 _        = (0, 0)
+{-# INLINE d2 #-}
+
+tangents :: Tower a -> Tower a
+tangents (Tower []) = Tower []
+tangents (Tower (_:xs)) = Tower xs
+{-# INLINE tangents #-}
+
+bundle :: a -> Tower a -> Tower a
+bundle a (Tower as) = Tower (a:as)
+{-# INLINE bundle #-}
+
+apply :: Num a => (AD Tower a -> b) -> a -> b
+apply f a = f (AD (Tower [a,1]))
+{-# INLINE apply #-}
+
+getADTower :: AD Tower a -> [a]
+getADTower (AD t) = getTower t
+{-# INLINE getADTower #-}
+
+instance Primal Tower where
+    primal (Tower (x:_)) = x
+    primal _ = 0
+
+instance Lifted Tower => Mode Tower where
+    lift a = Tower [a]
+    zero = Tower []
+
+    Tower [] <+> bs = bs
+    as <+> Tower [] = as
+    Tower (a:as) <+> Tower (b:bs) = Tower (c:cs)
+        where 
+            c = a + b
+            Tower cs = Tower as <+> Tower bs
+
+    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
+    type D Tower = Tower
+    unary f dadb b = bundle (f (primal b)) (tangents b *! dadb)
+    lift1 f df b   = bundle (f (primal b)) (tangents b *! df b)
+    lift1_ f df b = a where 
+        a = bundle (f (primal b)) (tangents b *! df a b)
+
+    binary f dadb dadc b c = bundle (f (primal b) (primal c)) (tangents b *! dadb +! tangents c *! dadc) 
+    lift2 f df b c = bundle (f (primal b) (primal c)) (tangents b *! dadb +! tangents c *! dadc) where
+        (dadb, dadc) = df b c 
+    lift2_ f df b c = a where 
+        a0 = f (primal b) (primal c)
+        da = tangents b *! dadb +! tangents c *! dadc
+        a = bundle a0 da 
+        (dadb, dadc) = df a b c
+
+deriveLifted (conT ''Tower)
diff --git a/Numeric/AD/Newton.hs b/Numeric/AD/Newton.hs
new file mode 100644
--- /dev/null
+++ b/Numeric/AD/Newton.hs
@@ -0,0 +1,66 @@
+{-# LANGUAGE Rank2Types, BangPatterns #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.AD.Newton
+-- Copyright   :  (c) Edward Kmett 2010
+-- License     :  BSD3
+-- Maintainer  :  ekmett@gmail.com
+-- Stability   :  experimental
+-- Portability :  GHC only 
+--
+-----------------------------------------------------------------------------
+
+module Numeric.AD.Newton
+    ( 
+    -- * Newton's Method (Forward AD)
+      findZero
+    , inverse
+    , fixedPoint
+    , extremum
+    -- * Gradient Descent (Reverse AD)
+    , gradientDescent
+    -- * Exposed Types
+    , AD(..)
+    , Mode(..)
+    ) where
+
+import Prelude hiding (all)
+import Numeric.AD.Classes
+import Numeric.AD.Internal
+import Data.Foldable (all)
+import Data.Traversable (Traversable)
+import Numeric.AD.Forward (diff, diff2) 
+import Numeric.AD.Reverse (grad2) 
+
+findZero :: Fractional a => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]
+findZero f x0 = iterate (\x -> let (y,y') = diff2 f x in x - y/y') x0
+{-# INLINE findZero #-}
+
+inverse :: Fractional a => (forall s. Mode s => AD s a -> AD s a) -> a -> a -> [a]
+inverse f x0 y = findZero (\x -> f x - lift y) x0
+{-# INLINE inverse  #-}
+
+fixedPoint :: Fractional a => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]
+fixedPoint f = findZero (\x -> f x - x)
+{-# INLINE fixedPoint #-}
+
+extremum :: Fractional a => (forall t s. (Mode t, Mode s) => AD t (AD s a) -> AD t (AD s a)) -> a -> [a]
+extremum f x0 = findZero (diff f) x0
+{-# INLINE extremum #-}
+
+gradientDescent :: (Traversable f, Fractional a, Ord a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> [f a]
+gradientDescent f x0 = go x0 fx0 gx0 0.1 (0 :: Int)
+    where
+        (fx0, gx0) = grad2 f x0
+        go x fx gx !eta !i 
+            | eta == 0     = [] -- step size is 0
+            | fx1 > fx     = go x fx gx (eta/2) 0 -- we stepped too far
+            | all (==0) gx = [] -- gradient is 0
+            | otherwise    = x1 : if i == 10
+                                  then go x1 fx1 gx1 (eta*2) 0
+                                  else go x1 fx1 gx1 eta (i+1)
+            where
+                -- should check gx = 0 here
+                x1 = zipWithT (\xi gxi -> xi - eta * gxi) x gx
+                (fx1, gx1) = grad2 f x1
+{-# INLINE gradientDescent #-}
diff --git a/Numeric/AD/Reverse.hs b/Numeric/AD/Reverse.hs
new file mode 100644
--- /dev/null
+++ b/Numeric/AD/Reverse.hs
@@ -0,0 +1,115 @@
+{-# LANGUAGE Rank2Types, TemplateHaskell, BangPatterns #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.AD.Reverse
+-- Copyright   :  (c) Edward Kmett 2010
+-- License     :  BSD3
+-- Maintainer  :  ekmett@gmail.com
+-- Stability   :  experimental
+-- Portability :  GHC only 
+--
+-- Mixed-Mode Automatic Differentiation.
+-- 
+-- For reverse mode AD we use 'System.Mem.StableName.StableName' to recover sharing information from 
+-- the tape to avoid combinatorial explosion, and thus run asymptotically faster
+-- than it could without such sharing information, but the use of side-effects
+-- contained herein is benign.
+--
+-----------------------------------------------------------------------------
+
+module Numeric.AD.Reverse
+    ( 
+    -- * Gradient
+      grad
+    , grad2
+    -- * Jacobian
+    , jacobian
+    , jacobian2
+    -- * Derivatives
+    , diffUU
+    , diff2UU
+    , diffFU
+    , diff2FU
+    , diffUF
+    , diff2UF
+    -- * Synonyms
+    , diff
+    , diff2
+    -- * Exposed Types
+    , AD(..)
+    , Mode(..)
+    ) where
+
+import Control.Applicative ((<$>))
+import Data.Traversable (Traversable)
+
+import Numeric.AD.Classes
+import Numeric.AD.Internal
+import Numeric.AD.Internal.Reverse
+
+-- | The 'grad' function calculates the gradient of a non-scalar-to-scalar function with 'Reverse' AD in a single pass.
+grad :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f a
+grad f as = unbind vs (partialArray bds $ f vs)
+    where (vs,bds) = bind as
+{-# INLINE grad #-}
+
+-- | The 'grad2' function calculates the result and gradient of a non-scalar-to-scalar function with 'Reverse' AD in a single pass.
+grad2 :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f a)
+grad2 f as = (primal r, unbind vs $ partialArray bds r)
+    where (vs, bds) = bind as
+          r = f vs 
+{-# INLINE grad2 #-}
+
+-- | The 'jacobian' function calculates the jacobian of a non-scalar-to-non-scalar function with reverse AD lazily in @m@ passes for @m@ outputs.
+jacobian :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (f a)
+jacobian f as = unbind vs . partialArray bds <$> f vs where 
+    (vs, bds) = bind as
+{-# INLINE jacobian #-}
+
+-- | The 'jacobian2' function calculates both the result and the Jacobian of a nonscalar-to-nonscalar function, using @m@ invocations of reverse AD,
+-- where @m@ is the output dimensionality. Applying @fmap snd@ to the result will recover the result of 'jacobian'
+jacobian2 :: (Traversable f, Functor g, Num a) => (forall s. Mode s => f (AD s a) -> g (AD s a)) -> f a -> g (a, f a)
+jacobian2 f as = row <$> f vs where 
+    (vs, bds) = bind as
+    row a = (primal a, unbind vs (partialArray bds a))
+{-# INLINE jacobian2 #-}
+
+diffUU :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> a
+diffUU f a = derivative $ f (var a 0)
+{-# INLINE diffUU #-}
+
+diffUF :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f a
+diffUF f a = derivative <$> f (var a 0)
+{-# INLINE diffUF #-}
+
+-- | The 'diff2UU' function calculates the value and derivative, as a
+-- pair, of a scalar-to-scalar function.
+diff2UU :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a)
+diff2UU f a = derivative2 $ f (var a 0)
+{-# INLINE diff2UU #-}
+
+diff2UF :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f (a, a)
+diff2UF f a = derivative2 <$> f (var a 0)
+{-# INLINE diff2UF #-}
+
+diffFU :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> f a
+diffFU f as = unbind vs $ partialArray bds (f vs)
+    where (vs, bds) = bind as
+{-# INLINE diffFU #-}
+
+diff2FU :: (Traversable f, Num a) => (forall s. Mode s => f (AD s a) -> AD s a) -> f a -> (a, f a)
+diff2FU f as = (primal result, unbind vs $ partialArray bds result)
+    where (vs, bds) = bind as
+          result = f vs
+{-# INLINE diff2FU #-}
+
+-- | The 'diff' function is a synonym for 'diffUU'.
+diff :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> a
+diff = diffUU 
+{-# INLINE diff #-}
+
+-- | The 'diff2' function is a synonym for 'diff2UU'.
+diff2 :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a)
+diff2 = diff2UU
+{-# INLINE diff2 #-}
+
diff --git a/Numeric/AD/Tower.hs b/Numeric/AD/Tower.hs
new file mode 100644
--- /dev/null
+++ b/Numeric/AD/Tower.hs
@@ -0,0 +1,98 @@
+{-# LANGUAGE Rank2Types, BangPatterns #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      : Numeric.AD.Tower
+-- Copyright   : (c) Edward Kmett 2010
+-- License     : BSD3
+-- Maintainer  : ekmett@gmail.com
+-- Stability   : experimental
+-- Portability : GHC only 
+--
+-- Higher order derivatives via a \"dual number tower\".
+--
+-----------------------------------------------------------------------------
+
+module Numeric.AD.Tower
+    ( 
+    -- * Taylor Series
+      taylor, taylor0
+    , maclaurin, maclaurin0
+    -- * Derivatives
+    , diffUU
+    , diff2UU
+    , diffsUU
+    , diffs0UU
+    , diffsUF
+    , diffs0UF 
+    -- * Synonyms
+    , diffs, diffs0
+    , diff, diff2
+    -- * Exposed Types
+    , Mode(..)
+    , AD(..)
+    ) where
+
+-- TODO: argminNaiveGradient
+
+import Control.Applicative ((<$>))
+import Numeric.AD.Classes
+import Numeric.AD.Internal
+import Numeric.AD.Internal.Tower
+
+diffsUU :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]
+diffsUU f a = getADTower $ apply f a 
+{-# INLINE diffsUU #-}
+
+diffs0UU :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]
+diffs0UU f a = zeroPad (diffsUU f a)
+{-# INLINE diffs0UU #-}
+
+diffs0UF :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f [a]
+diffs0UF f a = (zeroPad . getADTower) <$> apply f a
+{-# INLINE diffs0UF #-}
+
+diffsUF :: (Functor f, Num a) => (forall s. Mode s => AD s a -> f (AD s a)) -> a -> f [a]
+diffsUF f a = getADTower <$> apply f a
+{-# INLINE diffsUF #-}
+
+diffs :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]
+diffs = diffsUU
+{-# INLINE diffs #-}
+
+diffs0 :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]
+diffs0 = diffs0UU
+{-# INLINE diffs0 #-}
+
+taylor :: Fractional a => (forall s. Mode s => AD s a -> AD s a) -> a -> a -> [a]
+taylor f x dx = go 1 1 (diffs f x)
+    where
+        go !n !acc (a:as) = a * acc : go (n + 1) (acc * dx / n) as
+        go _ _ [] = []
+
+taylor0 :: Fractional a => (forall s. Mode s => AD s a -> AD s a) -> a -> a -> [a]
+taylor0 f x dx = zeroPad (taylor f x dx)
+{-# INLINE taylor0 #-}
+
+maclaurin :: Fractional a => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]
+maclaurin f = taylor f 0 
+{-# INLINE maclaurin #-}
+
+maclaurin0 :: Fractional a => (forall s. Mode s => AD s a -> AD s a) -> a -> [a]
+maclaurin0 f = taylor0 f 0 
+{-# INLINE maclaurin0 #-}
+
+diffUU :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> a
+diffUU f a = d $ diffs f a 
+{-# INLINE diffUU #-}
+
+diff2UU :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a)
+diff2UU f a = d2 $ diffs f a 
+{-# INLINE diff2UU #-}
+
+diff :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> a
+diff = diffUU
+{-# INLINE diff #-}
+
+diff2 :: Num a => (forall s. Mode s => AD s a -> AD s a) -> a -> (a, a)
+diff2 = diff2UU
+{-# INLINE diff2 #-}
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/TODO b/TODO
new file mode 100644
--- /dev/null
+++ b/TODO
@@ -0,0 +1,13 @@
+* Implement the diffMF, etc. functionality from Numeric.FAD
+
+* Allow the type to vary within our AD data type container, in the same fashion as Numeric.FAD.
+
+    Although, while Pearlmutter and Siskind provided the functionality to permit it in derivative combinator, they provided
+    no combinators to convert, for instance, a Dual tag Float to a Dual tag Double, so that extra functionality cannot
+    currently be leveraged.
+
+* Provide the ability to use reverse mode locally on FAD inputs, i.e.
+
+    reverseCheckpoint :: (forall s. AD s a -> AD s a) -> AD t a -> AD t a 
+
+* Provide forward-on-reverse computation of Hessians
diff --git a/ad.cabal b/ad.cabal
new file mode 100644
--- /dev/null
+++ b/ad.cabal
@@ -0,0 +1,37 @@
+Name:         ad
+Version:      0.12
+License:      BSD3
+License-File: LICENSE
+Copyright:    Edward Kmett 2010
+Author:       Edward Kmett 2010
+Maintainer:   ekmett@gmail.com
+Stability:    Experimental
+Category:     Math
+Homepage:     http://comonad.com/reader/
+Synopsis:     Mixed-Mode Automatic Differentiation.
+Description:  Forward, reverse, and higher-order automatic differentiation with a common API
+
+Build-Type:   Simple
+Build-Depends:       
+    base >= 4 && < 5,
+    data-reify >= 0.5 && < 0.6, 
+    containers >= 0.2 && < 0.4,
+    template-haskell >= 2.4 && < 2.5,
+    array >= 0.2 && < 0.4
+
+Exposed-Modules:
+    Numeric.AD
+    Numeric.AD.Forward
+    Numeric.AD.Reverse
+    Numeric.AD.Tower
+    Numeric.AD.Directed
+    Numeric.AD.Newton
+    Numeric.AD.Classes
+    Numeric.AD.Internal
+    Numeric.AD.Internal.Forward
+    Numeric.AD.Internal.Reverse
+    Numeric.AD.Internal.Tower
+
+Extra-Source-Files: TODO
+
+GHC-Options: -Wall