diff --git a/Goal/Geometry.hs b/Goal/Geometry.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Geometry.hs
@@ -0,0 +1,27 @@
+module Goal.Geometry
+    (
+    -- * Re-Exports
+      module Goal.Geometry.Set
+    , module Goal.Geometry.Manifold
+    , module Goal.Geometry.Linear
+    , module Goal.Geometry.Map
+    , module Goal.Geometry.Map.Multilinear
+    , module Goal.Geometry.Differential
+    , module Goal.Geometry.Differential.Convex
+    , module Goal.Geometry.Plot
+    ) where
+
+
+-- Imports --
+
+
+-- Re-exports --
+
+import Goal.Geometry.Set
+import Goal.Geometry.Manifold
+import Goal.Geometry.Linear
+import Goal.Geometry.Map
+import Goal.Geometry.Map.Multilinear
+import Goal.Geometry.Differential
+import Goal.Geometry.Differential.Convex
+import Goal.Geometry.Plot
diff --git a/Goal/Geometry/Differential.hs b/Goal/Geometry/Differential.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Geometry/Differential.hs
@@ -0,0 +1,236 @@
+-- | 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
+    ) 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
+import Goal.Geometry.Map.Multilinear
+
+-- Qualified --
+
+import qualified Data.Vector.Storable as C
+import qualified Numeric.LinearAlgebra.HMatrix as H
+
+--import Data.Vector.Storable.UnsafeSerialize
+
+
+--- Differentiable Manifolds ---
+
+
+-- | '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)
+
+-- | A 'Tangent' 'Bundle' is the original 'Manifold' combined with all its
+-- 'Tangent' spaces.
+newtype Bundle c m = Bundle { removeBundle :: m } deriving (Eq, Read, Show)
+
+-- | 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)
+
+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
+
+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
+
+
+-- 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)
+
+
+--- Riemannian Manifolds ---
+
+
+-- | '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
+
+
+--- 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
+
+-- Tangent Spaces --
+
+instance Manifold m => Manifold (Tangent c m) where
+    dimension (Tangent p) = dimension $ manifold p
+
+instance Manifold m => Manifold (Bundle c m) where
+    dimension (Bundle m) = 2 * dimension m
+
+-- Tanget Space Coordinates --
+
+instance Primal Partials where
+    type Dual Partials = Differentials
+
+instance Primal Differentials where
+    type Dual Differentials = Partials
+
+
+--- Graveyard ---
+
+
+
+{-
+--- Functions ---
+
+
+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
+
+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
+-}
diff --git a/Goal/Geometry/Differential/Convex.hs b/Goal/Geometry/Differential/Convex.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Geometry/Differential/Convex.hs
@@ -0,0 +1,52 @@
+-- | 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/Linear.hs b/Goal/Geometry/Linear.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Geometry/Linear.hs
@@ -0,0 +1,108 @@
+-- | 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
new file mode 100644
--- /dev/null
+++ b/Goal/Geometry/Manifold.hs
@@ -0,0 +1,284 @@
+-- | 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'.
+module Goal.Geometry.Manifold
+    ( -- * Manifolds
+      Manifold (dimension)
+    , Transition (transition)
+    -- ** Sets
+    , Embedded (Embedded, disembed)
+    -- ** Points
+    , Coordinates
+    , (:#:) (coordinates, manifold)
+    , coordinate
+    , chart
+    , breakChart
+    , alterChart
+    , listCoordinates
+    , alterCoordinates
+    , toPair
+    -- ** Charts
+    , Cartesian (Cartesian)
+    , Polar (Polar)
+    -- ** Constructors
+    , fromList
+    , fromCoordinates
+    , euclideanPoint
+    , realNumber
+    -- * Direct Sums
+    -- ** Replicated
+    , mapReplicated
+    , joinReplicated
+    , concatReplicated
+    -- ** DirectSum
+    , joinPair
+    , splitPair
+    , joinPair'
+    , splitPair'
+    , joinTriple
+    , splitTriple
+    , joinTriple'
+    , splitTriple'
+    ) where
+
+
+--- Imports ---
+
+
+-- Goal --
+
+import Goal.Core
+
+import Goal.Geometry.Set
+
+-- Qualified --
+
+import qualified Data.Vector.Storable as C
+
+
+--- 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 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)
+
+infixr 1 :#:
+
+coordinate :: Int -> c :#: m -> Double
+coordinate n (Point cs _) = cs C.! n
+
+data Embedded m c = Embedded { disembed :: m } deriving (Eq, Read, Show)
+
+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)
+
+alterChart :: Manifold m => d -> c :#: m -> d :#: m
+-- | Combines 'breakChart' and 'chart'.
+alterChart _ = breakChart
+
+toPair :: c :#: m -> (Double,Double)
+toPair p = (coordinate 0 p,coordinate 1 p)
+
+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
+
+listCoordinates :: c :#: m -> [Double]
+-- | Returns the 'Coordinates' of the point in list form.
+listCoordinates (Point cs _) = C.toList cs
+
+-- | 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
+
+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
+
+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) ++ ")."
+
+euclideanPoint :: [Double] -> Cartesian :#: Euclidean
+-- | A convenience function for building 'Euclidean' vectors.
+euclideanPoint xs = fromList (Euclidean $ length xs) xs
+
+realNumber :: Double -> Cartesian :#: Continuum
+-- | A convenience function for building elements of a 'Continuum'.
+realNumber x = fromList Continuum [x]
+
+--- Construction ---
+
+
+-- Euclidean --
+
+-- | The 'Cartesian' coordinate system.
+data Cartesian = Cartesian
+
+-- | The 'Polar' coordinate system.
+data Polar = Polar
+
+-- | 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 ] ]
+
+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)
+
+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)
+
+-- Direct Sums --
+
+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
+
+splitPair :: (Manifold m, Manifold n) => (c,d) :#: (m,n) -> (c :#: m, d :#: n)
+-- | Splits a direct sum pair.
+splitPair = unsafeSplitPair
+
+joinPair' :: (Manifold m, Manifold n) => c :#: m -> c :#: n -> c :#: (m,n)
+-- | Alternative version where we assume that the Charts are shared.
+joinPair' = unsafeJoinPair
+
+splitPair' :: (Manifold m, Manifold n) => c :#: (m,n) -> (c :#: m, c :#: n)
+-- | Alternative version where we assume that the Charts are shared.
+splitPair' = unsafeSplitPair
+
+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
+
+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)
+
+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
+
+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
+
+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
+
+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)
+
+
+--- Instances ---
+
+
+instance Transition c c m where
+    transition = id
+
+-- Embedded --
+
+instance Manifold m => Set (Embedded m c) where
+    type Element (Embedded m c) = c :#: m
+
+-- Euclidean --
+
+instance Manifold Euclidean where
+    dimension (Euclidean n) = n
+
+instance Manifold Continuum where
+    dimension _ = 1
+
+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 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])
+
+-- DirectSum --
+
+instance (Manifold m, Manifold n) => Manifold (m,n) where
+    dimension (m,n) = dimension m + dimension n
+
+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 m => Manifold (Replicated m) where
+    dimension (Replicated m rn) = dimension m * rn
diff --git a/Goal/Geometry/Map.hs b/Goal/Geometry/Map.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Geometry/Map.hs
@@ -0,0 +1,76 @@
+-- | 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 (
+    -- * Maps
+      Map (Domain, domain, Codomain, codomain)
+    , Apply ((>.>), (>$>))
+    -- * Map Charts
+    , Function (Function)
+    ) where
+
+
+--- Imports ---
+
+
+-- Goal --
+
+import Goal.Geometry.Manifold
+
+--- Maps between Manifolds ---
+
+-- 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
+-}
diff --git a/Goal/Geometry/Map/Multilinear.hs b/Goal/Geometry/Map/Multilinear.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Geometry/Map/Multilinear.hs
@@ -0,0 +1,211 @@
+-- | 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/Plot.hs b/Goal/Geometry/Plot.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Geometry/Plot.hs
@@ -0,0 +1,29 @@
+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
new file mode 100644
--- /dev/null
+++ b/Goal/Geometry/Set.hs
@@ -0,0 +1,134 @@
+-- | 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/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright (c) 2014, Sacha Sokoloski
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above
+      copyright notice, this list of conditions and the following
+      disclaimer in the documentation and/or other materials provided
+      with the distribution.
+
+    * Neither the name of Sacha Sokoloski nor the names of other
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/goal-geometry.cabal b/goal-geometry.cabal
new file mode 100644
--- /dev/null
+++ b/goal-geometry.cabal
@@ -0,0 +1,59 @@
+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
+license-file: LICENSE
+author: Sacha Sokoloski
+maintainer: sokolo@mis.mpg.de
+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.Map,
+        Goal.Geometry.Map.Multilinear,
+        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
+    build-depends:
+        base==4.*,
+        goal-core==0.1,
+        goal-geometry==0.1
+    default-language: Haskell2010
+
+
diff --git a/scripts/coordinates.hs b/scripts/coordinates.hs
new file mode 100644
--- /dev/null
+++ b/scripts/coordinates.hs
@@ -0,0 +1,52 @@
+--- 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
new file mode 100644
--- /dev/null
+++ b/scripts/gradient-descent.hs
@@ -0,0 +1,75 @@
+--- 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