diff --git a/Goal/Geometry.hs b/Goal/Geometry.hs
--- a/Goal/Geometry.hs
+++ b/Goal/Geometry.hs
@@ -1,14 +1,16 @@
+-- | The main module of @goal-geometry@. Import this module to use all the
+-- features provided by this library.
 module Goal.Geometry
     (
     -- * Re-Exports
-      module Goal.Geometry.Set
-    , module Goal.Geometry.Manifold
-    , module Goal.Geometry.Linear
+    module Goal.Geometry.Manifold
+    , module Goal.Geometry.Vector
     , module Goal.Geometry.Map
-    , module Goal.Geometry.Map.Multilinear
+    , module Goal.Geometry.Map.Linear
+    , module Goal.Geometry.Map.Linear.Convolutional
+    , module Goal.Geometry.Map.NeuralNetwork
     , module Goal.Geometry.Differential
-    , module Goal.Geometry.Differential.Convex
-    , module Goal.Geometry.Plot
+    , module Goal.Geometry.Differential.GradientPursuit
     ) where
 
 
@@ -17,11 +19,11 @@
 
 -- Re-exports --
 
-import Goal.Geometry.Set
 import Goal.Geometry.Manifold
-import Goal.Geometry.Linear
+import Goal.Geometry.Vector
 import Goal.Geometry.Map
-import Goal.Geometry.Map.Multilinear
+import Goal.Geometry.Map.Linear
+import Goal.Geometry.Map.Linear.Convolutional
+import Goal.Geometry.Map.NeuralNetwork
 import Goal.Geometry.Differential
-import Goal.Geometry.Differential.Convex
-import Goal.Geometry.Plot
+import Goal.Geometry.Differential.GradientPursuit
diff --git a/Goal/Geometry/Differential.hs b/Goal/Geometry/Differential.hs
--- a/Goal/Geometry/Differential.hs
+++ b/Goal/Geometry/Differential.hs
@@ -1,236 +1,212 @@
--- | This module provides tools for working with differential and Riemannian
--- geometry.
-module Goal.Geometry.Differential (
-    -- * Tangent Spaces
-    -- ** Types
-      Tangent (Tangent, removeTangent)
-    , Bundle (Bundle, removeBundle)
-    , Partials (Partials)
-    , Differentials (Differentials)
-    -- ** Functions
-    , gradientStep
-    , projectTangent
-    , tangentToBundle
-    , bundleToTangent
-    -- * Riemannian Manifolds
-    , Riemannian (metric, flat, sharp)
-    -- ** Gradient Pursuit
-    , gradientAscent
-    , vanillaGradientAscent
-    , gradientDescent
-    , vanillaGradientDescent
+{-# LANGUAGE UndecidableInstances,UndecidableSuperClasses #-}
+
+-- | Tools for modelling the differential and Riemannian geometry of a
+-- 'Manifold'.
+module Goal.Geometry.Differential
+    ( -- * Riemannian Manifolds
+      Riemannian (metric, flat, sharp)
+    , euclideanDistance
+    -- * Backpropagation
+    , Propagate (propagate)
+    , backpropagation
+    -- * Legendre Manifolds
+    , PotentialCoordinates
+    , Legendre (potential)
+    , DuallyFlat (dualPotential)
+    , canonicalDivergence
+    -- * Automatic Differentiation
+    , differential
+    , hessian
     ) where
 
 
 --- Imports ---
 
 
-import Prelude hiding (map,minimum,maximum)
-
--- Package --
+-- Goal --
 
 import Goal.Core
 
-import Goal.Geometry.Set
+import qualified Goal.Core.Vector.Storable as S
+import qualified Goal.Core.Vector.Boxed as B
+import qualified Goal.Core.Vector.Generic as G
+
 import Goal.Geometry.Manifold
-import Goal.Geometry.Linear
+import Goal.Geometry.Vector
 import Goal.Geometry.Map
-import Goal.Geometry.Map.Multilinear
+import Goal.Geometry.Map.Linear
 
 -- Qualified --
 
-import qualified Data.Vector.Storable as C
-import qualified Numeric.LinearAlgebra.HMatrix as H
+import qualified Numeric.AD as D
 
---import Data.Vector.Storable.UnsafeSerialize
 
+-- | Computes the differential of a function of the coordinates at a point using
+-- automatic differentiation.
+differential
+    :: Manifold x
+    => (forall a. RealFloat a => B.Vector (Dimension x) a -> a)
+    -> c # x
+    -> c #* x
+{-# INLINE differential #-}
+differential f = Point . G.convert . D.grad f . boxCoordinates
 
---- Differentiable Manifolds ---
+-- | Computes the Hessian of a function at a point with automatic differentiation.
+hessian
+    :: Manifold x
+    => (forall a. RealFloat a => B.Vector (Dimension x) a -> a)
+    -> c # x
+    -> c #* Tensor x x -- ^ The Hessian
+{-# INLINE hessian #-}
+hessian f p =
+    fromMatrix . S.fromRows . G.convert $ G.convert <$> D.hessian f (boxCoordinates p)
 
+-- | A class of 'Map's which can 'propagate' errors. That is, given an error
+-- derivative on the output, the input which caused the output, and a
+-- 'Map' to derive, return the derivative of the error with respect to the
+-- parameters of the 'Map', as well as the output of the 'Map'.
+class Map c f y x => Propagate c f y x where
+    propagate :: [c #* y] -- ^ The error differential
+              -> [c #* x] -- ^ A vector of inputs
+              -> c # f y x -- ^ The function to differentiate
+              -> (c #* f y x, [c # y]) -- ^ The derivative, and function output
 
--- | 'Tangent' spaces on 'Manifold's are the basis for differential geometry.
--- 'Tangent' spaces are defined at each point on a differentiable 'Manifold'.
-newtype Tangent c m = Tangent { removeTangent :: c :#: m } deriving (Eq, Read, Show)
+-- | Distance between two 'Point's based on the 'Euclidean' metric (l2 distance).
+euclideanDistance
+    :: Manifold x => c # x -> c # x -> Double
+{-# INLINE euclideanDistance #-}
+euclideanDistance xs ys = S.l2Norm (coordinates $ xs - ys)
 
--- | A 'Tangent' 'Bundle' is the original 'Manifold' combined with all its
--- 'Tangent' spaces.
-newtype Bundle c m = Bundle { removeBundle :: m } deriving (Eq, Read, Show)
+-- | An implementation of backpropagation using the 'Propagate' class. The first
+-- argument is a function which takes a generalized target output and function
+-- output and returns an error. The second argument is a list of target outputs
+-- and function inputs. The third argument is the parameteric function to be
+-- optimized, and its differential is what is returned.
+backpropagation
+    :: Propagate c f y x
+    => (a -> c # y -> c #* y)
+    -> [(a, c #* x)]
+    -> c # f y x
+    -> c #* f y x
+{-# INLINE backpropagation #-}
+backpropagation grd ysxs f =
+    let (yss,xs) = unzip ysxs
+        (df,yhts) = propagate dys xs f
+        dys = zipWith grd yss yhts
+     in df
 
--- | The 'Partials' coordinate system is defined as the partial derivatives of
--- the coordinate functions at a particular point.
-data Partials = Partials deriving (Eq, Read, Show)
 
--- | The 'Differentials' coordinate system represents the set of linear
--- functionals on the 'Tangent' space.
-data Differentials = Differentials deriving (Eq, Read, Show)
+--- Riemannian Manifolds ---
 
-gradientStep :: Manifold m => Double -> Partials :#: Tangent c m -> c :#: m
--- | 'gradientStep' follows takes a gradient in a particular tangent space and
--- transforms the point underlying the given tangent space by shifting it
--- slightly in the direction of the gradient.
-gradientStep eps f' =
-    let (Tangent p) = manifold f'
-        x' = coordinates $ eps .> f'
-     in fromCoordinates (manifold p) (coordinates p + x')
 
-projectTangent :: d :#: Tangent c m -> c :#: m
--- | Returns the underlying 'Point' from a 'Tangent' vector.
-projectTangent = removeTangent . manifold
+-- | 'Riemannian' 'Manifold's are differentiable 'Manifold's associated with a
+-- smoothly varying 'Tensor' known as the Riemannian 'metric'. 'flat' and
+-- 'sharp' correspond to applying this 'metric' to elements of the 'Primal' and
+-- 'Dual' spaces, respectively.
+class (Primal c, Manifold x) => Riemannian c x where
+    metric :: c # x -> c #* Tensor x x
+    flat :: c # x -> c # x -> c #* x
+    {-# INLINE flat #-}
+    flat p v = metric p >.> v
+    sharp :: c # x -> c #* x -> c # x
+    {-# INLINE sharp #-}
+    sharp p v = inverse (metric p) >.> v
 
-bundleToTangent :: Manifold m => c :#: Bundle d m -> c :#: Tangent d m
--- | Converts a 'Point' on a 'Tangent' 'Bundle' into a 'Tangent' vector.
-bundleToTangent p =
-    let (cs,dcs) = C.splitAt (div (dimension $ manifold p) 2) $ coordinates p
-        (Bundle m) = manifold p
-     in fromCoordinates (Tangent $ fromCoordinates m cs) dcs
 
-tangentToBundle :: Manifold m => c :#: Tangent d m -> c :#: Bundle d m
--- | Converts  a 'Tangent' vector into a 'Point' on a 'Tangent' 'Bundle'.
-tangentToBundle cm =
-    let (Tangent dm) = manifold cm
-        m = manifold dm
-     in fromCoordinates (Bundle m) $ coordinates dm C.++ coordinates cm
-
-replicatedTangents :: Manifold m => d :#: Tangent c (Replicated m) -> [d :#: Tangent c m]
--- | Converts a 'Tangent' vector on a 'Replicated' 'Manifold' into a list of
--- 'Tangent' vectors.
-replicatedTangents dp =
-    let (Tangent p) = manifold dp
-        ts = mapReplicated Tangent p
-        cs = listCoordinates dp
-     in zipWith fromList ts $ breakEvery (dimension $ head ts) cs
+--- Dually Flat Manifolds ---
 
 
--- Gradient Pursuit --
-
-gradientAscent :: (Riemannian c m, Manifold m)
-    => Double -- ^ Step size
-    -> (c :#: m -> Differentials :#: Tangent c m) -- ^ Gradient calculator
-    -> (c :#: m) -- ^ The initial point
-    -> [c :#: m] -- ^ The gradient ascent
-gradientAscent eps f' = iterate (gradientStep eps . sharp . f')
-
-vanillaGradientAscent :: Manifold m
-    => Double -- ^ Step size
-    -> (c :#: m -> Differentials :#: Tangent c m) -- ^ Gradient calculator
-    -> (c :#: m) -- ^ The initial point
-    -> [c :#: m] -- ^ The gradient ascent
-vanillaGradientAscent eps f' = iterate (gradientStep eps . breakChart . f')
-
-gradientDescent :: (Riemannian c m, Manifold m)
-    => Double -- ^ Step size
-    -> (c :#: m -> Differentials :#: Tangent c m) -- ^ Gradient calculator
-    -> (c :#: m) -- ^ The initial point
-    -> [c :#: m] -- ^ The gradient ascent
-gradientDescent eps = gradientAscent (-eps)
-
-vanillaGradientDescent :: Manifold m
-    => Double -- ^ Step size
-    -> (c :#: m -> Differentials :#: Tangent c m) -- ^ Gradient calculator
-    -> (c :#: m) -- ^ The initial point
-    -> [c :#: m] -- ^ The gradient ascent
-vanillaGradientDescent eps = vanillaGradientAscent (-eps)
-
+-- | Although convex analysis is usually developed seperately from differential
+-- geometry, it arises naturally out of the theory of dually flat 'Manifold's (<https://books.google.com/books?hl=en&lr=&id=vc2FWSo7wLUC&oi=fnd&pg=PR7&dq=methods+of+information+geometry&ots=4HsxHD_5KY&sig=gURe0tA3IEO-z-Cht_2TNsjjOG8#v=onepage&q=methods%20of%20information%20geometry&f=false Amari and Nagaoka, 2000>).
+--
+-- A 'Manifold' is 'Legendre' if it is associated with a particular convex
+-- function known as a 'potential'.
+class ( Primal (PotentialCoordinates x), Manifold x ) => Legendre x where
+    potential :: PotentialCoordinates x # x -> Double
 
---- Riemannian Manifolds ---
+-- | The (natural) coordinates of the given 'Manifold', on which the 'potential'
+-- is defined.
+type family PotentialCoordinates x :: Type
 
+-- | A 'Manifold' is 'DuallyFlat' when we can describe the 'dualPotential', which
+-- is the convex conjugate of 'potential'.
+class Legendre x => DuallyFlat x where
+    dualPotential :: PotentialCoordinates x #* x -> Double
 
--- | 'Riemannian' 'Manifold's are differentiable 'Manifold's where associated
--- with each point in the 'Manifold' is a 'Tangent' space with a smoothly
--- varying inner product. 'flat' and 'sharp' correspond to lowering and
--- raising the indices via the musical isomorphism determined by the metric
--- tensor.
---
--- A 'Riemannian' 'Manifold' should should satisfy the law
---
--- > flat $ sharp p = p
---
-class Manifold m => Riemannian c m where
-    metric :: c :#: m -> Function Partials Differentials :#: Tensor (Tangent c m) (Tangent c m)
-    flat :: Partials :#: Tangent c m -> Differentials :#: Tangent c m
-    flat p = matrixApply (metric $ projectTangent p) p
-    sharp :: Differentials :#: Tangent c m -> Partials :#: Tangent c m
-    sharp p = matrixApply (matrixInverse . metric $ projectTangent p) p
+-- | Computes the 'canonicalDivergence' between two points. Note that relative
+-- to the typical definition of the KL-Divergence/relative entropy, the
+-- arguments of this function are flipped.
+canonicalDivergence
+    :: DuallyFlat x => PotentialCoordinates x # x -> PotentialCoordinates x #* x -> Double
+{-# INLINE canonicalDivergence #-}
+canonicalDivergence pp dq = potential pp + dualPotential dq - (pp <.> dq)
 
 
 --- Instances ---
 
 
--- Replicated --
-
-instance (Manifold m, Riemannian c m) => Riemannian c (Replicated m) where
-    metric p =
-        let mtxs = mapReplicated (toHMatrix . metric) p
-         in fromHMatrix (Tensor (Tangent p) (Tangent p)) $ H.diagBlock mtxs
-    flat dp =
-        fromCoordinates (manifold dp) . C.concat $ coordinates . flat <$> replicatedTangents dp
-    sharp dp =
-        fromCoordinates (manifold dp) . C.concat $ coordinates . sharp <$> replicatedTangents dp
-
 -- Euclidean --
 
-instance Riemannian Cartesian Continuum where
-    metric p = fromList (Tensor (Tangent p) (Tangent p)) [1]
-    flat = breakChart
-    sharp = breakChart
-
-instance Riemannian Cartesian Euclidean where
-    metric p = fromHMatrix (Tensor (Tangent p) (Tangent p)) . H.ident . dimension $ manifold p
-    flat = breakChart
-    sharp = breakChart
-
--- Trivial higher order spaces --
-
-instance (Manifold m, Riemannian c m) => Riemannian Partials (Tangent c m) where
-    metric dp =
-        fromCoordinates (Tensor (Tangent dp) (Tangent dp)) . coordinates . metric $ projectTangent dp
-    sharp ddp = fromCoordinates (manifold ddp) . coordinates
-        . sharp . fromCoordinates (manifold $ projectTangent ddp) $ coordinates ddp
-    flat pdd = fromCoordinates (manifold pdd) . coordinates
-        . flat . fromCoordinates (manifold $ projectTangent pdd) $ coordinates pdd
+instance KnownNat k => Riemannian Cartesian (Euclidean k) where
+    {-# INLINE metric #-}
+    metric _ = fromMatrix S.matrixIdentity
+    {-# INLINE flat #-}
+    flat _ = breakPoint
+    {-# INLINE sharp #-}
+    sharp _ = breakPoint
 
--- Tangent Spaces --
+-- Replicated Riemannian Manifolds --
 
-instance Manifold m => Manifold (Tangent c m) where
-    dimension (Tangent p) = dimension $ manifold p
+--instance {-# OVERLAPPABLE #-} (Riemannian c x, KnownNat k) => Riemannian c (Replicated k x) where
+--    metric = error "Do not call metric on a replicated manifold"
+--    {-# INLINE flat #-}
+--    flat = S.map flat
+--    {-# INLINE sharp #-}
+--    sharp = S.map sharp
 
-instance Manifold m => Manifold (Bundle c m) where
-    dimension (Bundle m) = 2 * dimension m
+-- Backprop --
 
--- Tanget Space Coordinates --
+instance (Bilinear Tensor y x, Primal c) => Propagate c Tensor y x where
+    {-# INLINE propagate #-}
+    propagate dps qs pq = (dps >$< qs, pq >$> qs)
 
-instance Primal Partials where
-    type Dual Partials = Differentials
+--instance (Bilinear Tensor y x, Primal c) => Propagate c Tensor y x where
+--    {-# INLINE propagate #-}
+--    propagate dps qs pq =
+--        let foldfun (dp,q) (k,dpq) = (k+1,(dp >.< q) + dpq)
+--         in (uncurry (/>) . foldr foldfun (0,0) $ zip dps qs, pq >$> qs)
 
-instance Primal Differentials where
-    type Dual Differentials = Partials
+instance (Translation z y, Map c (Affine f y) z x, Propagate c f y x)
+  => Propagate c (Affine f y) z x where
+    {-# INLINE propagate #-}
+    propagate dzs xs fzx =
+        let z :: c # z
+            yx :: c # f y x
+            (z,yx) = split fzx
+            dys = anchor <$> dzs
+            (dyx,ys) = propagate dys xs yx
+         in (join (average dzs) dyx, (z >+>) <$> ys)
 
 
---- Graveyard ---
-
+-- Sums --
 
+type instance PotentialCoordinates (x,y) = PotentialCoordinates x
 
-{-
---- Functions ---
+instance (Legendre x, Legendre y, PotentialCoordinates x ~ PotentialCoordinates y)
+  => Legendre (x,y) where
+      {-# INLINE potential #-}
+      potential pmn =
+          let (pm,pn) = split pmn
+           in potential pm + potential pn
 
+type instance PotentialCoordinates (Replicated k x) = PotentialCoordinates x
 
-pushForward :: (Manifold m, Manifold n)
-    => Function c d :#: Tensor n m
-    -> c :#: m
-    -> Function Partials Partials :#: Tensor (Tangent d n) (Tangent c m)
--- | 'pushForward' takes a 'Map' between 'Manifold's and turns it into a map
--- between the 'Tangent' spaces of the 'Manifold's. Although this ought to be a
--- class, right now it's simply the trivial 'pushForward' as applied to linear
--- maps.
-pushForward pq q = fromCoordinates (Tensor (Tangent $ matrixApply pq q) (Tangent q)) $ coordinates pq
+instance (Legendre x, KnownNat k) => Legendre (Replicated k x) where
+    {-# INLINE potential #-}
+    potential ps =
+        S.sum $ mapReplicated potential ps
 
-pushForward0 :: (Manifold m, Manifold n)
-    => Function c d :#: Tensor n m
-    -> c :#: m
-    -> d :#: n
-    -> Function Partials Partials :#: Tensor (Tangent d n) (Tangent c m)
--- | 'pushForward0' takes a 'Map' between 'Manifold's and turns it into a map
--- between the 'Tangent' spaces of the 'Manifold's. In this version we can
--- specify the target space more directly.
-pushForward0 pq q p = fromCoordinates (Tensor (Tangent p) (Tangent q)) $ coordinates pq
--}
+instance (DuallyFlat x, KnownNat k) => DuallyFlat (Replicated k x) where
+    {-# INLINE dualPotential #-}
+    dualPotential ps =
+        S.sum $ mapReplicated dualPotential ps
diff --git a/Goal/Geometry/Differential/Convex.hs b/Goal/Geometry/Differential/Convex.hs
deleted file mode 100644
--- a/Goal/Geometry/Differential/Convex.hs
+++ /dev/null
@@ -1,52 +0,0 @@
--- | Tools are also provided for convex analysis, as the dual structures of
--- convex analysis are equivalent to Riemannian manifolds with certain
--- properties.
-module Goal.Geometry.Differential.Convex where
-
-
---- Imports ---
-
-
--- Goal --
-
-import Goal.Geometry.Set
-import Goal.Geometry.Manifold
-import Goal.Geometry.Linear
-import Goal.Geometry.Differential
-
---- Dually Flat Manifolds ---
-
--- | Although convex analysis is usually developed seperately from differential
--- geometry, it arrises naturally out of the theory of dually flat 'Manifold's.
---
--- A 'Manifold' is 'Legendre' for a particular coordinated system if it is
--- associated with a particular convex function on points of the manifold known
--- as a 'potential'.
-class (Primal c, Manifold m) => Legendre c m where
-    potential :: (c :#: m) -> Double
-    potentialDifferentials :: (c :#: m) -> Differentials :#: Tangent c m
-
-potentialMapping :: Legendre c m => (c :#: m) -> Dual c :#: m
-potentialMapping p = fromCoordinates (manifold p) . coordinates $ potentialDifferentials p
-
--- | Computes the 'divergence' between two points.
-divergence :: (Primal c, Legendre c m, Legendre (Dual c) m) => (c :#: m) -> (Dual c :#: m) -> Double
-divergence pp dq = potential pp + potential dq - (pp <.> dq)
-
-legendreFlat :: (Legendre c m, Riemannian c m) => c :#: m -> c :#: m -> Dual c :#: m
--- | Applies 'flat' to the second input, based on the tangent space at the first input.
-legendreFlat mp err = fromCoordinates (manifold mp) . coordinates . flat . fromCoordinates (Tangent mp) $ coordinates err
-
-
---- Instances ---
-
-
--- Generic --
-
--- Direct Sums --
-
-instance Legendre c m => Legendre c (Replicated m) where
-    potential ps = sum $ mapReplicated potential ps
-    potentialDifferentials ps =
-        let dps = mapReplicated potentialDifferentials ps
-        in fromCoordinates (Tangent ps) . coordinates $ joinReplicated dps
diff --git a/Goal/Geometry/Differential/GradientPursuit.hs b/Goal/Geometry/Differential/GradientPursuit.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Geometry/Differential/GradientPursuit.hs
@@ -0,0 +1,211 @@
+-- | Gradient pursuit-based optimization on manifolds.
+
+module Goal.Geometry.Differential.GradientPursuit
+    ( -- * Cauchy Sequences
+      cauchyLimit
+    , cauchySequence
+    -- * Gradient Pursuit
+    , vanillaGradient
+    , gradientStep
+    -- ** Algorithms
+    , GradientPursuit (Classic,Momentum,Adam)
+    , gradientPursuitStep
+    , gradientSequence
+    , vanillaGradientSequence
+    , gradientCircuit
+    , vanillaGradientCircuit
+    -- *** Defaults
+    , defaultMomentumPursuit
+    , defaultAdamPursuit
+    ) where
+
+
+--- Imports ---
+
+
+-- Goal --
+
+import Goal.Core
+
+import Goal.Geometry.Manifold
+import Goal.Geometry.Vector
+import Goal.Geometry.Differential
+
+import qualified Goal.Core.Vector.Storable as S
+
+
+--- Cauchy Sequences ---
+
+
+-- | Attempts to calculate the limit of a sequence by finding the iteration
+-- with a sufficiently small distance from its previous iteration.
+cauchyLimit
+    :: (c # x -> c # x -> Double) -- ^ Distance (divergence) from previous to next
+    -> Double -- ^ Epsilon
+    -> [c # x] -- ^ Input sequence
+    -> c # x
+{-# INLINE cauchyLimit #-}
+cauchyLimit f eps ps = last $ cauchySequence f eps ps
+
+-- | Attempts to calculate the limit of a sequence. Returns the list up to the limit.
+cauchySequence
+    :: (c # x -> c # x -> Double) -- ^ Distance (divergence) from previous to next
+    -> Double -- ^ Epsilon
+    -> [c # x] -- ^ Input list
+    -> [c # x] -- ^ Truncated list
+{-# INLINE cauchySequence #-}
+cauchySequence f eps ps =
+    let pps = takeWhile taker . zip ps $ tail ps
+     in head ps : fmap snd pps
+       where taker (p1,p2) = eps < f p1 p2
+
+
+--- Gradient Pursuit ---
+
+-- | Ignore the Riemannian metric, and convert a 'Point' from a 'Dual' space to
+-- its 'Primal' space.
+vanillaGradient :: Manifold x => c #* x -> c # x
+{-# INLINE vanillaGradient #-}
+vanillaGradient = breakPoint
+
+-- | 'gradientStep' takes a step size, a 'Point', a tangent vector at that
+-- point, and returns a 'Point' with coordinates that have moved in the
+-- direction of the tangent vector.
+gradientStep
+    :: Manifold x
+    => Double
+    -> c # x -- ^ Point
+    -> c # x -- ^ Tangent Vector
+    -> c # x -- ^ Stepped point
+{-# INLINE gradientStep #-}
+gradientStep eps (Point xs) pd =
+    Point $ xs + coordinates (eps .> pd)
+
+
+-- | An ADT reprenting three basic gradient descent algorithms.
+data GradientPursuit
+    = Classic
+    | Momentum (Int -> Double)
+    | Adam Double Double Double
+
+-- | A standard momentum schedule.
+defaultMomentumPursuit :: Double -> GradientPursuit
+{-# INLINE defaultMomentumPursuit #-}
+defaultMomentumPursuit mxmu = Momentum fmu
+    where fmu k = min mxmu $ 1 - 2**((negate 1 -) . logBase 2 . fromIntegral $ div k 250 + 1)
+
+-- | Standard Adam parameters.
+defaultAdamPursuit :: GradientPursuit
+{-# INLINE defaultAdamPursuit #-}
+defaultAdamPursuit = Adam 0.9 0.999 1e-8
+
+-- | A single step of a gradient pursuit algorithm.
+gradientPursuitStep
+    :: Manifold x
+    => Double -- ^ Learning Rate
+    -> GradientPursuit -- ^ Gradient pursuit algorithm
+    -> Int -- ^ Algorithm step
+    -> c # x -- ^ The point
+    -> c # x -- ^ The derivative
+    -> [c # x] -- ^ The velocities
+    -> (c # x, [c # x]) -- ^ The updated point and velocities
+{-# INLINE gradientPursuitStep #-}
+gradientPursuitStep eps Classic _ cp dp _ = (gradientStep eps cp dp,[])
+gradientPursuitStep eps (Momentum fmu) k cp dp (v:_) =
+    let (p,v') = momentumStep eps (fmu k) cp dp v
+     in (p,[v'])
+gradientPursuitStep eps (Adam b1 b2 rg) k cp dp (m:v:_) =
+    let (p,m',v') = adamStep eps b1 b2 rg k cp dp m v
+     in (p,[m',v'])
+gradientPursuitStep _ _ _ _ _ _ = error "Momentum list length mismatch in gradientPursuitStep"
+
+-- | Gradient ascent based on the 'Riemannian' metric.
+gradientSequence
+    :: Riemannian c x
+    => (c # x -> c #* x)  -- ^ Differential calculator
+    -> Double -- ^ Step size
+    -> GradientPursuit  -- ^ Gradient pursuit algorithm
+    -> c # x -- ^ The initial point
+    -> [c # x] -- ^ The gradient ascent
+{-# INLINE gradientSequence #-}
+gradientSequence f eps gp p0 =
+    fst <$> iterate iterator (p0,(repeat 0,0))
+        where iterator (p,(vs,k)) =
+                  let dp = sharp p $ f p
+                      (p',vs') = gradientPursuitStep eps gp k p dp vs
+                   in (p',(vs',k+1))
+
+-- | Gradient ascent which ignores the 'Riemannian' metric.
+vanillaGradientSequence
+    :: Manifold x
+    => (c # x -> c #* x)  -- ^ Differential calculator
+    -> Double -- ^ Step size
+    -> GradientPursuit  -- ^ Gradient pursuit algorithm
+    -> c # x -- ^ The initial point
+    -> [c # x] -- ^ The gradient ascent
+{-# INLINE vanillaGradientSequence #-}
+vanillaGradientSequence f eps gp p0 =
+    fst <$> iterate iterator (p0,(repeat 0,0))
+        where iterator (p,(vs,k)) =
+                  let dp = vanillaGradient $ f p
+                      (p',vs') = gradientPursuitStep eps gp k p dp vs
+                   in (p',(vs',k+1))
+
+-- | A 'Circuit' for gradient descent.
+gradientCircuit
+    :: (Monad m, Manifold x)
+    => Double -- ^ Learning Rate
+    -> GradientPursuit -- ^ Gradient pursuit algorithm
+    -> Circuit m (c # x, c # x) (c # x) -- ^ (Point, Gradient) to Updated Point
+{-# INLINE gradientCircuit #-}
+gradientCircuit eps gp = accumulateFunction (repeat 0,0) $ \(p,dp) (vs,k) -> do
+    let (p',vs') = gradientPursuitStep eps gp k p dp vs
+    return (p',(vs',k+1))
+
+-- | A 'Circuit' for gradient descent.
+vanillaGradientCircuit
+    :: (Monad m, Manifold x)
+    => Double -- ^ Learning Rate
+    -> GradientPursuit -- ^ Gradient pursuit algorithm
+    -> Circuit m (c # x, c #* x) (c # x) -- ^ (Point, Gradient) to Updated Point
+{-# INLINE vanillaGradientCircuit #-}
+vanillaGradientCircuit eps gp = second (arr vanillaGradient) >>> gradientCircuit eps gp
+
+--- Internal ---
+
+
+momentumStep
+    :: Manifold x
+    => Double -- ^ The learning rate
+    -> Double -- ^ The momentum decay
+    -> c # x -- ^ The subsequent TangentPair
+    -> c # x -- ^ The subsequent TangentPair
+    -> c # x -- ^ The current velocity
+    -> (c # x, c # x) -- ^ The (subsequent point, subsequent velocity)
+{-# INLINE momentumStep #-}
+momentumStep eps mu p fd v =
+    let v' = eps .> fd + mu .> v
+     in (gradientStep 1 p v', v')
+
+adamStep
+    :: Manifold x
+    => Double -- ^ The learning rate
+    -> Double -- ^ The first momentum rate
+    -> Double -- ^ The second momentum rate
+    -> Double -- ^ Second moment regularizer
+    -> Int -- ^ Algorithm step
+    -> c # x -- ^ The subsequent gradient
+    -> c # x -- ^ The subsequent gradient
+    -> c # x -- ^ First order velocity
+    -> c # x -- ^ Second order velocity
+    -> (c # x, c # x, c # x) -- ^ Subsequent (point, first velocity, second velocity)
+{-# INLINE adamStep #-}
+adamStep eps b1 b2 rg k0 p fd m v =
+    let k = k0+1
+        fd' = S.map (^(2 :: Int)) $ coordinates fd
+        m' = (1-b1) .> fd + b1 .> m
+        v' = (1-b2) .> Point fd' + b2 .> v
+        mhat = (1-b1^k) /> m'
+        vhat = (1-b2^k) /> v'
+        fd'' = S.zipWith (/) (coordinates mhat) . S.map ((+ rg) . sqrt) $ coordinates vhat
+     in (gradientStep eps p $ Point fd'', m',v')
diff --git a/Goal/Geometry/Linear.hs b/Goal/Geometry/Linear.hs
deleted file mode 100644
--- a/Goal/Geometry/Linear.hs
+++ /dev/null
@@ -1,108 +0,0 @@
--- | The 'Linear' module provides the tools for treating a given 'Manifold' as a
--- linear space.
-module Goal.Geometry.Linear (
-    -- * Vector Spaces
-      (<+>)
-    , (.>)
-    , (<->)
-    , (/>)
-    , meanPoint
-    -- * Dual Spaces
-    , Primal
-    , Dual
-    , (<.>)
-    ) where
-
---- Imports ---
-
-import Prelude hiding (map,minimum,maximum)
-
--- Package --
-
-import Goal.Core hiding (dot)
-import Goal.Geometry.Manifold
-
--- Unqualified --
-
-import Numeric.LinearAlgebra.HMatrix hiding (Field,(><),(<>),(<.>))
-
---import Data.Vector.Storable.UnsafeSerialize
-
-
---- Vector Spaces on Manifolds ---
-
-
-infixl 6 <+>
-(<+>) :: Manifold m => c :#: m -> c :#: m -> c :#: m
--- | Vector addition of points on a manifold.
-(<+>) p p' = fromCoordinates (manifold p) (coordinates p' + coordinates p)
-{-
-  | m == manifold p' = fromCoordinates m (coordinates p' + coordinates p)
-  | otherwise = error "Attempting to add points from distinct manifolds."
-    where m = manifold p
-          -}
-
-infixl 6 <->
-(<->) :: Manifold m => c :#: m -> c :#: m -> c :#: m
--- | Vector subtraction of points on a manifold.
-(<->) p p' = fromCoordinates (manifold p) (coordinates p - coordinates p')
-{-
-  | m == manifold p' = fromCoordinates m (coordinates p - coordinates p')
-  | otherwise = error "Attempting to subtract points from distinct manifolds."
-    where m = manifold p
-          -}
-
-
-infix 7 .>
-(.>) :: Manifold m => Double -> c :#: m -> c :#: m
--- | Scalar multiplication of points on a manifold.
-(.>) a = alterCoordinates (*a)
-
-infix 7 />
-(/>) :: Manifold m => Double -> c :#: m -> c :#: m
--- | Scalar division of points on a manifold.
-(/>) a v = recip a .> v
-
-
---- Dual Spaces ---
-
-
--- | 'Primal' charts have a 'Dual' coordinate system. The 'Dual' coordinate
--- system is the system which determines the dual basis of the dual vector
--- space via the restriction that the inner product '<.>' be the dot product.
---
--- Since finite dimensional vector spaces are isomorphic to their dual spaces
--- through the dual basis,  vector space duality is handled purely at the level
--- of coordinates in Goal -- that is, 'Primal' and 'Dual' coordinates are
--- considered different ways of describing the same fundamental objects. In
--- practice, encoding this relationship purely at the level of Charts saves a
--- great deal of computational effort.
-class (Dual (Dual c)) ~ c => Primal c where
-    type Dual c :: *
-
-infix 7 <.>
-(<.>) :: c :#: m -> Dual c :#: m -> Double
--- | '<.>' is the inner product between a dual pair of 'Point's. The defining
--- property of 'Dual' coordinate systems is that the inner product can be
--- expressed as a dot product.
-(<.>) p q = dot (coordinates p) (coordinates q)
-
--- Utility --
-
-meanPoint :: Manifold m => [c :#: m] -> c :#: m
--- | Finds the midpoint amongst a set of vectors in a convex set.
-meanPoint ps = fromCoordinates (manifold $ head ps) . mean $ coordinates <$> ps
-  {-
-  | all (== m) (manifold <$> ps) = fromCoordinates m . mean $ coordinates <$> ps
-  | otherwise = error "Attempting to add points from distinct manifolds."
-    where m = manifold $ head ps
-          -}
-
-
---- Instances ---
-
-
--- Cartesian Spaces --
-
-instance Primal Cartesian where
-    type Dual Cartesian = Cartesian
diff --git a/Goal/Geometry/Manifold.hs b/Goal/Geometry/Manifold.hs
--- a/Goal/Geometry/Manifold.hs
+++ b/Goal/Geometry/Manifold.hs
@@ -1,57 +1,47 @@
--- | This module provides the core mathematical definitions used by the rest of
--- Goal. In Goal, all mathematical structures are 'Manifold's, even when they are
--- not especially complicated ones; 'Manifold's may indicate highly articulated
--- structures, but may also indicate simpler concepts such as (vector) spaces.
---
--- 'Manifold's are sets of points which can be described locally as 'Euclidean'
--- spaces. In geometry, a point is typically a member of the actual 'Manifold'.
--- However, arbitrary types of points will often be difficult to represent
--- directly, and so points in Goal are always represented in terms of their
--- 'Coordinates' in terms of a given chart.
---
--- Charts are in turn represented by phantom types. Mathematically, charts are
--- maps between the 'Manifold' and the relevant 'Cartesian' coordinate system.
--- However, since we do not represent the points of a 'Manifold' explicility,
--- we also cannot represent Charts explicitly. As such, Atlases merely index a
--- point so as to indicate how to interpret its particular 'Coordinates'.
+{-# OPTIONS_GHC -fplugin=GHC.TypeLits.KnownNat.Solver -fplugin=GHC.TypeLits.Normalise -fconstraint-solver-iterations=10 #-}
+{-# LANGUAGE
+    UndecidableInstances,
+    StandaloneDeriving,
+    GeneralizedNewtypeDeriving
+    #-}
+-- | The core mathematical definitions used by the rest of Goal. The central
+-- object is a 'Point' on a 'Manifold'. A 'Manifold' is an object with a
+-- 'Dimension', and a 'Point' represents an element of the 'Manifold' in a
+-- particular coordinate system, represented by a chart.
 module Goal.Geometry.Manifold
     ( -- * Manifolds
-      Manifold (dimension)
-    , Transition (transition)
-    -- ** Sets
-    , Embedded (Embedded, disembed)
-    -- ** Points
-    , Coordinates
-    , (:#:) (coordinates, manifold)
-    , coordinate
-    , chart
-    , breakChart
-    , alterChart
+    Manifold (Dimension)
+    , dimension
+    -- ** Combinators
+    , Replicated
+    , R
+    -- * Points
+    , Point (Point,coordinates)
+    , type (#)
+    , breakPoint
     , listCoordinates
-    , alterCoordinates
-    , toPair
-    -- ** Charts
-    , Cartesian (Cartesian)
-    , Polar (Polar)
+    , boxCoordinates
     -- ** Constructors
-    , fromList
-    , fromCoordinates
-    , euclideanPoint
-    , realNumber
-    -- * Direct Sums
-    -- ** Replicated
-    , mapReplicated
+    , singleton
+    , fromTuple
+    , fromBoxed
+    , Product (First,Second,split,join)
+    -- ** Reshaping Points
+    , splitReplicated
     , joinReplicated
-    , concatReplicated
-    -- ** DirectSum
-    , joinPair
-    , splitPair
-    , joinPair'
-    , splitPair'
-    , joinTriple
-    , splitTriple
-    , joinTriple'
-    , splitTriple'
+    , joinBoxedReplicated
+    , mapReplicated
+    , mapReplicatedPoint
+    , splitReplicatedProduct
+    , joinReplicatedProduct
+    -- * Euclidean Manifolds
+    , Euclidean
+    -- ** Charts
+    , Cartesian
+    , Polar
+    -- ** Transition
+    , Transition (transition)
+    , transition2
     ) where
 
 
@@ -61,224 +51,279 @@
 -- Goal --
 
 import Goal.Core
-
-import Goal.Geometry.Set
+import qualified Goal.Core.Vector.Generic as G
+import qualified Goal.Core.Vector.Storable as S
+import qualified Goal.Core.Vector.Boxed as B
 
--- Qualified --
+-- Unqualified --
 
-import qualified Data.Vector.Storable as C
+import Foreign.Storable
+import Data.IndexedListLiterals
+--import Control.Parallel.Strategies
 
 
 --- Manifolds ---
 
 
--- | A geometric object with a certain 'dimension'. We assume that a 'Manifold'
--- somehow represents all the geometric, coordinate independent structure under
--- consideration. 'Manifold's should satisfy
---
--- > dimension m = length $ coordinates (Point m cs)
---
-class Eq m => Manifold m where
-    dimension :: m -> Int
+-- | A geometric object with a certain 'Dimension'.
+class KnownNat (Dimension x) => Manifold x where
+    type Dimension x :: Nat
 
--- | A point is an element of a 'Manifold' 'm' in terms of a particular
--- chart 'c'.
-data c :#: m = Point
-    { coordinates :: !Coordinates
-    , manifold :: m } deriving (Eq, Read, Show)
+dimension0 :: Manifold x => Proxy (Dimension x) -> Proxy x -> Int
+{-# INLINE dimension0 #-}
+dimension0 prxy _ = natValInt prxy
 
-infixr 1 :#:
+-- | The 'Dimension' of the given 'Manifold'.
+dimension :: Manifold x => Proxy x -> Int
+{-# INLINE dimension #-}
+dimension = dimension0 Proxy
 
-coordinate :: Int -> c :#: m -> Double
-coordinate n (Point cs _) = cs C.! n
 
-data Embedded m c = Embedded { disembed :: m } deriving (Eq, Read, Show)
+--- Points ---
 
-chart :: Manifold m => c -> c :#: m -> c :#: m
--- | 'chart' allows one to specify the Atlas of a new point. This is often
--- necessary when typeclass methods are used to generate points under a
--- variety of coordinate systems.
-chart _ = id
 
-breakChart :: Manifold m => c :#: m -> d :#: m
-breakChart p = Point (coordinates p) (manifold p)
+-- | A 'Point' on a 'Manifold'. The phantom type @m@ represents the 'Manifold', and the phantom type
+-- @c@ represents the coordinate system, or chart, in which the 'Point' is represented.
+newtype Point c x =
+    Point { coordinates :: S.Vector (Dimension x) Double }
+    deriving (Eq,Ord,Show,NFData)
 
-alterChart :: Manifold m => d -> c :#: m -> d :#: m
--- | Combines 'breakChart' and 'chart'.
-alterChart _ = breakChart
+deriving instance (KnownNat (Dimension x)) => Storable (Point c x)
+deriving instance (Manifold x, KnownNat (Dimension x)) => Floating (Point c x)
+deriving instance (Manifold x, KnownNat (Dimension x)) => Fractional (Point c x)
 
-toPair :: c :#: m -> (Double,Double)
-toPair p = (coordinate 0 p,coordinate 1 p)
+-- | An infix version of 'Point', where @x@ is assumed to be of type 'Double'.
+type (c # x) = Point c x
+infix 3 #
 
-alterCoordinates :: Manifold m => (Double -> Double) -> c :#: m -> c :#: m
--- | 'alterCoordinates' allows one to map a function over the 'coordinates' of a
--- point without changing the chart.
-alterCoordinates f (Point cs m) = Point (C.map f cs) m
+-- | Returns the coordinates of the point in list form.
+listCoordinates :: c # x -> [Double]
+{-# INLINE listCoordinates #-}
+listCoordinates = S.toList . coordinates
 
-listCoordinates :: c :#: m -> [Double]
--- | Returns the 'Coordinates' of the point in list form.
-listCoordinates (Point cs _) = C.toList cs
+-- | Returns the coordinates of the point as a boxed vector.
+boxCoordinates :: c # x -> B.Vector (Dimension x) Double
+{-# INLINE boxCoordinates #-}
+boxCoordinates =  G.convert . coordinates
 
--- | A 'transition' involves taking a point represented by the chart 'c',
--- and re-representing in terms of the chart 'd'. This will usually require
--- recomputation of the 'Coordinates'. 'Transition's should satisfy the law
---
--- > transition $ transition p = p
---
-class Transition c d m where
-    transition :: c :#: m -> d :#: m
+-- | Constructs a point with coordinates given by a boxed vector.
+fromBoxed :: B.Vector (Dimension x) Double -> c # x
+{-# INLINE fromBoxed #-}
+fromBoxed =  Point . G.convert
 
-fromList :: Manifold m => m -> [Double] -> c :#: m
--- | 'fromList' builds points without the need to work with vectors.
-fromList m cs = fromCoordinates m $ C.fromList cs
+-- | Throws away the type-level information about the chart and manifold of the
+-- given 'Point'.
+breakPoint :: Dimension x ~ Dimension y => c # x -> Point d y
+{-# INLINE breakPoint #-}
+breakPoint (Point xs) = Point xs
 
-fromCoordinates :: Manifold m => m -> Coordinates -> c :#: m
-fromCoordinates m cs -- = Point cs m
-    | dimension m == C.length cs = Point cs m
-    | otherwise = error
-        $ "Coordinate dimension (" ++ show (C.length cs) ++ ") does not match Manifold dimension (" ++ show (dimension m) ++ ")."
+-- | Constructs a 'Point' with 'Dimension' 1.
+singleton :: Dimension x ~ 1 => Double -> c # x
+{-# INLINE singleton #-}
+singleton = Point . S.singleton
 
-euclideanPoint :: [Double] -> Cartesian :#: Euclidean
--- | A convenience function for building 'Euclidean' vectors.
-euclideanPoint xs = fromList (Euclidean $ length xs) xs
+-- | Constructs a 'Point' from a tuple.
+fromTuple
+    :: ( IndexedListLiterals ds (Dimension x) Double, KnownNat (Dimension x) )
+    => ds -> c # x
+{-# INLINE fromTuple #-}
+fromTuple = Point . S.fromTuple
 
-realNumber :: Double -> Cartesian :#: Continuum
--- | A convenience function for building elements of a 'Continuum'.
-realNumber x = fromList Continuum [x]
 
---- Construction ---
+-- Manifold Combinators --
 
+-- | A 'Product' 'Manifold' is one that is produced out of the
+-- sum/product/concatenation of two source 'Manifold's.
+class ( Manifold (First z), Manifold (Second z), Manifold z
+      , Dimension z ~ (Dimension (First z) + Dimension (Second z)) )
+      => Product z where
+    -- | The 'First' 'Manifold'.
+    type First z :: Type
+    -- | The 'Second 'Manifold'.
+    type Second z :: Type
+    -- | Combine 'Point's from the 'First' and 'Second' 'Manifold' into a
+    -- 'Point' on the 'Product' 'Manifold'.
+    join :: c # First z -> c # Second z -> c # z
+    -- | Split a 'Point' on the 'Product' 'Manifold' into 'Point's from the
+    -- 'First' and 'Second' 'Manifold'.
+    split :: c # z -> (c # First z, c # Second z)
 
--- Euclidean --
+-- | A Sum type for repetitions of the same 'Manifold'.
+data Replicated (k :: Nat) m
 
--- | The 'Cartesian' coordinate system.
-data Cartesian = Cartesian
+-- | An abbreviation for 'Replicated'.
+type R k x = Replicated k x
 
--- | The 'Polar' coordinate system.
-data Polar = Polar
+-- | Splits a 'Point' on a 'Replicated' 'Manifold' into a Vector of of 'Point's.
+splitReplicated
+    :: (KnownNat k, Manifold x)
+    => c # Replicated k x
+    -> S.Vector k (c # x)
+{-# INLINE splitReplicated #-}
+splitReplicated = S.map Point . S.breakEvery . coordinates
 
--- | A function to map functions over a point on a 'Replicated' 'Manifold'.
-mapReplicated :: Manifold m => (c :#: m -> x) -> c :#: Replicated m -> [x]
-mapReplicated pf ps =
-    let (Replicated m k) = manifold ps
-        cs = coordinates ps
-        b = dimension m
-     in [ pf . fromCoordinates m $ C.slice (i * b) b cs | i <- [0.. k -1 ] ]
+-- | Joins a Vector of of 'Point's into a 'Point' on a 'Replicated' 'Manifold'.
+joinReplicated
+    :: (KnownNat k, Manifold x)
+    => S.Vector k (c # x)
+    -> c # Replicated k x
+{-# INLINE joinReplicated #-}
+joinReplicated ps = Point $ S.concatMap coordinates ps
 
-joinReplicated :: Manifold m => [c :#: m] -> c :#: Replicated m
--- | Joins a list of distributions into a 'Replicated' 'Manifold'. Be advised that this function assumes
--- that the families of the individual distributions are equal.
-joinReplicated ps =
-    Point (foldl1' (C.++) (coordinates <$> ps)) $ Replicated (manifold $ head ps) (length ps)
+-- | Joins a Vector of of 'Point's into a 'Point' on a 'Replicated' 'Manifold'.
+joinBoxedReplicated
+    :: (KnownNat k, Manifold x)
+    => B.Vector k (c # x)
+    -> c # Replicated k x
+{-# INLINE joinBoxedReplicated #-}
+joinBoxedReplicated ps = Point . S.concatMap coordinates $ G.convert ps
 
-concatReplicated :: c :#: Replicated m -> c :#: Replicated m -> c :#: Replicated m
--- | Joins two 'Replicated' 'Manifold's.
-concatReplicated (Point cs (Replicated m x)) (Point cs' (Replicated _ y)) = Point (cs C.++ cs') $ Replicated m (x + y)
+-- | A combination of 'splitReplicated' and 'fmap'.
+mapReplicated
+    :: (Storable a, KnownNat k, Manifold x)
+    => (c # x -> a) -> c # Replicated k x -> S.Vector k a
+{-# INLINE mapReplicated #-}
+mapReplicated f rp = f `S.map` splitReplicated rp
 
--- Direct Sums --
+-- | A combination of 'splitReplicated' and 'fmap', where the value of the mapped function is also a point.
+mapReplicatedPoint
+    :: (KnownNat k, Manifold x, Manifold y)
+    => (c # x -> Point d y) -> c # Replicated k x -> Point d (Replicated k y)
+{-# INLINE mapReplicatedPoint #-}
+mapReplicatedPoint f rp = Point . S.concatMap (coordinates . f) $ splitReplicated rp
 
-joinPair :: (Manifold m, Manifold n) => c :#: m -> d :#: n -> (c,d) :#: (m,n)
--- | Joins a pair of Points into a Point on the the direct sum of the underlying Charts and 'Manifold's.
-joinPair = unsafeJoinPair
+-- | Splits a 'Replicated' 'Product' 'Manifold' into a pair of 'Replicated' 'Manifold's.
+splitReplicatedProduct
+    :: (KnownNat k, Product x)
+    => c # Replicated k x
+    -> (c # Replicated k (First x), c # Replicated k (Second x))
+{-# INLINE splitReplicatedProduct #-}
+splitReplicatedProduct xys =
+    let (xs,ys) = B.unzip . B.map split . G.convert $ splitReplicated xys
+     in (joinBoxedReplicated xs, joinBoxedReplicated ys)
 
-splitPair :: (Manifold m, Manifold n) => (c,d) :#: (m,n) -> (c :#: m, d :#: n)
--- | Splits a direct sum pair.
-splitPair = unsafeSplitPair
+-- | joins a 'Replicated' 'Product' 'Manifold' out of a pair of 'Replicated' 'Manifold's.
+joinReplicatedProduct
+    :: (KnownNat k, Product x)
+    => c # Replicated k (First x)
+    -> c # Replicated k (Second x)
+    -> c # Replicated k x
+{-# INLINE joinReplicatedProduct #-}
+joinReplicatedProduct xs0 ys0 =
+    let xs = splitReplicated xs0
+        ys = splitReplicated ys0
+    in joinReplicated $ S.zipWith join xs ys
 
-joinPair' :: (Manifold m, Manifold n) => c :#: m -> c :#: n -> c :#: (m,n)
--- | Alternative version where we assume that the Charts are shared.
-joinPair' = unsafeJoinPair
+-- Charts on Euclidean Space --
 
-splitPair' :: (Manifold m, Manifold n) => c :#: (m,n) -> (c :#: m, c :#: n)
--- | Alternative version where we assume that the Charts are shared.
-splitPair' = unsafeSplitPair
+-- | @n@-dimensional Euclidean space.
+data Euclidean (n :: Nat)
 
-unsafeJoinPair :: (Manifold m, Manifold n) => c :#: m -> d :#: n -> e :#: (m,n)
-unsafeJoinPair cm dn =
-    fromCoordinates (manifold cm,manifold dn) $ coordinates cm C.++ coordinates dn
+-- | 'Cartesian' coordinates on 'Euclidean' space.
+data Cartesian
 
-unsafeSplitPair :: (Manifold m, Manifold n) => c :#: (m,n) -> (d :#: m, e :#: n)
-unsafeSplitPair cmn =
-    let (m,n) = manifold cmn
-        cs = coordinates cmn
-        (mcs,ncs) = C.splitAt (dimension m) cs
-     in (fromCoordinates m mcs, fromCoordinates n ncs)
+-- | 'Polar' coordinates on 'Euclidean' space.
+data Polar
 
-joinTriple :: (Manifold m, Manifold n, Manifold o) => c :#: m -> d :#: n -> e :#: o -> (c,d,e) :#: (m,n,o)
--- | Joins a triple of Points into a Point on the the direct sum of the underlying Charts and 'Manifold's.
-joinTriple = unsafeJoinTriple
+-- | A 'transition' involves taking a point represented by the chart c,
+-- and re-representing in terms of the chart d.
+class Transition c d x where
+    transition :: c # x -> d # x
 
-splitTriple :: (Manifold m, Manifold n, Manifold o) => (c,d,e) :#: (m,n,o) -> (c :#: m, d :#: n, e :#: o)
--- | Splits a direct sum triple.
-splitTriple = unsafeSplitTriple
+-- | Generalizes a function of two points in given coordinate systems to a
+-- function on arbitrary coordinate systems.
+transition2
+    :: (Transition cx dx x, Transition cy dy y)
+    => (dx # x -> dy # y -> a)
+    -> cx # x
+    -> cy # y
+    -> a
+{-# INLINE transition2 #-}
+transition2 f p q =
+   f (transition p) (transition q)
 
-joinTriple' :: (Manifold m, Manifold n, Manifold o) => c :#: m -> c :#: n -> c :#: o -> c :#: (m,n,o)
--- | Alternative version where we assume that the Charts are shared.
-joinTriple' = unsafeJoinTriple
 
-splitTriple' :: (Manifold m, Manifold n, Manifold o) => c :#: (m,n,o) -> (c :#: m, c :#: n, c :#: o)
--- | Alternative version where we assume that the Charts are shared.
-splitTriple' = unsafeSplitTriple
+--- Instances ---
 
-unsafeJoinTriple :: (Manifold m, Manifold n, Manifold o) => c :#: m -> d :#: n -> e :#: o -> f :#: (m,n,o)
-unsafeJoinTriple cm dn eo =
-    fromCoordinates (manifold cm, manifold dn, manifold eo) $ coordinates cm C.++ coordinates dn C.++ coordinates eo
 
-unsafeSplitTriple :: (Manifold m, Manifold n, Manifold o) => c :#: (m,n,o) -> (d :#: m, e :#: n, f :#: o)
-unsafeSplitTriple cmno =
-    let (m,n,o) = manifold cmno
-        (mcs,cs') = C.splitAt (dimension m) $ coordinates cmno
-        (ncs,ocs) = C.splitAt (dimension n) cs'
-     in (fromCoordinates m mcs, fromCoordinates n ncs, fromCoordinates o ocs)
+-- Transition --
 
 
---- Instances ---
+-- Combinators --
 
+instance Manifold x => Manifold [x] where
+    -- | The list 'Manifold' represents identical copies of the given 'Manifold'.
+    type Dimension [x] = Dimension x
 
-instance Transition c c m where
-    transition = id
+instance (Manifold x, Manifold y) => Manifold (x,y) where
+    type Dimension (x,y) = Dimension x + Dimension y
 
--- Embedded --
+instance (KnownNat k, Manifold x) => Manifold (Replicated k x) where
+    type Dimension (Replicated k x) = k * Dimension x
 
-instance Manifold m => Set (Embedded m c) where
-    type Element (Embedded m c) = c :#: m
+instance (Manifold x, Manifold y) => Product (x,y) where
+    type First (x,y) = x
+    type Second (x,y) = y
+    {-# INLINE split #-}
+    split (Point xs) =
+        let (xms,xns) = S.splitAt xs
+         in (Point xms, Point xns)
+    {-# INLINE join #-}
+    join (Point xms) (Point xns) =
+        Point $ xms S.++ xns
 
--- Euclidean --
 
-instance Manifold Euclidean where
-    dimension (Euclidean n) = n
+-- Euclidean Space --
 
-instance Manifold Continuum where
-    dimension _ = 1
+instance (KnownNat k) => Manifold (Euclidean k) where
+    type Dimension (Euclidean k) = k
 
-instance Transition Polar Cartesian Euclidean where
-    transition p =
-        let r:phis = listCoordinates p
-            phiss = reverse . tails $ reverse phis
-            m = manifold p
-            xs = [ r * cos phi * product (sin <$> phis') | (phi,phis') <- zip phis phiss ]
-         in fromList m $ xs ++ [r * product (sin <$> phis)]
+instance Transition Polar Cartesian (Euclidean 2) where
+    {-# INLINE transition #-}
+    transition rphi =
+        let [r,phi] = listCoordinates rphi
+            x = r * cos phi
+            y = r * sin phi
+         in fromTuple (x,y)
 
-instance Transition Cartesian Polar Euclidean where
-    transition p =
-        let (Euclidean n) = manifold p
-            xs = listCoordinates p
-            xs2 = listCoordinates $ alterCoordinates (^2) p
-            r = sqrt $ sum xs2
-            (phis,phin0:_) = splitAt (n-2) [ acos $ xi / sqrt (sum xs2i) | (xi,xs2i) <- zip xs (tails xs2) ]
-            xn = last xs
-            phin = if xn > 0 then phin0 else 2*pi - phin0
-         in fromList (Euclidean n) $ r : (phis ++ [phin])
+instance Transition Cartesian Polar (Euclidean 2) where
+    {-# INLINE transition #-}
+    transition xy =
+        let [x,y] = listCoordinates xy
+            r = sqrt $ (x*x) + (y*y)
+            phi = atan2 y x
+         in fromTuple (r,phi)
 
--- DirectSum --
 
-instance (Manifold m, Manifold n) => Manifold (m,n) where
-    dimension (m,n) = dimension m + dimension n
+--- Transitions ---
 
-instance (Manifold m, Manifold n, Manifold o) => Manifold (m,n,o) where
-    dimension (m,n,o) = dimension m + dimension n + dimension o
 
--- Replicated --
+instance (Manifold x, Manifold y, Transition c d x, Transition c d y)
+  => Transition c d (x,y) where
+    {-# INLINE transition #-}
+    transition cxy =
+        let (cx,cy) = split cxy
+         in join (transition cx) (transition cy)
 
-instance Manifold m => Manifold (Replicated m) where
-    dimension (Replicated m rn) = dimension m * rn
+instance (KnownNat k, Manifold x, Transition c d x) => Transition c d (Replicated k x) where
+    {-# INLINE transition #-}
+    transition = mapReplicatedPoint transition
+
+
+--- Numeric Classes ---
+
+
+instance (Manifold x, KnownNat (Dimension x)) => Num (c # x) where
+    {-# INLINE (+) #-}
+    (+) (Point xs) (Point xs') = Point $ S.add xs xs'
+    {-# INLINE (*) #-}
+    (*) (Point xs) (Point xs') = Point $ xs * xs'
+    {-# INLINE negate #-}
+    negate (Point xs) = Point $ S.scale (-1) xs
+    {-# INLINE abs #-}
+    abs (Point xs) = Point $ abs xs
+    {-# INLINE signum #-}
+    signum (Point xs) = Point $ signum xs
+    {-# INLINE fromInteger #-}
+    fromInteger x = Point . S.replicate $ fromInteger x
+
diff --git a/Goal/Geometry/Map.hs b/Goal/Geometry/Map.hs
--- a/Goal/Geometry/Map.hs
+++ b/Goal/Geometry/Map.hs
@@ -1,16 +1,8 @@
--- | The Map module provides tools for developing function space 'Manifold's.
--- A map is a 'Manifold' where the 'Point's of the Manifold represent
--- parametric functions between 'Manifold's. The defining feature of 'Map's is
--- that they have a particular 'Domain' and 'Codomain', which themselves are
--- 'Manifold's.
+-- | Definitions for working with manifolds of functions, a.k.a. function spaces.
 
 module Goal.Geometry.Map (
-    -- * Maps
-      Map (Domain, domain, Codomain, codomain)
-    , Apply ((>.>), (>$>))
-    -- * Map Charts
-    , Function (Function)
-    ) where
+     Map ((>.>),(>$>))
+     ) where
 
 
 --- Imports ---
@@ -19,58 +11,16 @@
 -- Goal --
 
 import Goal.Geometry.Manifold
-
---- Maps between Manifolds ---
+import Goal.Geometry.Vector
 
 -- Charts on Maps --
 
-data Function c d = Function c d
--- | 'Function' Charts help track Charts on the 'Domain' and 'Codomain'. The
--- first Chart corresponds to the 'Domain's chart.
-
-class Manifold m => Map m where
-    type Domain m :: *
-    domain :: m -> Domain m
-    type Codomain m :: *
-    codomain :: m -> Codomain m
-
-class Map m => Apply c d m where
-    -- | 'Map' application.
-    (>.>) :: Function c d :#: m -> c :#: Domain m -> d :#: Codomain m
-    (>.>) f x = head $ f >$> [x]
-    -- | 'Map' list application. May sometimes have a more efficient implementation
-    -- than simply list-mapping (>.>).
-    (>$>) :: Function c d :#: m -> [c :#: Domain m] -> [d :#: Codomain m]
-    (>$>) f = map (f >.>)
-
-infix 8 >.>
-infix 8 >$>
-
-
-
-{-
---- Tables ---
-
-
-newtype Table s = Table s deriving (Eq, Read, Show)
-
-
---- Instances ---
-
-
--- Table --
-
-instance Discrete s => Manifold (Table s) where
-    dimension (Table s) = length $ elements s
-
-instance Discrete s => Function Cartesian (Table s) where
-    type Domain Cartesian (Table s) = s
-    domain cm = let (Table s) = manifold cm in s
-    type Codomain Cartesian (Table s) = Continuum
-    codomain _ = Continuum
-    (>.>) cm k =
-        let ctgs = listCoordinates cm
-            Just (ctg,_) = find ((==k) . snd) . zip ctgs . elements $ domain cm
-         in ctg
-    (>$>) cm ks = (cm >.>) <$> ks
--}
+-- | A 'Manifold' is a 'Map' if it is a binary type-function of two `Manifold's, and can transforms 'Point's on the first 'Manifold' into 'Point's on the second 'Manifold'.
+class (Manifold x, Manifold y, Manifold (f y x)) => Map c f y x where
+    -- | 'Map' application restricted.
+    (>.>) :: c # f y x -> c #* x -> c # y
+    -- | 'Map' vector application. May sometimes have a more efficient implementation
+    -- than simply mapping (>.>).
+    (>$>) :: c # f y x
+          -> [c #* x]
+          -> [c # y]
diff --git a/Goal/Geometry/Map/Linear.hs b/Goal/Geometry/Map/Linear.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Geometry/Map/Linear.hs
@@ -0,0 +1,215 @@
+{-# OPTIONS_GHC -fplugin=GHC.TypeLits.KnownNat.Solver -fplugin=GHC.TypeLits.Normalise -fconstraint-solver-iterations=10 #-}
+{-# LANGUAGE UndecidableInstances,UndecidableSuperClasses #-}
+-- | This module provides tools for working with linear and affine
+-- transformations.
+
+module Goal.Geometry.Map.Linear
+    ( -- * Bilinear Forms
+    Bilinear ((>$<),(>.<),transpose)
+    , (<.<)
+    , (<$<)
+    -- * Tensors
+    , Tensor
+    -- ** Matrix Construction
+    , toMatrix
+    , fromMatrix
+    , toRows
+    , toColumns
+    , fromRows
+    , fromColumns
+    -- ** Computation
+    --, (<#>)
+    , inverse
+    , determinant
+    -- * Affine Functions
+    , Affine (Affine)
+    , Translation ((>+>),anchor)
+    , (>.+>)
+    , (>$+>)
+    , type (<*)
+    ) where
+
+--- Imports ---
+
+-- Package --
+
+import Goal.Core
+
+import Goal.Geometry.Manifold
+import Goal.Geometry.Vector
+import Goal.Geometry.Map
+
+import qualified Goal.Core.Vector.Storable as S
+import qualified Goal.Core.Vector.Generic as G
+
+
+-- Bilinear Forms --
+
+
+-- | A 'Manifold' is 'Bilinear' if its elements are bilinear forms.
+class (Bilinear f x y, Manifold x, Manifold y, Manifold (f x y)) => Bilinear f y x where
+    -- | Tensor outer product.
+    (>.<) :: c # y -> c # x -> c # f y x
+    -- | Average of tensor outer products.
+    (>$<) :: [c # y] -> [c # x] -> c # f y x
+    -- | Tensor transpose.
+    transpose :: c # f y x -> c # f x y
+
+-- | Transposed application.
+(<.<) :: (Map c f x y, Bilinear f y x) => c #* y -> c # f y x -> c # x
+{-# INLINE (<.<) #-}
+(<.<) dy f = transpose f >.> dy
+
+-- | Mapped transposed application.
+(<$<) :: (Map c f x y, Bilinear f y x) => [c #* y] -> c # f y x -> [c # x]
+{-# INLINE (<$<) #-}
+(<$<) dy f = transpose f >$> dy
+
+
+-- Tensor Products --
+
+-- | 'Manifold' of 'Tensor's given by the tensor product of the underlying pair of 'Manifold's.
+data Tensor y x
+
+-- | The inverse of a tensor.
+inverse
+    :: (Manifold x, Manifold y, Dimension x ~ Dimension y)
+    => c # Tensor y x -> c #* Tensor x y
+{-# INLINE inverse #-}
+inverse p = fromMatrix . S.pseudoInverse $ toMatrix p
+
+-- | The determinant of a tensor.
+determinant
+    :: (Manifold x, Manifold y, Dimension x ~ Dimension y)
+    => c # Tensor y x
+    -> Double
+{-# INLINE determinant #-}
+determinant = S.determinant . toMatrix
+
+-- | Converts a point on a 'Tensor manifold into a Matrix.
+toMatrix :: (Manifold x, Manifold y) => c # Tensor y x -> S.Matrix (Dimension y) (Dimension x) Double
+{-# INLINE toMatrix #-}
+toMatrix (Point xs) = G.Matrix xs
+
+-- | Converts a point on a 'Tensor' manifold into a a vector of rows.
+toRows :: (Manifold x, Manifold y) => c # Tensor y x -> S.Vector (Dimension y) (c # x)
+{-# INLINE toRows #-}
+toRows tns = S.map Point . S.toRows $ toMatrix tns
+
+-- | Converts a point on a 'Tensor' manifold into a a vector of rows.
+toColumns :: (Manifold x, Manifold y) => c # Tensor y x -> S.Vector (Dimension x) (c # y)
+{-# INLINE toColumns #-}
+toColumns tns = S.map Point . S.toColumns $ toMatrix tns
+
+-- | Converts a vector of rows into a 'Tensor'.
+fromRows :: (Manifold x, Manifold y) => S.Vector (Dimension y) (c # x) -> c # Tensor y x
+{-# INLINE fromRows #-}
+fromRows rws = fromMatrix . S.fromRows $ S.map coordinates rws
+
+-- | Converts a vector of rows into a 'Tensor'.
+fromColumns :: (Manifold x, Manifold y) => S.Vector (Dimension x) (c # y) -> c # Tensor y x
+{-# INLINE fromColumns #-}
+fromColumns rws = fromMatrix . S.fromColumns $ S.map coordinates rws
+
+-- | Converts a Matrix into a 'Point' on a 'Tensor 'Manifold'.
+fromMatrix :: S.Matrix (Dimension y) (Dimension x) Double -> c # Tensor y x
+{-# INLINE fromMatrix #-}
+fromMatrix (G.Matrix xs) = Point xs
+
+
+--- Affine Functions ---
+
+
+-- | An 'Affine' 'Manifold' represents linear transformations followed by a
+-- translation. The 'First' component is the translation, and the 'Second'
+-- component is the linear transformation.
+newtype Affine f y z x = Affine (z,f y x)
+
+deriving instance (Manifold z, Manifold (f y x)) => Manifold (Affine f y z x)
+deriving instance (Manifold z, Manifold (f y x)) => Product (Affine f y z x)
+
+-- | Infix synonym for simple 'Affine' transformations.
+type (y <* x) = Affine Tensor y y x
+infixr 6 <*
+
+-- | The 'Translation' class is used to define translations where we only want
+-- to translate a subset of the parameters of the given object.
+class (Manifold y, Manifold z) => Translation z y where
+    -- | Translates the the first argument by the second argument.
+    (>+>) :: c # z -> c # y -> c # z
+    -- | Returns the subset of the parameters of the given 'Point' that are
+    -- translated in this instance.
+    anchor :: c # z -> c # y
+
+-- | Operator that applies a 'Map' to a subset of an input's parameters.
+(>.+>) :: (Map c f y x, Translation z x) => c # f y x -> c #* z -> c # y
+(>.+>) f w = f >.> anchor w
+
+-- | Operator that maps a 'Map' over a subset of the parameters of a list of inputs.
+(>$+>) :: (Map c f y x, Translation z x) => c # f y x -> [c #* z] -> [c # y]
+(>$+>) f w = f >$> (anchor <$> w)
+
+
+--- Instances ---
+
+-- Tensors --
+
+instance (Manifold x, Manifold y) => Manifold (Tensor y x) where
+    type Dimension (Tensor y x) = Dimension x * Dimension y
+
+instance (Manifold x, Manifold y) => Map c Tensor y x where
+    {-# INLINE (>.>) #-}
+    (>.>) pq (Point xs) = Point $ S.matrixVectorMultiply (toMatrix pq) xs
+    {-# INLINE (>$>) #-}
+    (>$>) pq qs = Point <$> S.matrixMap (toMatrix pq) (coordinates <$> qs)
+
+instance (Manifold x, Manifold y) => Bilinear Tensor y x where
+    {-# INLINE (>.<) #-}
+    (>.<) (Point pxs) (Point qxs) = fromMatrix $ pxs `S.outerProduct` qxs
+    {-# INLINE (>$<) #-}
+    (>$<) ps qs = fromMatrix . S.averageOuterProduct $ zip (coordinates <$> ps) (coordinates <$> qs)
+    {-# INLINE transpose #-}
+    transpose (Point xs) = fromMatrix . S.transpose $ G.Matrix xs
+
+
+-- Affine Maps --
+
+instance Manifold z => Translation z z where
+    (>+>) z1 z2 = z1 + z2
+    anchor = id
+
+instance (Manifold z, Manifold y) => Translation (y,z) y where
+    (>+>) yz y' =
+        let (y,z) = split yz
+         in join (y + y') z
+    anchor = fst . split
+
+instance (Translation z y, Map c f y x) => Map c (Affine f y) z x where
+    {-# INLINE (>.>) #-}
+    (>.>) fyzx x =
+        let (yz,yx) = split fyzx
+         in   yz >+> (yx >.> x)
+    (>$>) fyzx xs =
+        let (yz,yx) = split fyzx
+         in (yz >+>) <$> yx >$> xs
+
+--instance (KnownNat n, Translation w z)
+--  => Translation (Replicated n w) (Replicated n z) where
+--      {-# INLINE (>+>) #-}
+--      (>+>) w z =
+--          let ws = splitReplicated w
+--              zs = splitReplicated z
+--           in joinReplicated $ S.zipWith (>+>) ws zs
+--      {-# INLINE anchor #-}
+--      anchor = mapReplicatedPoint anchor
+
+
+--instance (Map c f z x) => Map c (Affine f z) z x where
+--    {-# INLINE (>.>) #-}
+--    (>.>) ppq q =
+--        let (p,pq) = split ppq
+--         in p + pq >.> q
+--    {-# INLINE (>$>) #-}
+--    (>$>) ppq qs =
+--        let (p,pq) = split ppq
+--         in (p +) <$> (pq >$> qs)
diff --git a/Goal/Geometry/Map/Linear/Convolutional.hs b/Goal/Geometry/Map/Linear/Convolutional.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Geometry/Map/Linear/Convolutional.hs
@@ -0,0 +1,172 @@
+{-# OPTIONS_GHC -fplugin=GHC.TypeLits.KnownNat.Solver -fplugin=GHC.TypeLits.Normalise -fconstraint-solver-iterations=10 #-}
+{-# LANGUAGE ConstraintKinds,TypeApplications,UndecidableInstances #-}
+
+-- | Manifolds of 'Convolutional' operators. This is hardly used, but could in
+-- theory power conv nets. One day.
+module Goal.Geometry.Map.Linear.Convolutional
+    ( -- * Convolutional Manifolds
+      Convolutional
+    , KnownConvolutional
+    ) where
+
+
+--- Imports ---
+
+
+-- Goal --
+
+import Goal.Core
+import Goal.Geometry.Manifold
+import Goal.Geometry.Map
+import Goal.Geometry.Vector
+import Goal.Geometry.Map.Linear
+import Goal.Geometry.Differential
+
+import qualified Goal.Core.Vector.Generic as G
+import qualified Goal.Core.Vector.Storable as S
+
+
+-- Convolutional Layers --
+
+-- | A 'Manifold' of correlational/convolutional transformations, defined by the
+-- number of kernels, their radius, the depth of the input, and its number of
+-- rows and columns.
+data Convolutional (rd :: Nat) (r :: Nat) (c :: Nat) :: Type -> Type -> Type
+
+-- | A convenience type for ensuring that all the type-level Nats of a
+-- 'Convolutional' 'Manifold's are 'KnownNat's.
+type KnownConvolutional rd r c z x
+  = ( KnownNat rd, KnownNat r, KnownNat c, 1 <= r*c
+    , Dimension x ~ (Div (Dimension x) (r*c) * r*c)
+    , Dimension z ~ (Div (Dimension z) (r*c) * r*c)
+    , Manifold (Convolutional rd r c z x)
+    , Manifold x, Manifold z
+    , KnownNat (Div (Dimension x) (r*c))
+    , KnownNat (Div (Dimension z) (r*c))
+    )
+
+inputToImage
+    :: (KnownConvolutional rd r c z x)
+    => a # Convolutional rd r c z x
+    -> a #* x
+    -> S.Matrix (Div (Dimension x) (r*c)) (r*c) Double
+{-# INLINE inputToImage #-}
+inputToImage _ (Point img) = G.Matrix img
+
+outputToImage
+    :: (KnownConvolutional rd r c z x)
+    => a # Convolutional rd r c z x
+    -> a #* z
+    -> S.Matrix (Div (Dimension z) (r*c)) (r*c) Double
+{-# INLINE outputToImage #-}
+outputToImage _ (Point img) = G.Matrix img
+
+layerToKernels
+    :: ( KnownConvolutional rd r c z x)
+    => a # Convolutional rd r c z x
+    -> S.Matrix (Div (Dimension z) (r*c)) (Div (Dimension x) (r*c) * (2*rd+1)*(2*rd+1)) Double
+{-# INLINE layerToKernels #-}
+layerToKernels (Point krns) = G.Matrix krns
+
+convolveApply
+    :: forall a rd r c z x
+    . KnownConvolutional rd r c z x
+    => a # Convolutional rd r c z x
+    -> a #* x
+    -> a # z
+{-# INLINE convolveApply #-}
+convolveApply cnv imp =
+    let img :: S.Matrix (Div (Dimension x) (r*c)) (r*c) Double
+        img = inputToImage cnv imp
+        krns :: S.Matrix (Div (Dimension z) (r*c)) (Div (Dimension x) (r*c) * (2*rd+1)*(2*rd+1)) Double
+        krns = layerToKernels cnv
+     in Point . G.toVector
+         $ S.crossCorrelate2d (Proxy @ rd) (Proxy @ rd) (Proxy @ r) (Proxy @ c) krns img
+
+convolveTranspose
+    :: forall a rd r c z x
+    . KnownConvolutional rd r c z x
+    => a # Convolutional rd r c z x
+    -> a # Convolutional rd r c x z
+{-# INLINE convolveTranspose #-}
+convolveTranspose cnv =
+    let krns = layerToKernels cnv
+        pnk = Proxy :: Proxy (Div (Dimension z) (r*c))
+        pmd = Proxy :: Proxy (Div (Dimension x) (r*c))
+        krn' :: S.Matrix (Div (Dimension x) (r*c)) (Div (Dimension z) (r*c)*(2*rd+1)*(2*rd+1)) Double
+        krn' = S.kernelTranspose pnk pmd (Proxy @ rd) (Proxy @ rd) krns
+     in Point $ G.toVector krn'
+
+--convolveTransposeApply
+--    :: forall a rd r c z x . KnownConvolutional rd r c z x
+--    => Dual a # z
+--    -> a #> Convolutional rd r c z x
+--    -> a # x
+--{-# INLINE convolveTransposeApply #-}
+--convolveTransposeApply imp cnv =
+--    let img = outputToImage cnv imp
+--        krns = layerToKernels cnv
+--     in Point . G.toVector
+--         $ S.convolve2d (Proxy @ rd) (Proxy @ rd) (Proxy @ r) (Proxy @ c) krns img
+
+convolutionalOuterProduct
+    :: forall a rd r c z x
+    . KnownConvolutional rd r c z x
+      => a # z
+      -> a # x
+      -> a # Convolutional rd r c z x
+{-# INLINE convolutionalOuterProduct #-}
+convolutionalOuterProduct (Point oimg) (Point iimg) =
+    let omtx = G.Matrix oimg
+        imtx = G.Matrix iimg
+     in Point . G.toVector $ S.kernelOuterProduct (Proxy @ rd) (Proxy @ rd) (Proxy @ r) (Proxy @ c) omtx imtx
+
+convolvePropagate
+    :: forall a rd r c z x . KnownConvolutional rd r c z x
+      => [a #* z]
+      -> [a #* x]
+      -> a # Convolutional rd r c z x
+      -> (a #* Convolutional rd r c z x, [a # z])
+{-# INLINE convolvePropagate #-}
+convolvePropagate omps imps cnv =
+    let prdkr = Proxy :: Proxy rd
+        prdkc = Proxy :: Proxy rd
+        pmr = Proxy :: Proxy r
+        pmc = Proxy :: Proxy c
+        foldfun (omp,imp) (k,dkrns) =
+            let img = inputToImage cnv imp
+                dimg = outputToImage cnv omp
+                dkrns' = Point . G.toVector $ S.kernelOuterProduct prdkr prdkc pmr pmc dimg img
+             in (k+1,dkrns' + dkrns)
+     in (uncurry (/>) . foldr foldfun (0,0) $ zip omps imps, cnv >$> imps)
+
+
+--- Instances ---
+
+
+-- Convolutional Manifolds --
+
+instance ( 1 <= r*c, Manifold x, Manifold y, KnownNat r, KnownNat c, KnownNat rd
+         , KnownNat (Div (Dimension x) (r*c)) , KnownNat (Div (Dimension y) (r*c)) )
+  => Manifold (Convolutional rd r c y x) where
+      type Dimension (Convolutional rd r c y x)
+        = (Div (Dimension y) (r * c) * ((Div (Dimension x) (r * c) * (2 * rd + 1)) * (2 * rd + 1)))
+
+
+instance KnownConvolutional rd r c z x => Map a (Convolutional rd r c) z x where
+      {-# INLINE (>.>) #-}
+      (>.>) = convolveApply
+      {-# INLINE (>$>) #-}
+      (>$>) cnv = map (convolveApply cnv)
+
+instance KnownConvolutional rd r c z x => Bilinear (Convolutional rd r c) z x where
+    {-# INLINE (>.<) #-}
+    (>.<) = convolutionalOuterProduct
+    {-# INLINE (>$<) #-}
+    (>$<) ps qs = sum $ zipWith convolutionalOuterProduct ps qs
+    {-# INLINE transpose #-}
+    transpose = convolveTranspose
+
+instance KnownConvolutional rd r c z x => Propagate a (Convolutional rd r c) z x where
+    {-# INLINE propagate #-}
+    propagate = convolvePropagate
diff --git a/Goal/Geometry/Map/Multilinear.hs b/Goal/Geometry/Map/Multilinear.hs
deleted file mode 100644
--- a/Goal/Geometry/Map/Multilinear.hs
+++ /dev/null
@@ -1,211 +0,0 @@
--- | The Map module provides tools for developing function space 'Manifold's.
--- A map is a 'Manifold' where the 'Point's of the Manifold represent
--- parametric functions between 'Manifold's. The defining feature of 'Map's is
--- that they have a particular 'Domain' and 'Codomain', which themselves are
--- 'Manifold's.
-
-module Goal.Geometry.Map.Multilinear (
-    -- * Tensors
-      Tensor (Tensor)
-    -- ** Construction
-    , (>.<)
-    -- ** Matrix Operations
-    , (<#>)
-    , matrixRank
-    , matrixInverse
-    , matrixTranspose
-    , matrixSquareRoot
-    , matrixApply
-    , matrixMap
-    , matrixDiagonalConcatenate
-    -- ** Cartesian
-    , coordinateTransform
-    , linearProjection
-    -- ** HMatrix Conversion
-    , toHMatrix
-    , fromHMatrix
-    -- * Affine Functions
-    , Affine (Affine)
-    , splitAffine
-    , joinAffine
-    ) where
-
---- Imports ---
-
-import Prelude hiding (map,minimum,maximum)
-
--- Package --
-
-import Goal.Core
-
-import Goal.Geometry.Set
-import Goal.Geometry.Manifold
-import Goal.Geometry.Linear
-import Goal.Geometry.Map
-
--- Qualified --
-
-import qualified Data.Vector.Storable as C
-import qualified Numeric.LinearAlgebra.HMatrix as H
-
---import Data.Vector.Storable.UnsafeSerialize
-
-
-
---- Affine Functions ---
-
-
--- | 'Manifold's of 'Affine' functions.
-data Affine m n = Affine m n deriving (Eq, Read, Show)
-
-splitAffine :: (Manifold m, Manifold n) => Function c d :#: Affine m n -> (d :#: m, Function c d :#: Tensor m n)
--- | Splits an 'Point' on an 'Affine' space into a matrix and a constant.
-splitAffine aff =
-    let (Affine m n) = manifold aff
-        tns = Tensor m n
-        css = coordinates aff
-        (mcs,mtxcs) = C.splitAt (dimension m) css
-     in (fromCoordinates m mcs, fromCoordinates tns mtxcs)
-
-joinAffine :: (Manifold m, Manifold n) => d :#: m -> Function c d :#: Tensor m n -> Function c d :#: Affine m n
--- | Combines a matrix and a constant into 'Point' on an 'Affine' space.
-joinAffine dm mtx =
-    let (Tensor m n) = manifold mtx
-     in fromCoordinates (Affine m n) $ coordinates dm C.++ coordinates mtx
-
--- Tensor Products --
-
--- | 'Manifold' of 'Tensor's given by the tensor product of the underlying pair of 'Manifold's.
-data Tensor m n = Tensor m n deriving (Eq, Read, Show)
-
-toHMatrix :: Manifold n => c :#: Tensor m n -> H.Matrix Double
--- | Converts a point on a 'Tensor' product manifold to a matrix for snappy
--- calculation.
-toHMatrix pq =
-    let (Tensor _ m) = manifold pq
-     in H.reshape (dimension m) $ coordinates pq
-
-fromHMatrix :: (Manifold m, Manifold n) => Tensor m n -> H.Matrix Double -> c :#: Tensor m n
-fromHMatrix tns = fromCoordinates tns . H.flatten
-
-matrixRank :: (Manifold m, Manifold n) => c :#: Tensor m n -> Int
-matrixRank = H.rank . toHMatrix
-
-(>.<) :: (Manifold m, Manifold n) => d :#: m -> c :#: n -> Function (Dual c) d :#: Tensor m n
--- | '>.<' denotes the outer product between two points. It provides a way of
--- constructing matrices of the 'Tensor' product space.
-(>.<) p q = fromHMatrix (Tensor (manifold p) $ manifold q) $ coordinates p `H.outer` coordinates q
-
-(<#>) :: (Manifold m, Manifold n, Manifold o)
-      => Function d e :#: Tensor m n -> Function c d :#: Tensor n o -> Function c e :#: Tensor m o
--- | Tensor product composition.
-(<#>) p q =
-    let (Tensor m _) = manifold p
-        (Tensor _ o) = manifold q
-     in fromHMatrix (Tensor m o) $ toHMatrix p <> toHMatrix q
-
-matrixSquareRoot :: Manifold m => c :#: Tensor m m -> c :#: Tensor m m
--- | The square root of a matrix.
-matrixSquareRoot pq = fromHMatrix (manifold pq) . H.sqrtm $ toHMatrix pq
-
-matrixInverse :: (Manifold n, Manifold m) => Function c d :#: Tensor m n -> Function d c :#: Tensor n m
--- | The inverse of a given 'Tensor' point.
-matrixInverse pq =
-    let Tensor m n = manifold pq
-     in fromHMatrix (Tensor n m) . H.inv $ toHMatrix pq
-
-matrixTranspose :: (Manifold m, Manifold n) => Function c d :#: Tensor m n -> Function (Dual d) (Dual c) :#: Tensor n m
--- | The transpose of a given 'Tensor' point.
-matrixTranspose pq =
-    let Tensor m n = manifold pq
-     in fromHMatrix (Tensor n m) . H.tr $ toHMatrix pq
-
-matrixDiagonalConcatenate :: (Manifold m, Manifold n, Manifold o, Manifold p)
-    => Function c d :#: Tensor m n
-    -> Function e f :#: Tensor o p
-    -> Function (c,e) (d,f) :#: Tensor (m,o) (n,p)
--- | Creates a block diagonal matrix.
-matrixDiagonalConcatenate cdmn efop =
-    let (Tensor m n) = manifold cdmn
-        (Tensor o p) = manifold efop
-     in fromHMatrix (Tensor (m,o) (n,p)) $ H.diagBlock [toHMatrix cdmn, toHMatrix efop]
-
-
-coordinateTransform :: Manifold m => [c :#: m] -> Function Cartesian c :#: Tensor m Euclidean
--- | Returns the coordinate transformation from 'Euclidean' space into the space
--- defined by the given basis vectors. This is a glorified fromColumns function.
-coordinateTransform bss =
-    fromHMatrix (Tensor (manifold $ head bss) . Euclidean $ length bss) . H.fromColumns $ coordinates <$> bss
-
-linearProjection :: Manifold m => [Cartesian :#: m] -> Function Cartesian Cartesian :#: Tensor m m
--- | Returns the linear projection operator for the given subset of basis vectors.
-linearProjection bss =
-    let mtx = coordinateTransform bss
-        mtxt = matrixTranspose mtx
-     in mtx <#> matrixInverse (mtxt <#> mtx) <#> mtxt
-
-matrixApply :: (Manifold m, Manifold n) => (Function c d :#: Tensor n m) -> (c :#: m) -> d :#: n
--- | Matrix vector multiplication.
-matrixApply pq p =
-    let (Tensor n _) = manifold pq
-     in fromCoordinates n $ toHMatrix pq H.#> coordinates p
-    {-
-    let (Tensor n m) = manifold pq
-     in if m == manifold p
-          then fromCoordinates n $ toHMatrix pq H.#> coordinates p
-          else error "matrix applied to wrong Manifold"
-          -}
-
-matrixMap :: (Manifold m, Manifold n) => (Function c d :#: Tensor m n) -> [c :#: n] -> [d :#: m]
--- | Mapped matrix vector multiplication, where we first turn the input vectors into a matrix itself (this can greatly improve computation time).
-matrixMap pq ps =
-    let (Tensor n _) = manifold pq
-        mtx = toHMatrix pq
-        xs = H.fromColumns $ coordinates <$> ps
-     in map (fromCoordinates n) . H.toColumns $ mtx <> xs
-    {-
-    let (Tensor n m) = manifold pq
-        mtx = toHMatrix pq
-        xs = H.fromColumns $ coordinates <$> ps
-     in if all (== m) $ manifold <$> ps
-           then map (fromCoordinates n) . H.toColumns $ mtx <> xs
-           else error "matrix applied to wrong Manifold"
-           -}
-
-
---- Instances ---
-
-
--- Tensor Products --
-
-instance (Manifold m, Manifold n) => Manifold (Tensor n m) where
-    dimension (Tensor n m) = dimension m * dimension n
-
-instance (Manifold m, Manifold n) => Map (Tensor m n) where
-    type Domain (Tensor m n) = n
-    domain (Tensor _ n) = n
-    type Codomain (Tensor m n) = m
-    codomain (Tensor m _) = m
-
-instance (Manifold m, Manifold n) => Apply c d (Tensor m n) where
-    (>.>) = matrixApply
-    (>$>) = matrixMap
-
--- Affine Map --
-
-instance (Manifold m, Manifold n) => Manifold (Affine m n) where
-    dimension (Affine m n) = dimension m * dimension n + dimension m
-
-instance (Manifold m, Manifold n) => Map (Affine m n) where
-    type Domain (Affine m n) = n
-    domain (Affine _ n) = n
-    type Codomain (Affine m n) = m
-    codomain (Affine m _) = m
-
-instance (Manifold m, Manifold n) => Apply c d (Affine m n) where
-    (>.>) p x =
-        let (b,mtx) = splitAffine p
-         in mtx >.> x <+> b
-    (>$>) p xs =
-        let (b,mtx) = splitAffine p
-         in map (<+> b) $ mtx >$> xs
diff --git a/Goal/Geometry/Map/NeuralNetwork.hs b/Goal/Geometry/Map/NeuralNetwork.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Geometry/Map/NeuralNetwork.hs
@@ -0,0 +1,122 @@
+{-# OPTIONS_GHC -fplugin=GHC.TypeLits.KnownNat.Solver -fplugin=GHC.TypeLits.Normalise -fconstraint-solver-iterations=10 #-}
+{-# LANGUAGE UndecidableInstances #-}
+
+-- | Multilayer perceptrons which instantiate backpropagation through laziness.
+-- Right now the structure is simplier than it could be, but it leads to nice
+-- types. If anyone ever wants to use a DNN with super-Affine biases, the code
+-- is willing.
+module Goal.Geometry.Map.NeuralNetwork
+    ( -- * Neural Networks
+      NeuralNetwork
+    ) where
+
+
+--- Imports ---
+
+
+-- Goal --
+
+import Goal.Core
+
+import Goal.Geometry.Manifold
+import Goal.Geometry.Map
+import Goal.Geometry.Vector
+import Goal.Geometry.Map.Linear
+import Goal.Geometry.Differential
+
+import qualified Goal.Core.Vector.Storable as S
+
+--- Multilayer ---
+
+
+-- | A multilayer, artificial neural network.
+data NeuralNetwork (gys :: [(Type -> Type -> Type,Type)])
+    (f :: (Type -> Type -> Type)) z x
+
+
+--- Instances ---
+
+
+instance Manifold (Affine f z z x) => Manifold (NeuralNetwork '[] f z x) where
+      type Dimension (NeuralNetwork '[] f z x) = Dimension (Affine f z z x)
+
+instance (Manifold (Affine f z z y), Manifold (NeuralNetwork gys g y x))
+  => Manifold (NeuralNetwork ('(g,y) : gys) f z x) where
+      type Dimension (NeuralNetwork ('(g,y) : gys) f z x)
+        = Dimension (Affine f z z y) + Dimension (NeuralNetwork gys g y x)
+
+
+fromSingleLayerNetwork :: c # NeuralNetwork '[] f z x -> c # Affine f z z x
+{-# INLINE fromSingleLayerNetwork #-}
+fromSingleLayerNetwork = breakPoint
+
+toSingleLayerNetwork :: c # Affine f z z x -> c # NeuralNetwork '[] f z x
+{-# INLINE toSingleLayerNetwork #-}
+toSingleLayerNetwork = breakPoint
+
+-- | Seperates a 'NeuralNetwork' into the final layer and the rest of the network.
+splitNeuralNetwork
+    :: (Manifold (Affine f z z y), Manifold (NeuralNetwork gys g y x))
+    => c # NeuralNetwork ('(g,y):gys) f z x
+    -> (c # Affine f z z y, c # NeuralNetwork gys g y x)
+{-# INLINE splitNeuralNetwork #-}
+splitNeuralNetwork (Point xs) =
+    let (xys,xns) = S.splitAt xs
+     in (Point xys, Point xns)
+
+-- | Joins a layer onto the end of a 'NeuralNetwork'.
+joinNeuralNetwork
+    :: (Manifold (Affine f z z y), Manifold (NeuralNetwork gys g y x))
+    => c # Affine f z z y
+    -> c # NeuralNetwork gys g y x
+    -> c # NeuralNetwork ('(g,y):gys) f z x
+{-# INLINE joinNeuralNetwork #-}
+joinNeuralNetwork (Point xys) (Point xns) =
+    Point $ xys S.++ xns
+
+instance (Manifold (Affine f z z y), Manifold (NeuralNetwork gys g y x))
+  => Product (NeuralNetwork ('(g,y) : gys) f z x) where
+      type First (NeuralNetwork ('(g,y) : gys) f z x)
+        = Affine f z z y
+      type Second (NeuralNetwork ('(g,y) : gys) f z x)
+        = NeuralNetwork gys g y x
+      join = joinNeuralNetwork
+      split = splitNeuralNetwork
+
+instance (Map c f z y, Map c (NeuralNetwork gys g) y x, Transition c (Dual c) y)
+  => Map c (NeuralNetwork ('(g,y) : gys) f) z x where
+    {-# INLINE (>.>) #-}
+    (>.>) fg x =
+        let (f,g) = split fg
+         in f >.> transition (g >.> x)
+    {-# INLINE (>$>) #-}
+    (>$>) fg xs =
+        let (f,g) = split fg
+         in f >$> map transition (g >$> xs)
+
+instance Map c f z x => Map c (NeuralNetwork '[] f) z x where
+    {-# INLINE (>.>) #-}
+    (>.>) f x = fromSingleLayerNetwork f >.> x
+    {-# INLINE (>$>) #-}
+    (>$>) f xs = fromSingleLayerNetwork f >$> xs
+
+instance (Propagate c f z x) => Propagate c (NeuralNetwork '[] f) z x where
+    {-# INLINE propagate #-}
+    propagate dps qs f =
+        let (df,ps) = propagate dps qs $ fromSingleLayerNetwork f
+         in (toSingleLayerNetwork df,ps)
+
+instance
+    ( Propagate c f z y, Propagate c (NeuralNetwork gys g) y x, Map c f y z
+    , Transition c (Dual c) y, Legendre y, Riemannian c y, Bilinear f z y)
+  => Propagate c (NeuralNetwork ('(g,y) : gys) f) z x where
+      {-# INLINE propagate #-}
+      propagate dzs xs fg =
+          let (f,g) = split fg
+              fmtx = snd $ split f
+              mys = transition <$> ys
+              (df,zhts) = propagate dzs mys f
+              (dg,ys) = propagate dys xs g
+              dys0 = dzs <$< fmtx
+              dys = zipWith flat ys dys0
+           in (join df dg, zhts)
diff --git a/Goal/Geometry/Plot.hs b/Goal/Geometry/Plot.hs
deleted file mode 100644
--- a/Goal/Geometry/Plot.hs
+++ /dev/null
@@ -1,29 +0,0 @@
-module Goal.Geometry.Plot where
-
-
---- Imports ---
-
-
--- Goal --
-
-import Goal.Core
-import Goal.Geometry.Set
-
-import qualified Data.Vector.Storable as C
-
-coordinateLogHistogram :: Int -> String -> [String] -> [Coordinates] -> Layout Double Double
-coordinateLogHistogram nbns ttl ttls css =
-    let bplt = plot_bars_titles .~ ttls $ logHistogramPlot0 nbns (C.toList <$> css) def
-     in layout_title .~ ttl
-        $ layout_y_axis . laxis_override .~ axisGridHide
-        $ layout_x_axis . laxis_override .~ axisGridHide
-        $ logHistogramLayout bplt def
-
-coordinateHistogram :: Int -> String -> [String] -> [Coordinates] -> Layout Double Double
-coordinateHistogram nbns ttl ttls css =
-    let bplt = plot_bars_titles .~ ttls
-            $ histogramPlot0 nbns (C.toList <$> css) def
-     in layout_title .~ ttl
-        $ layout_y_axis . laxis_override .~ axisGridHide
-        $ layout_x_axis . laxis_override .~ axisGridHide
-        $ histogramLayout bplt def
diff --git a/Goal/Geometry/Set.hs b/Goal/Geometry/Set.hs
deleted file mode 100644
--- a/Goal/Geometry/Set.hs
+++ /dev/null
@@ -1,134 +0,0 @@
--- | A module for describing 'Set's of 'Element's. Necessary in a few cases (such as discrete sets) that 'Manifold's don't handle well.
-module Goal.Geometry.Set
-    ( -- * Sets
-      Set
-    , Element
-    , Discrete (elements)
-    -- * Instances
-    -- ** Discrete
-    , Boolean (Boolean)
-    , NaturalNumbers (NaturalNumbers)
-    , Integers (Integers)
-    -- ** Continuous
-    , Coordinates
-    , Euclidean (Euclidean)
-    , Continuum (Continuum)
-    -- * Combinators
-    -- ** Replicated
-    , Replicated (Replicated)
-    ) where
-
-
---- Imports ---
-
-
--- Goal --
-
-import Goal.Core
-
--- Qualified --
-
-import qualified Data.Vector.Storable as C
-
-
---- Classes ---
-
-
--- | 'Set's are collections of distinguishable 'Element's.
-class (Eq s, Eq (Element s)) => Set s where
-    type Element s :: *
-
-
--- | A 'Discrete' 'Set' is one where we can list its elements. The
--- returned list should satisfy the law
---
--- > elements s = nub $ elements s
---
-class Set s => Discrete s where
-    elements :: s -> [Element s]
-
-
---- Types ---
-
-
--- Discrete --
-
--- | The set of natural numbers.
-data NaturalNumbers = NaturalNumbers deriving (Eq,Read,Show)
-
--- | The set of integers.
-data Integers = Integers deriving (Eq,Read,Show)
-
--- | 'True' and 'False'.
-data Boolean = Boolean deriving (Eq,Read,Show)
-
--- Continuous  --
-
--- | 'Euclidean' space.
-newtype Euclidean = Euclidean Int deriving (Eq,Read,Show)
-
--- | One dimensional 'Euclidean' space.
-data Continuum = Continuum deriving (Eq,Read,Show)
-
--- | 'Element's of 'Euclidean' spaces are referred to as 'Coordinates'.
-type Coordinates = C.Vector Double
-
--- Replicated --
-
--- | A 'Replicated' set is a single set multiplied a specified number of times
--- via the Cartesian product.
-data Replicated m = Replicated !m !Int deriving (Eq,Read,Show)
-
-
---- Instances ---
-
-
--- Discrete --
-
-instance Set NaturalNumbers where
-    type Element NaturalNumbers = Int
-
-instance Discrete NaturalNumbers where
-    elements _ = [0..]
-
-instance Set Integers where
-    type Element Integers = Int
-
-instance Discrete Integers where
-    elements _ = (0:) $ concat [ [-k,k] | k <- [1..] ]
-
-instance Set Boolean where
-    type Element Boolean = Bool
-
-instance Discrete Boolean where
-    elements _ = [True,False]
-
-instance Eq k => Set [k] where
-    type Element [k] = k
-
-instance Eq k => Discrete [k] where
-    elements = id
-
--- Continuous --
-
-instance Set Continuum where
-    type Element Continuum = Double
-
-instance Set Euclidean where
-    type Element Euclidean = Coordinates
-
-
--- Replicated --
-
-instance Set s => Set (Replicated s) where
-    type Element (Replicated s) = [Element s]
-
-
-instance Discrete s => Discrete (Replicated s) where
-    elements (Replicated s n) = replicateM n $ elements s
-
--- Direct Sums --
-
-instance (Set s, Set r) => Set (s,r) where
-    type Element (s,r) = (Element s,Element r)
-
diff --git a/Goal/Geometry/Vector.hs b/Goal/Geometry/Vector.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Geometry/Vector.hs
@@ -0,0 +1,72 @@
+{-# LANGUAGE UndecidableSuperClasses #-}
+-- | The Linear module provides the tools for treating a locally Euclidean patch
+-- of a manifold as a linear space.
+module Goal.Geometry.Vector
+    ( -- * Vector Spaces
+      (.>)
+    , (/>)
+    , convexCombination
+    -- * Dual Spaces
+    , Primal (Dual)
+    , type (#*)
+    , (<.>)
+    , dotMap
+    ) where
+
+--- Imports ---
+
+-- Package --
+
+import Goal.Core
+import Goal.Geometry.Manifold
+
+import qualified Goal.Core.Vector.Storable as S
+
+--- Vector Spaces on Manifolds ---
+
+
+-- | Scalar multiplication of points on a manifold.
+(.>) :: Double -> c # x -> c # x
+{-# INLINE (.>) #-}
+(.>) a (Point xs) = Point $ S.scale a xs
+infix 7 .>
+
+-- | Scalar division of points on a manifold.
+(/>) :: Double -> c # x -> c # x
+{-# INLINE (/>) #-}
+(/>) a (Point xs) = Point $ S.scale (recip a) xs
+infix 7 />
+
+-- | Combination of two 'Point's. Takes the first argument of the second
+-- argument, and (1-first argument) of the third argument.
+convexCombination :: Manifold x => Double -> c # x -> c # x -> c # x
+convexCombination x p1 p2 = x .> p1 + (1-x) .> p2
+
+
+--- Dual Spaces ---
+
+
+-- | 'Primal' charts have a 'Dual' coordinate system.
+class (Dual (Dual c) ~ c, Primal (Dual c)) => Primal c where
+    type Dual c :: Type
+
+-- | A 'Point' on a 'Manifold' in the 'Dual' coordinates of c.
+type (c #* x) = Point (Dual c) x
+infix 3 #*
+
+-- | '<.>' is the inner product between a dual pair of 'Point's.
+(<.>) :: c # x -> c #* x -> Double
+{-# INLINE (<.>) #-}
+(<.>) p q = S.dotProduct (coordinates p) (coordinates q)
+
+infix 7 <.>
+
+-- | 'dotMap' computes the inner product over a list of dual elements.
+dotMap :: Manifold x => c # x -> [c #* x] -> [Double]
+{-# INLINE dotMap #-}
+dotMap p qs = S.dotMap (coordinates p) (coordinates <$> qs)
+
+-- Cartesian Spaces --
+
+instance Primal Cartesian where
+    type Dual Cartesian = Cartesian
diff --git a/LICENSE b/LICENSE
--- a/LICENSE
+++ b/LICENSE
@@ -1,4 +1,4 @@
-Copyright (c) 2014, Sacha Sokoloski
+Copyright (c) 2017, Sacha Sokoloski
 
 All rights reserved.
 
diff --git a/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,95 @@
+In this package we find all the basic types and classes which drive the
+manifold/geometry based approach of Goal. In particular, points and manifolds,
+dual spaces, function spaces and multilayer neural networks, and generic
+optimization routines such as gradient pursuit. What follows is very brief
+introduction to how we define points on a manifold in Goal.
+
+The fundamental class in Goal is the `Manifold`
+```haskell
+class KnownNat (Dimension x) => Manifold x where
+    type Dimension x :: Nat
+```
+`Manifold`s have an associated type, which is the `Dimension` of the `Manifold`.
+The `Dimension` of a `Manifold` tells us the size required of vector to
+represent a 'Point's on the given `Manifold`. In turn a `Point` is defined as
+```haskell
+newtype Point c x =
+    Point { coordinates :: S.Vector (Dimension x) Double }
+```
+At the value level, a `Point` is a wrapper around an `S.Vector`, which is a
+storable vector from the
+[vector-sized](https://hackage.haskell.org/package/vector-sized) package, with
+size `Dimension x`. In general, numerical operations in Goal are defined in
+terms of [vector-sized](https://hackage.haskell.org/package/vector-sized) and
+[hmatrix](https://hackage.haskell.org/package/hmatrix), with specific functions
+for applying operations in bulk. Although I make no promises, Goal should be
+quite efficient, at least for a CPU-based numerical library.
+
+To continue, a `Point` is defined at the type-level by a `Manifold` `x`, and the
+mysterious phantom type `c`.  In Goal `c` is referred to as a coordinate system,
+or more succinctly as a chart.  A coordinate system describes how the abstract
+elements of a `Manifold` may be uniquely represented by a vector of numbers. In
+Goal we usually refer to `Point`s with the following infix type synonym
+```haskell
+type (c # x) = Point c x
+```
+which we may read as a `Point` in `c` coordinates on the `x` `Manifold`. I chose
+the `#` symbol because it is reminiscent of the grid of a coordinate system.
+
+Finally, with the notion of a coordinate system in hand, we may definition
+`transition` functions for re-representing `Point`s in alternative coordinate
+systems
+```haskell
+class Transition c d x where
+    transition :: c # x -> d # x
+```
+
+As an example, where we define `Euclidean` space
+```haskell
+data Euclidean (n :: Nat)
+
+instance (KnownNat n) => Manifold (Euclidean n) where
+    type Dimension (Euclidean n) = n
+```
+and two coordinate systems on Euclidean space with an appropriate transition function
+```haskell
+data Cartesian
+data Polar
+
+instance Transition Cartesian Polar (Euclidean 2) where
+    {-# INLINE transition #-}
+    transition p =
+        let [x,y] = listCoordinates p
+            r = sqrt $ (x*x) + (y*y)
+            phi = atan2 y x
+         in fromTuple (r,phi)
+```
+we may create a `Point` in `Cartesian` coordinates an easily convert it to `Polar` coordinates
+```haskell
+xcrt :: Cartesian # Euclidean 2
+xcrt = fromTuple (1,2)
+
+xplr :: Polar # Euclidean 2
+xplr = transition xcrt
+```
+
+So what has this bought us? Why would we make use of not only one, but
+essentially two phantom types for describing vectors? Intuitively, the
+`Manifold` under investigation is what we care about. If, for example, we
+consider a `Manifold` of probability distributions, it is the distributions
+themselves we care about. But distributions are abstract things, and so we
+represent them in various coordinate systems (e.g. mean and variance) to handle
+them numerically.
+
+The charts available for a given `Manifold` are thus different (but isomorphic)
+representations of the same thing. In particular, many coordinate systems have a
+dual coordinate system that describes function differentials, which is critical
+for numerical optimization. In general, many optimization problems can be
+greatly simplified by finding the right coordinate system, and many complex
+optimization problems can be solved by sequence of coordinate transformations
+and simple numerical operations. Numerically the resulting computation is not
+trivial, but theoretically it becomes an intuitive thing.
+
+For in-depth tutorials visit my
+[blog](https://sacha-sokoloski.gitlab.io/website/pages/blog.html).
+
diff --git a/goal-geometry.cabal b/goal-geometry.cabal
--- a/goal-geometry.cabal
+++ b/goal-geometry.cabal
@@ -1,59 +1,48 @@
+cabal-version: 3.0
 name: goal-geometry
-version: 0.1
-synopsis: Scientific computing on geometric objects
-description: This library provides all the types and classes essential for
-    defining manifolds. This includes definitions and algorithms for sets,
-    points, linear and multilinear algebra, function spaces, basic differential
-    geometry, and convex analysis.
-license: BSD3
+version: 0.20
+synopsis: The basic geometric type system of Goal
+description: goal-geometry provides the basic types and classes which drive the manifold/geometry based approach of Goal. Points and manifolds, dual spaces, function spaces and multilayer neural networks, and generic optimization routines are defined here.
+license: BSD-3-Clause
 license-file: LICENSE
+extra-source-files: README.md
 author: Sacha Sokoloski
-maintainer: sokolo@mis.mpg.de
+maintainer: sacha.sokoloski@mailbox.org
+homepage: https://gitlab.com/sacha-sokoloski/goal
 category: Math
 build-type: Simple
-cabal-version: >=1.10
 
 library
     exposed-modules:
         Goal.Geometry,
-        Goal.Geometry.Set,
         Goal.Geometry.Manifold,
-        Goal.Geometry.Linear,
+        Goal.Geometry.Vector,
         Goal.Geometry.Map,
-        Goal.Geometry.Map.Multilinear,
+        Goal.Geometry.Map.Linear,
+        Goal.Geometry.Map.Linear.Convolutional,
+        Goal.Geometry.Map.NeuralNetwork,
         Goal.Geometry.Differential,
-        Goal.Geometry.Differential.Convex,
-        Goal.Geometry.Plot
-    build-depends:
-        base==4.*,
-        goal-core==0.1,
-        vector==0.11.*,
-        hmatrix==0.17.*
-    default-extensions: TypeOperators, TypeFamilies, MultiParamTypeClasses,
-        FlexibleInstances, FlexibleContexts
-    default-language: Haskell2010
-    ghc-options: -O2 -Wall -fno-warn-type-defaults -fno-warn-missing-signatures
-
-executable coordinates
-    main-is: coordinates.hs
-    hs-source-dirs: scripts
-    ghc-options: -Wall -O2 -threaded -rtsopts -fno-warn-type-defaults
-        -fno-warn-missing-signatures -fno-warn-unused-do-bind
-    build-depends:
-        base==4.*,
-        goal-core==0.1,
-        goal-geometry==0.1
-    default-language: Haskell2010
-
-executable gradient-descent
-    main-is: gradient-descent.hs
-    hs-source-dirs: scripts
-    ghc-options: -Wall -O2 -threaded -rtsopts -fno-warn-type-defaults
-        -fno-warn-missing-signatures -fno-warn-unused-do-bind
+        Goal.Geometry.Differential.GradientPursuit
     build-depends:
-        base==4.*,
-        goal-core==0.1,
-        goal-geometry==0.1
+        base >= 4.13 && < 4.15,
+        goal-core,
+        ad,
+        indexed-list-literals,
+        ghc-typelits-natnormalise,
+        ghc-typelits-knownnat
     default-language: Haskell2010
-
-
+    default-extensions:
+        ScopedTypeVariables,
+        ExplicitNamespaces,
+        TypeOperators,
+        KindSignatures,
+        DataKinds,
+        RankNTypes,
+        TypeFamilies,
+        NoStarIsType,
+        FlexibleContexts,
+        MultiParamTypeClasses,
+        GeneralizedNewtypeDeriving,
+        StandaloneDeriving,
+        FlexibleInstances
+    ghc-options: -Wall -O2
diff --git a/scripts/coordinates.hs b/scripts/coordinates.hs
deleted file mode 100644
--- a/scripts/coordinates.hs
+++ /dev/null
@@ -1,52 +0,0 @@
---- Imports ---
-
-
-import Goal.Core
-import Goal.Geometry
-
-
---- Program ---
-
-
--- Globals --
-
-mxx = pi
-mnx = -pi
-mxy = pi
-mny = -pi
-nstps = 10
-npnts = 50
-
-axprms = LinearAxisParams (show . round) 5 5
-
-hlns = [ [ [x,y] | x <- range mnx mxx npnts ] | y <- range mnx mxx nstps ]
-vlns = [ [ [x,y] | y <- range mny mxy npnts ] | x <- range mny mxy nstps ]
-lns0 = hlns ++ vlns
-
-eclds = map euclideanPoint <$> lns0
-plrs = map (chart Cartesian . transition . chart Polar . fromList (Euclidean 2)) <$> lns0
-
-layoutMaker lns = execEC $ do
-
-    layout_x_axis . laxis_override .= axisGridHide
-    layout_x_axis . laxis_generate .= scaledAxis axprms (-2,2)
-    layout_x_axis . laxis_title .= "x"
-
-    layout_y_axis . laxis_override .= axisGridHide
-    layout_y_axis . laxis_generate .= scaledAxis axprms (-2,2)
-    layout_y_axis . laxis_title .= "y"
-
-    plot . liftEC $ do
-
-        plot_lines_values .= (map toPair <$> lns)
-        plot_lines_style .= solidLine 3 (opaque black)
-
-
--- Main --
-
-main = do
-    let lyt1 = layoutMaker eclds
-        lyt2 = layoutMaker plrs
-        rnbl = toRenderable $ StackedLayouts [StackedLayout lyt1, StackedLayout lyt2] False
-    --renderableToAspectWindow False 400 800 rnbl
-    void $ renderableToFile (FileOptions (200,400) PDF) "coordinates.pdf" rnbl
diff --git a/scripts/gradient-descent.hs b/scripts/gradient-descent.hs
deleted file mode 100644
--- a/scripts/gradient-descent.hs
+++ /dev/null
@@ -1,75 +0,0 @@
---- Imports ---
-
-
-import Goal.Core
-import Goal.Geometry
-
-
---- Globals ---
-
-
--- Functions --
-
-f p = let (x,y) = toPair p in x^2 + 2*y^2 + (x-y)^2
-
-df p =
-    let (x,y) = toPair p
-        x' = 2*x + 2*(x-y)
-        y' = 4*y - 2*(x-y)
-     in fromList (Tangent p) [x',y']
-
--- Plot --
-
-res = 400
-mn = -4
-mx = 4
-niso = 10
-cntrf x y = f $ euclideanPoint [x,y]
-rng = (mn,mx,res)
-clrs = rgbaGradient (0.9,0,0,1) (0,0,0,1) niso
-axprms = LinearAxisParams (show . round) 5 5
-
--- Gradient Descent --
-
-p0 = euclideanPoint [-4,2]
-eps = 0.01
-nstps = 200
-grds = take nstps $ gradientDescent eps df p0
-
-
---- Main ---
-
-
-main = do
-
-    -- Contour plots
-    let rnbl = toRenderable . execEC $ do
-
-            let cntrs = contours rng rng niso cntrf
-
-            sequence_ $ do
-
-                ((_,cntr),clr) <- zip cntrs clrs
-
-                return . plot . liftEC $ do
-
-                    plot_lines_style .= solidLine 3 clr
-                    plot_lines_values .= cntr
-
-            layout_x_axis . laxis_generate .= scaledAxis axprms (mn,mx)
-            layout_x_axis . laxis_override .= axisGridHide
-            layout_x_axis . laxis_title .= "x"
-            layout_y_axis . laxis_generate .= scaledAxis axprms (mn,mx)
-            layout_y_axis . laxis_override .= axisGridHide
-            layout_y_axis . laxis_title .= "y"
-
-            plot . liftEC $ do
-                plot_points_style .= filledCircles 5 (opaque red)
-                plot_points_values .= [(0,0)]
-
-            plot . liftEC $ do
-                plot_lines_style .= solidLine 3 (opaque black)
-                plot_lines_values .= [toPair <$> grds]
-
-    --void $ renderableToAspectWindow False 800 800 rnbl
-    void $ renderableToFile (FileOptions (200,200) PDF) "gradient-descent.pdf" rnbl