diff --git a/Goal/Core.hs b/Goal/Core.hs
--- a/Goal/Core.hs
+++ b/Goal/Core.hs
@@ -1,39 +1,43 @@
+-- | This module re-exports a number of (compatible) modules from across base
+-- and other libraries, as well as most of the modules in @goal-core@. It does
+-- not re-export the Vector modules, which should be imported with
+-- qualification.
 module Goal.Core
     ( -- * Module Exports
-      module Goal.Core.Plot
+      module Goal.Core.Util
+    , module Goal.Core.Project
+    , module Goal.Core.Circuit
     , module Data.Function
+    , module Data.Functor
+    , module Data.Foldable
+    , module Data.Traversable
     , module Data.Ord
-    , module Data.Monoid
-    , module Data.List
     , module Data.Maybe
     , module Data.Either
-    , module Data.Default.Class
+    , module Data.Finite
+    , module Data.Csv
+    , module Data.Proxy
+    , module Data.Kind
+    , module Data.Functor.Identity
+    , module Data.Type.Equality
     , module Control.Applicative
     , module Control.Monad
+    , module Control.Monad.Primitive
     , module Control.Monad.ST
     , module Control.Arrow
-    , module Control.Lens.Type
-    , module Control.Lens.Getter
-    , module Control.Lens.Setter
-    , module Control.Lens.TH
     , module Control.Concurrent
+    , module Control.DeepSeq
     , module Numeric
+    , module Numeric.SpecFunctions
+    , module Options.Applicative
+    , module GHC.TypeNats
+    , module GHC.Generics
     , module Debug.Trace
-    -- * Lists
-    , takeEvery
-    , breakEvery
-    -- * Low-Level
-    , traceGiven
-    -- * Numeric
-    , roundSD
-    , toPi
-    -- ** Functions
-    , logistic
-    , logit
-    -- ** Lists
-    , mean
-    , range
-    , discretizeFunction
+    , module System.Directory
+    -- * (Re)names
+    , NatNumber
+    , ByteString
+    , orderedHeader
     ) where
 
 
@@ -42,140 +46,48 @@
 
 -- Re-exports --
 
-import Goal.Core.Plot hiding (empty,over)
+import Goal.Core.Util
+import Goal.Core.Project
+import Goal.Core.Circuit
 
+import Data.Csv hiding (Parser,Field,header)
+import qualified Data.Csv as CSV
+import Data.Functor
+import Data.Foldable
+import Data.Traversable
 import Data.Ord
-import Data.Function
-import Data.Monoid hiding (Dual)
-import Data.List hiding (sum)
+import Data.Function hiding ((&))
 import Data.Maybe
 import Data.Either
+import Data.Proxy
+import Data.Finite
+import Data.Kind (Type)
+import Data.Functor.Identity
+import Data.Type.Equality
 
-import Control.Applicative
+import Control.Applicative hiding (empty)
 import Control.Arrow hiding ((<+>))
-import Control.Monad
+import Control.Monad hiding (join)
 import Control.Monad.ST
-import Control.Lens.Type
-import Control.Lens.Getter
-import Control.Lens.Setter hiding (Identity)
-import Control.Lens.TH
 import Control.Concurrent
-
-import Debug.Trace
-import Data.Default.Class
-import Numeric
-
-
---- General Functions ---
-
-
-takeEvery :: Int -> [x] -> [x]
--- | Takes every nth element, starting with the head of the list.
-takeEvery m = map snd . filter (\(x,_) -> mod x m == 0) . zip [0..]
-
-breakEvery :: Int -> [x] -> [[x]]
--- | Break the list up into lists of length n.
-breakEvery _ [] = []
-breakEvery n xs = take n xs : breakEvery n (drop n xs)
-
-traceGiven :: Show a => a -> a
--- | Runs traceShow on the given element.
-traceGiven a = traceShow a a
-
-
---- Numeric ---
-
-
-roundSD :: (Floating x, RealFrac x) => Int -> x -> x
--- | Roundest the number to the specified significant digit.
-roundSD n x = (/10^n) . fromIntegral . round $ 10^n * x
-
-toPi :: (Floating x, RealFrac x) => x -> x
--- | Modulo pi thingy.
-toPi x =
-    let xpi = x / pi
-        n = floor xpi
-        f = xpi - fromIntegral n
-    in if even n then pi * f else -(pi * (1 - f))
-
-logistic :: Floating x => x -> x
--- | A standard sigmoid function.
-logistic x = 1 / (1 + exp(negate x))
-
-logit :: Floating x => x -> x
--- | The inverse of the logistic.
-logit x = log $ x / (1 - x)
-
--- Lists --
-
-mean :: Fractional x => [x] -> x
--- | Average value of a list of numbers.
-mean = uncurry (/) . foldr (\e (s,c) -> (e+s,c+1)) (0,0)
-
-range :: Double -> Double -> Int -> [Double]
--- | Returns n  numbers from mn to mx.
-range _ _ 0 = []
-range mn mx 1 = [(mn + mx) / 2]
-range mn mx n =
-    [ x * mx + (1 - x) * mn | x <- (/ (fromIntegral n - 1)) . fromIntegral <$> [0 .. n-1] ]
+import Control.DeepSeq hiding (force)
+import Control.Monad.Primitive hiding (internal)
 
-discretizeFunction :: Double -> Double -> Int -> (Double -> Double) -> [(Double,Double)]
--- | Takes range information in the form of a minimum, maximum, and sample count,
--- a function to sample, and returns a list of pairs (x,f(x)) over the specified
--- range.
-discretizeFunction mn mx n f =
-    let rng = range mn mx n
-    in zip rng $ f <$> rng
+import Options.Applicative
 
--- Graveyard --
+import GHC.TypeNats hiding (Mod)
+import GHC.Generics (Generic)
 
-{-
-parMapWH :: (a -> b) -> [a] -> [b]
--- | ParMap using rseq. WH stands for Weak Head (normal form).
-parMapWH = parMap rseq
+import Debug.Trace
+import System.Directory
+import Numeric hiding (log1p,expm1)
+import Numeric.SpecFunctions
 
-parMapDS :: NFData b => (a -> b) -> [a] -> [b]
-parMapDS = parMap rdeepseq
+import Numeric.Natural
+import Data.ByteString (ByteString)
 
-gridSearch
-    :: (a -> Double) -- ^ The error function on the model
-    -> (x -> a) -- ^ A constructor from the parameter being tested to the model
-    -> [x] -- ^ The list of parameter values to test
-    -> (x,[(x,Double)]) -- ^ The best parameter and model, and accompanying statistics
+type NatNumber = Natural
 
--- | A general implementation of a grid search. Returns a triple where the first
--- element is the best parameter, the second is the best model, and the third is a list
--- of the errors calculated at each parameter value.
-gridSearch errorfun constructor xs =
-    let as = constructor <$> xs
-        errs = parMapDS errorfun as
-        x = fst . minimumBy (comparing snd) $ zip xs errs
-    in (x,zip xs errs)
+orderedHeader :: [ByteString] -> Header
+orderedHeader = CSV.header
 
-iterativeOptimization
-    :: (a -> a -> Double) -- ^ Difference measure
-    -> Double -- ^ Differential error threshold
-    -> Int -- ^ Maximum number of iterations (< 1 is interpreted as infinity)
-    -> (a -> a) -- ^ Iterator
-    -> a -- ^ Initial Value
-    -> (a,[(a,Double)]) -- ^ The final value and associated descent
--- | Iterates a value, stopping when a stopping criterion is satisfied, or the maximum
--- number of iterations is reached. The algorithm returns the resulting value as well as
--- the descent preceding it.
---
--- The error accompanying each element of the descent corresponds to the error between
--- that element and the following element. The last error in the descent will therefore
--- correspond to the measured difference between the last element of the descent, and the
--- returned singleton in the pair.
---
--- If the last error in the returned descent is less then the given threshold, then the
--- threshold was reached. Otherwise, n will have been reached as the terminal condition.
-iterativeOptimization difference thrsh n iterator a =
-    let itrs = iterate iterator a
-        itrs' = if n > 0 then take (n+1) itrs else itrs
-        zps = zip itrs' $ tail itrs'
-        zps' = case break ((< thrsh) . snd) . zip zps $ uncurry difference <$> zps of
-                   (zps0,[]) -> zps0
-                   (zps0,rst) -> zps0 ++ [head rst]
-    in (snd . fst . last $ zps', first fst <$> zps')
-    -}
diff --git a/Goal/Core/Circuit.hs b/Goal/Core/Circuit.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Core/Circuit.hs
@@ -0,0 +1,215 @@
+{-# LANGUAGE Arrows,LambdaCase #-}
+-- | A set of functions for working with the 'Arrow' known as a Mealy automata,
+-- here referred to as 'Circuit's. Circuits are essentialy a way of building
+-- composable fold and iterator operations, where some of the values being
+-- processed can be hidden.
+module Goal.Core.Circuit
+    ( -- * Circuits
+    Circuit (Circuit, runCircuit)
+    , accumulateFunction
+    , accumulateCircuit
+    , streamCircuit
+    , iterateCircuit
+    , loopCircuit
+    , loopAccumulator
+    , arrM
+    -- * Chains
+    , Chain
+    , chain
+    , chainCircuit
+    , streamChain
+    , iterateChain
+    , skipChain
+    , skipChain0
+    -- ** Recursive Computations
+    , iterateM
+    , iterateM'
+    ) where
+
+--- Imports ---
+
+
+-- Unqualified --
+
+import Control.Arrow
+
+-- Qualified --
+
+import qualified Control.Category as C
+
+--- Circuits ---
+
+-- | An arrow which takes an input, monadically produces an output, and updates
+-- an (inaccessable) internal state.
+newtype Circuit m a b = Circuit
+    { runCircuit :: a -> m (b, Circuit m a b) }
+
+-- | Takes a function from a value and an accumulator (e.g. just a sum value or
+-- an evolving set of parameters for some model) to a value and an accumulator.
+-- The accumulator is then looped back into the function, returning a Circuit
+-- from a to b, which updates the accumulator every step.
+accumulateFunction :: Monad m => acc -> (a -> acc -> m (b,acc)) -> Circuit m a b
+{-# INLINE accumulateFunction #-}
+accumulateFunction acc f = Circuit $ \a -> do
+    (b,acc') <- f a acc
+    return (b,accumulateFunction acc' f)
+
+-- | accumulateCircuit takes a 'Circuit' and an inital value and loops it.
+accumulateCircuit :: Monad m => acc -> Circuit m (a,acc) (b,acc) -> Circuit m a b
+{-# INLINE accumulateCircuit #-}
+accumulateCircuit acc0 mly0 = accumulateFunction (acc0,mly0) $ \a (acc,Circuit crcf) -> do
+    ((b,acc'),mly') <- crcf (a,acc)
+    return (b,(acc',mly'))
+
+-- | Takes a Circuit and an inital value and loops it, but continues
+-- to return both the output and the accumulated value.
+loopCircuit :: Monad m => acc -> Circuit m (a,acc) (b,acc) -> Circuit m a (b,acc)
+{-# INLINE loopCircuit #-}
+loopCircuit acc0 mly0 = accumulateFunction (acc0,mly0) $ \a (acc,Circuit crcf) -> do
+    ((b,acc'),mly') <- crcf (a,acc)
+    return ((b,acc'),(acc',mly'))
+
+-- | Takes a Circuit which only produces an accumulating value, and loops it.
+loopAccumulator :: Monad m => acc -> Circuit m (a,acc) acc -> Circuit m a acc
+{-# INLINE loopAccumulator #-}
+loopAccumulator acc0 mly0 = accumulateFunction (acc0,mly0) $ \a (acc,Circuit crcf) -> do
+    (acc',mly') <- crcf (a,acc)
+    return (acc',(acc',mly'))
+
+-- | Feeds a list of inputs into a 'Circuit' and returns the (monadic) list of outputs.
+streamCircuit :: Monad m => Circuit m a b -> [a] -> m [b]
+{-# INLINE streamCircuit #-}
+streamCircuit _ [] = return []
+streamCircuit (Circuit mf) (a:as) = do
+    (b,crc') <- mf a
+    (b :) <$> streamCircuit crc' as
+
+-- | Feeds a list of inputs into a Circuit automata and returns the final
+-- monadic output. Throws an error on the empty list.
+iterateCircuit :: Monad m => Circuit m a b -> [a] -> m b
+{-# INLINE iterateCircuit #-}
+iterateCircuit _ [] = error "Empty list fed to iterateCircuit"
+iterateCircuit (Circuit mf) [a] = fst <$> mf a
+iterateCircuit (Circuit mf) (a:as) = do
+    (_,crc') <- mf a
+    iterateCircuit crc' as
+
+-- | Turn a monadic function into a circuit.
+arrM :: Monad m => (a -> m b) -> Circuit m a b
+{-# INLINE arrM #-}
+arrM mf = Circuit $ \a -> do
+    b <- mf a
+    return (b, arrM mf)
+
+
+--- Chains ---
+
+
+-- | A 'Chain' is an iterator built on a 'Circuit'. 'Chain' constructors are
+-- designed to ensure that the first value returned is the initial value of the
+-- iterator (this is not entirely trivial).
+type Chain m x = Circuit m () x
+
+-- | Creates a 'Chain' from an initial state and a transition function. The
+-- first step of the chain returns the initial state, and then continues with
+-- generated states.
+chain
+    :: Monad m
+    => x -- ^ The initial state
+    -> (x -> m x) -- ^ The transition function
+    -> Chain m x -- ^ The resulting 'Chain'
+{-# INLINE chain #-}
+chain x0 mf = accumulateFunction x0 $ \() x -> do
+    x' <- mf x
+    return (x,x')
+
+-- | Creates a 'Chain' from an initial state and a transition circuit. The
+-- first step of the chain returns the initial state, and then continues with
+-- generated states.
+chainCircuit
+    :: Monad m
+    => x -- ^ The initial state
+    -> Circuit m x x -- ^ The transition circuit
+    -> Chain m x -- ^ The resulting 'Chain'
+{-# INLINE chainCircuit #-}
+chainCircuit x0 crc = accumulateCircuit x0 $ proc ((),x) -> do
+    x' <- crc -< x
+    returnA -< (x,x')
+
+-- | Returns the specified number of the given 'Chain's output.
+streamChain :: Monad m => Int -> Chain m x -> m [x]
+{-# INLINE streamChain #-}
+streamChain n chn = streamCircuit chn $ replicate (n+1) ()
+
+-- | Returns the given 'Chain's output at the given index.
+iterateChain :: Monad m => Int -> Chain m x -> m x
+{-# INLINE iterateChain #-}
+iterateChain 0 (Circuit mf) = fst <$> mf ()
+iterateChain k (Circuit mf) = mf () >>= iterateChain (k-1) . snd
+
+-- | Modify the given 'Chain' so that it returns the initial value, and then
+-- skips the specified number of outputs before producing each subsequent output.
+skipChain :: Monad m => Int -> Chain m x -> Chain m x
+{-# INLINE skipChain #-}
+skipChain n (Circuit mf) = Circuit $ \() -> do
+    (x',crc') <- mf ()
+    return (x', skipChain0 n crc')
+
+-- | Modify the given 'Chain' so that it skips the specified number of outputs
+-- before producing each subsequent output (this skips the initial output too).
+skipChain0 :: Monad m => Int -> Chain m x -> Chain m x
+{-# INLINE skipChain0 #-}
+skipChain0 n crc = Circuit $ \() -> do
+    (Circuit mf) <- iterateM' n iterator crc
+    (x',crc') <- mf ()
+    return (x', skipChain0 n crc')
+        where iterator (Circuit mf') = snd <$> mf' ()
+
+
+-- | Iterate a monadic action the given number of times, returning the complete
+-- sequence of values.
+iterateM :: Monad m => Int -> (x -> m x) -> x -> m [x]
+{-# INLINE iterateM #-}
+iterateM n mf x0 = streamChain n $ chain x0 mf
+
+-- | Iterate a monadic action the given number of times, returning the final value.
+iterateM' :: Monad m => Int -> (x -> m x) -> x -> m x
+{-# INLINE iterateM' #-}
+iterateM' n mf x0 = iterateChain n $ chain x0 mf
+
+
+
+--- Instances ---
+
+
+instance Monad m => C.Category (Circuit m) where
+    --id :: Circuit a a
+    {-# INLINE id #-}
+    id = Circuit $ \a -> return (a,C.id)
+    --(.) :: Circuit b c -> Circuit a b -> Circuit a c
+    {-# INLINE (.) #-}
+    (.) = dot
+        where dot (Circuit crc1) (Circuit crc2) = Circuit $ \a -> do
+                  (b, crcA2') <- crc2 a
+                  (c, crcA1') <- crc1 b
+                  return (c, crcA1' `dot` crcA2')
+
+instance Monad m => Arrow (Circuit m) where
+    --arr :: (a -> b) -> Circuit a b
+    {-# INLINE arr #-}
+    arr f = Circuit $ \a -> return (f a, arr f)
+    --first :: Circuit a b -> Circuit (a,c) (b,c)
+    {-# INLINE first #-}
+    first (Circuit crcf) = Circuit $ \(a,c) -> do
+        (b, crcA') <- crcf a
+        return ((b,c), first crcA')
+
+instance Monad m => ArrowChoice (Circuit m) where
+    --left :: Circuit a b -> Circuit (Either a c) (Either b c)
+    {-# INLINE left #-}
+    left crcA@(Circuit crcf) = Circuit $
+        \case
+          Left a -> do
+              (b,crcA') <- crcf a
+              return (Left b,left crcA')
+          Right c -> return (Right c,left crcA)
diff --git a/Goal/Core/Plot.hs b/Goal/Core/Plot.hs
deleted file mode 100644
--- a/Goal/Core/Plot.hs
+++ /dev/null
@@ -1,264 +0,0 @@
-module Goal.Core.Plot
-    ( -- * Module Exports
-      module Graphics.Rendering.Chart
-    , module Data.Colour
-    , module Data.Colour.Names
-    , module Data.Colour.SRGB.Linear
-    , module Graphics.Rendering.Chart.Backend.Cairo
-    , module Graphics.Rendering.Chart.Grid
-    , module Graphics.Rendering.Chart.Gtk
-    , module Graphics.Rendering.Chart.State
-    , module Goal.Core.Plot.Contour
-    -- * Plots
-    -- ** PixMap
-    , pixMapPlot
-    -- ** Histograms
-    , histogramPlot
-    , histogramPlot0
-    , logHistogramPlot
-    , logHistogramPlot0
-    -- * Layouts
-    -- ** PixMap
-    , pixMapLayout
-    -- ** Histogram
-    , histogramLayout
-    , logHistogramLayout
-    , histogramLayoutLR
-    -- * Util
-    , rgbaGradient
-    -- * Rendering
-    , renderableToAspectWindow
-    ) where
-
-
---- Imports ---
-
-import Data.List hiding (sum)
-
-import Control.Monad
-import Control.Lens.Getter
-import Control.Lens.Setter hiding (Identity)
-
-import Data.Default.Class
-import Numeric
-
-
--- Re-exports --
-
-import Graphics.Rendering.Chart hiding (x0,y0,Point)
-import Data.Colour
-import Data.Colour.Names
-import Data.Colour.SRGB.Linear
-import Graphics.Rendering.Chart.Backend.Cairo
-import Graphics.Rendering.Chart.State
-import Graphics.Rendering.Chart.Grid
-import Graphics.Rendering.Chart.Gtk
-
--- Scientific --
-
-import Goal.Core.Plot.Contour
-
--- Qualified --
-
-import qualified Graphics.UI.Gtk as G
-
-
--- Unqualified --
-
-import Graphics.Rendering.Cairo (liftIO)
-
-
---- Util ---
-
-rgbaGradient :: (Double, Double, Double, Double) -> (Double, Double, Double, Double) -> Int
-    -> [AlphaColour Double]
--- | Returns an ordered list of colours useful for plotting.
-rgbaGradient (rmn,gmn,bmn,amn) (rmx,gmx,bmx,amx) n =
-    zipWith (flip withOpacity) [amn,amn + astp .. amx]
-    $ zipWith3 rgb [rmn,rmn + rstp .. rmx] [gmn,gmn + gstp .. gmx] [bmn,bmn + bstp .. bmx]
-    where rstp = (rmx - rmn) / fromIntegral n
-          gstp = (gmx - gmn) / fromIntegral n
-          bstp = (bmx - bmn) / fromIntegral n
-          astp = (amx - amn) / fromIntegral n
-
-
---- Plots ---
-
-
--- PixMap --
-
-pixMapPlot :: (Double,Double) -> [[AlphaColour Double]] -> Plot Double Double
--- | Returns a pixmap representation of a matrix style set of doubles. Based on the
--- defaults, the list of colours are assumed to be in (y,x) coordinates, where the
--- origin is at the lower left of the image. If matrix style coordinates are desired,
--- The containing layout should be given a reversed y axis, so that the origin is at the
--- top left of the image.
-
---  The pair of doubles indicates where corner of the given image should be located.
-pixMapPlot (x0,y0) pss = foldr1 joinPlot
-    $ concat [ [ boxPlot r c . solidFillStyle $ p | (c,p) <- zip [x0..] ps ] | (r,ps) <- zip [y0..] pss ]
-    where boxPlot y x stl = toPlot
-              $ plot_fillbetween_style .~ stl
-              $ plot_fillbetween_values .~ [(x-0.5,(y-0.5,y+0.5)),(x+0.5,(y-0.5,y+0.5))]
-              $ def
-
-pixMapLayout :: Int -> Int -> Layout Double Double -> Layout Double Double
--- | A nice base layout for a pixMap, with a box around the pixmap with one 'pixel'
--- padding, and a reversed y axis for matrix style coordinates.
-pixMapLayout rws0 cls0 lyt =
-    layout_top_axis_visibility .~ AxisVisibility True True False
-    $ layout_left_axis_visibility .~ AxisVisibility True True False
-    $ layout_right_axis_visibility .~ AxisVisibility True True False
-    $ layout_bottom_axis_visibility .~ AxisVisibility True True False
-    $ layout_y_axis . laxis_reverse .~ True
-    $ layout_y_axis . laxis_generate .~
-        const (makeAxis (const "") ([-1.5,rws+0.5],[-1.5,rws+0.5],[-1.5,rws+0.5]))
-    $ layout_x_axis . laxis_generate .~
-        const (makeAxis (const "") ([-1.5,cls+0.5],[-1.5,cls+0.5],[-1.5,cls+0.5]))
-    $ lyt
-    where cls = fromIntegral cls0
-          rws = fromIntegral rws0
-
--- Histogram --
-
-histogramPlot0 :: (Num a,BarsPlotValue a) =>Int -> [[Double]] -> PlotBars Double a -> PlotBars Double a
--- | Generates a histogram plot where the min and max bin value is taken from the data set.
-histogramPlot0 n xss plt =
-    let mx = maximum $ maximum <$> xss
-        mn = minimum $ minimum <$> xss
-    in histogramPlot n mn mx xss plt
-
-logHistogramPlot0 :: Int -> [[Double]] -> PlotBars Double Double -> PlotBars Double Double
--- | Generates a histogram plot where the min and max bin value is taken from the data set.
-logHistogramPlot0 n xss =
-    let mx = maximum $ maximum <$> xss
-        mn = minimum $ minimum <$> xss
-    in logHistogramPlot n mn mx xss
-
-histogramPlot
-    :: (Num a,BarsPlotValue a)
-    => Int -- ^ Number of bins
-    -> Double -- ^ Min range
-    -> Double -- ^ Max range
-    -> [[Double]] -- ^ Data set
-    -> PlotBars Double a -- ^ Plot
-    -> PlotBars Double a -- ^ New Plot
--- | Creates a histogram out of a data set. The data set is a list of list of values, where
--- each sublist is a collection of data along an axis. Under and overflow is put into the
--- first and last bin, respectively. The bars are centered at the mid point between each
--- pair of bins.
-histogramPlot n mn mx xss plt =
-    let bns = range mn mx (n+1)
-        vls = transpose $ toHistogram (tail $ take n bns) . sort <$> xss
-        stp = (head (tail bns) - head bns) / 2
-        bns' = (+ stp) <$> take n bns
-     in plot_bars_alignment .~ BarsCentered $ plot_bars_values .~ zip bns' vls $ plt
-    where toHistogram _ [] = repeat 0
-          toHistogram [] xs = [genericLength xs]
-          toHistogram (bn:bns') xs =
-              let (hds,xs') = span (< bn) xs
-              in genericLength hds : toHistogram bns' xs'
-
-range _ _ 0 = []
-range mn mx 1 = [(mn + mx) / 2]
-range mn mx n =
-    [ x * mx + (1 - x) * mn | x <- (/ (fromIntegral n - 1)) . fromIntegral <$> [0 .. n-1] ]
-
-logHistogramPlot
-    :: Int -- ^ Number of bins
-    -> Double -- ^ Min range
-    -> Double -- ^ Max range
-    -> [[Double]] -- ^ Data set
-    -> PlotBars Double Double -- ^ Plot
-    -> PlotBars Double Double -- ^ New Plot
-logHistogramPlot n mn mx xss plt =
-    let bplt = histogramPlot n mn mx xss plt
-        vls = modVals <$> bplt ^. plot_bars_values
-        --cl = fromIntegral . ceiling . maximum . concat $ snd <$> vls
-    in plot_bars_values .~ vls $ bplt
-    where lbs = 10
-          modVals (x,ys) = (x, logBase lbs . (+1) <$> ys)
-
-histogramLayout :: BarsPlotValue a => PlotBars Double a -> Layout Double a -> Layout Double a
--- | The base layout for a histogram.
-histogramLayout pbrs lyt =
-    let bns = fst <$> pbrs ^. plot_bars_values
-        stp = (head (tail bns) - head bns) / 2
-        bns' = (head bns - stp) : map (+stp) bns
-        rng = abs $ maximum bns'
-        labelFun x
-            | rng >= 1000 = showEFloat (Just 2) x ""
-            | rng <= 0.01 = showEFloat (Just 2) x ""
-            | otherwise = reverse . dropWhile (== '.') . dropWhile (== '0') . reverse $ showFFloat (Just 2) x ""
-    in layout_x_axis . laxis_generate .~ const (makeAxis labelFun (bns',[],bns'))
-       $ layout_plots %~ (plotBars pbrs:)
-       $ lyt
-
-logHistogramLayout :: PlotBars Double Double -> Layout Double Double -> Layout Double Double
--- | Base layout for a log-histogram.
-logHistogramLayout pbrs lyt =
-    let vls = pbrs ^. plot_bars_values
-        bns = fst <$> vls
-        stp = (head (tail bns) - head bns) / 2
-        bns' = (head bns - stp) : map (+stp) bns
-        cl = fromIntegral . ceiling . maximum . concat $ snd <$> vls
-        rng = abs $ maximum bns'
-        xLabelFun x
-            | rng >= 1000 = showEFloat (Just 2) x ""
-            | rng <= 0.01 = showEFloat (Just 2) x ""
-            | otherwise = reverse . dropWhile (== '.') . dropWhile (== '0') . reverse $ showFFloat (Just 2) x ""
-    in layout_plots %~ (plotBars pbrs:)
-        $ layout_y_axis . laxis_generate .~ const (makeAxis yLabelFun ([0..cl],[],[0..cl]))
-        $ layout_x_axis . laxis_generate .~ const (makeAxis xLabelFun (bns',[],bns'))
-        $ lyt
-    where lbs = 10
-          yLabelFun 0 = show 0
-          yLabelFun x = show (round lbs) ++ "e" ++ show (round x) ++ "-1"
-
-histogramLayoutLR :: (BarsPlotValue a,PlotValue b) => PlotBars Double a -> LayoutLR Double a b -> LayoutLR Double a b
--- | The base layout for a histogramLR.
-histogramLayoutLR pbrs lyt =
-    let bns = fst <$> pbrs ^. plot_bars_values
-        stp = (head (tail bns) - head bns) / 2
-        bns' = (head bns - stp) : map (+stp) bns
-        rng = abs $ maximum bns'
-        labelFun x
-            | rng >= 1000 = showEFloat (Just 2) x ""
-            | rng <= 0.01 = showEFloat (Just 2) x ""
-            | otherwise = reverse . dropWhile (== '.') . dropWhile (== '0') . reverse $ showFFloat (Just 2) x ""
-    in layoutlr_x_axis . laxis_generate .~ const (makeAxis labelFun (bns',[],bns'))
-       $ layoutlr_plots %~ (Left (plotBars pbrs):)
-       $ lyt
-
---- IO ---
-
-
-renderableToAspectWindow
-    :: Bool -- ^ Display Full Screen
-    -> Int -- ^ Image width
-    -> Int -- ^ Image height
-    -> Renderable a -- ^ The Renderable
-    -> IO () -- ^ Renders the renderable to the screen
-
--- | Displays a renderable in a GTK aspect window.
-renderableToAspectWindow fs wdth hght rnbl = do
-
-    G.initGUI
-
-    win <- G.windowNew
-    afrm <- G.aspectFrameNew 0.5 0.5 . Just $ realToFrac (fromIntegral wdth / fromIntegral hght)
-    da <- G.drawingAreaNew
-
-    G.set afrm [ G.containerChild G.:= da ]
-    G.set win [ G.containerChild G.:= afrm ]
-    when fs $ G.windowFullscreen win
-    G.onDestroy win G.mainQuit
-
-    (da `G.on` G.exposeEvent) . liftIO $ do
-
-        updateCanvas rnbl da
-        return True
-
-    G.widgetShowAll win
-    G.mainGUI
-
diff --git a/Goal/Core/Plot/Contour.hs b/Goal/Core/Plot/Contour.hs
deleted file mode 100644
--- a/Goal/Core/Plot/Contour.hs
+++ /dev/null
@@ -1,163 +0,0 @@
-module Goal.Core.Plot.Contour (contours) where
-
-
---- Imports ---
-
-
--- General --
-
-import Control.Monad
-import Data.List
-
--- Qualified --
-
-import qualified Data.Map as M
-
-
--- Unqualified --
-
-import Data.Tuple (swap)
-
-
---- Contour Plot ---
-
-
-contours
-    :: (Double,Double,Int) -- ^ The range along the x axis.
-    -> (Double,Double,Int) -- ^ The range along the y axis.
-    -> Int -- ^ The number of isolevels.
-    -> (Double -> Double -> Double) -- ^ The function to contour.
-    -> [(Double,[[(Double,Double)]])] -- ^ (isolevel, (x,y))
-
--- | Given various parameters and a function to analyze, contourPlot returns a
--- list of PlotLines each of which is an isoline.
-contours (xmn,xmx,nx) (ymn,ymx,ny) nlvls f =
-    let lsts = functionToLists f (xmn,ymn) (xmx,ymx) stps
-        stps = ((xmx - xmn) / fromIntegral nx,(ymx - ymn) / fromIntegral ny)
-        (mn,mx) = (minimum $ minimum <$> lsts,maximum $ maximum <$> lsts)
-        isostp = (mx - mn) / fromIntegral nlvls
-        isolvls = take nlvls [mn + (isostp / 2),mn + (isostp * 1.5)..]
-     in contour lsts (xmx,ymx) stps <$> isolvls
-
-
---- Internal ---
-
-
-contour :: [[Double]] -> (Double,Double) -> (Double,Double) -> Double -> (Double,[[(Double,Double)]])
-contour lsts mxs stps isolvl = (isolvl,linksToLines mxs stps <$> listsToLinks lsts isolvl)
-
-
---- Contour Plot Internal ---
-
-
--- Types --
-
-type Link = (Int,Int)
-type ContourPair = (Link,Link)
-data ContourBox =
-    {- Styled like the indicies of a 3x3 Matrix. The lines are drawn from
-       left to right, and if need be, top to bottom -}
-    Empty
-    | Line ContourPair
-    | DLine ContourPair ContourPair
-    deriving (Show,Eq)
-
-data PairMap = PM (M.Map Link Link) (M.Map Link Link) deriving Show
-
--- Type Functions --
-
-contourBoxesToPairMap :: [ContourBox] -> PairMap
-contourBoxesToPairMap clns =
-    let cprs = concatMap toPairs clns
-    in PM (M.fromList cprs) (M.fromList $ map swap cprs)
-    where toPairs (Line prpr) = [prpr]
-          toPairs (DLine lprpr rprpr) = [lprpr,rprpr]
-
-deletePair :: ContourPair -> PairMap -> PairMap
-{- Deletes a left oriented line segmented -}
-deletePair (lpr,rpr) (PM lmp rmp) =
-    let lmp' = if M.lookup lpr lmp == Just rpr then M.delete lpr lmp else lmp
-        rmp' = if M.lookup rpr rmp == Just lpr then M.delete rpr rmp else rmp
-    in PM lmp' rmp'
-
-popLink :: Link -> PairMap -> (Maybe Link,PairMap)
-popLink lnk pmp@(PM lmp rmp) =
-    let lft = M.lookup lnk lmp
-        rgt = M.lookup lnk rmp
-    in case (lft,rgt) of
-           (Just rlnk,_) -> (Just rlnk,deletePair (lnk,rlnk) pmp)
-           (_,Just llnk) -> (Just llnk,deletePair (llnk,lnk) pmp)
-           _ -> (Nothing,pmp)
-
--- Core Algorithm --
-
-functionToLists :: (Double -> Double -> Double) -> (Double,Double) -> (Double,Double)
-    -> (Double,Double) -> [[Double]]
-functionToLists f (xmn,ymn) (xmx,ymx) (xstp,ystp) =
-    map (uncurry f) <$> [ [ (x,y) | x <- [xmx,(xmx - xstp)..xmn] ] | y <- [ymx,(ymx - ystp)..ymn] ]
-
-situationTable :: Double -> (((Bool,Double),(Bool,Double)),((Bool,Double),(Bool,Double))) -> ContourBox
-situationTable _ (((True,_),(True,_)),((True,_),(True,_))) = Empty
-situationTable _ (((True,_),(True,_)),((False,_),(True,_))) = Line ((1,0),(2,1))
-situationTable _ (((True,_),(True,_)),((True,_),(False,_))) = Line ((2,1),(1,2))
-situationTable _ (((True,_),(True,_)),((False,_),(False,_))) = Line ((1,0),(1,2))
-situationTable _ (((True,_),(False,_)),((True,_),(True,_))) = Line ((0,1),(1,2))
-situationTable isolvl (((True,ul),(False,ur)),((False,ll),(True,lr)))
-    | sum [ul,ur,ll,lr] / 4 < isolvl = DLine ((1,0),(0,1)) ((2,1),(1,2))
-    | otherwise = DLine ((1,0),(2,1)) ((0,1),(1,2))
-situationTable _ (((True,_),(False,_)),((True,_),(False,_))) = Line ((0,1),(2,1))
-situationTable _ (((True,_),(False,_)),((False,_),(False,_))) = Line ((1,0),(0,1))
-situationTable _ (((False,_),(True,_)),((True,_),(True,_))) = Line ((1,0),(0,1))
-situationTable _ (((False,_),(True,_)),((False,_),(True,_))) = Line ((0,1),(2,1))
-situationTable isolvl (((False,ul),(True,ur)),((True,ll),(False,lr)))
-    | sum [ul,ur,ll,lr] / 4 < isolvl = DLine ((1,0),(2,1)) ((0,1),(1,2))
-    | otherwise = DLine ((1,0),(0,1)) ((2,1),(1,2))
-situationTable _ (((False,_),(True,_)),((False,_),(False,_))) = Line ((0,1),(1,2))
-situationTable _ (((False,_),(False,_)),((True,_),(True,_))) = Line ((1,0),(1,2))
-situationTable _ (((False,_),(False,_)),((False,_),(True,_))) = Line ((2,1),(1,2))
-situationTable _ (((False,_),(False,_)),((True,_),(False,_))) = Line ((1,0),(2,1))
-situationTable _ (((False,_),(False,_)),((False,_),(False,_))) = Empty
-
-listsToContourBoxes :: [[Double]] -> Double -> [ContourBox]
-listsToContourBoxes lsts isolvl = do
-    (stnrw,r) <- zip stnlsts [0..]
-    (stn,c) <- zip stnrw [0..]
-    let bx = situationTable isolvl stn
-    guard (bx /= Empty)
-    return $ repositionBox (r,c) bx
-    where stnlsts = rowZipper $ elementPairs . map threshold <$> lsts
-          threshold x = (x >= isolvl,x)
-          elementPairs lst = zip lst $ tail lst
-          rowZipper rws = uncurry zip <$> elementPairs rws
-          repositionContourPair (r,c) ((lr,lc),(rr,rc)) =
-              ((lr + 2 * r, lc + 2 * c),(rr + 2 * r, rc + 2 * c))
-          repositionBox rc (Line pr) =
-              Line $ repositionContourPair rc pr
-          repositionBox rc (DLine pr1 pr2) =
-              DLine (repositionContourPair rc pr1) (repositionContourPair rc pr2)
-
-traceLinks :: ContourPair -> PairMap -> ([Link],PairMap)
-traceLinks (llnk,rlnk) pmp =
-    let (llnks,pmp') = tracer [] (Just llnk,pmp)
-        (rlnks,pmp'') = tracer [] (Just rlnk,pmp')
-    in (llnks ++ reverse rlnks,pmp'')
-    where tracer lnks (Just lnk,pmp') =
-              tracer (lnk:lnks) $ popLink lnk pmp'
-          tracer lnks (Nothing,pmp') = (lnks,pmp')
-
-popLinks :: PairMap -> Maybe ([Link],PairMap)
-popLinks pmp@(PM lmp _)
-    | M.null lmp = Nothing
-    | otherwise =
-        let cpr = M.findMin lmp
-        in Just (traceLinks cpr $ deletePair cpr pmp)
-
-listsToLinks :: [[Double]] -> Double -> [[(Int,Int)]]
-listsToLinks lsts isolvl =
-    unfoldr popLinks .  contourBoxesToPairMap $ listsToContourBoxes lsts isolvl
-
-linksToLines :: (Double,Double) -> (Double,Double) -> [(Int,Int)] -> [(Double,Double)]
-linksToLines (xmx,ymx) (xstp,ystp) lnks =
-   (\(r,c) -> (xmx - fromIntegral c * xstp / 2,ymx - fromIntegral r * ystp / 2)) <$> lnks
-
-
diff --git a/Goal/Core/Project.hs b/Goal/Core/Project.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Core/Project.hs
@@ -0,0 +1,185 @@
+{-# LANGUAGE OverloadedStrings #-}
+-- | This module provides functions for incorporating Goal into a
+-- data-processing project. In particular, this module provides tools for
+-- managing CSV files, and connecting them with gnuplot scripts for plotting.
+-- CSV management is powered by @cassava@.
+module Goal.Core.Project
+    (
+    -- * CSV
+      goalImport
+    , goalImportNamed
+    , goalExport
+    , goalExportLines
+    , goalExportNamed
+    , goalExportNamedLines
+    -- ** CSV Instances
+    , goalCSVParser
+    , goalCSVNamer
+    , goalCSVOrder
+    -- * Util
+    , runGnuplot
+    , runGnuplotWithVariables
+    ) where
+
+
+--- Imports ---
+
+
+-- Unqualified --
+
+import System.Process
+import System.Directory
+import Data.Csv
+import Data.Char
+import GHC.Generics
+
+-- Qualified --
+
+import qualified Data.Vector as V
+import qualified Data.ByteString.Lazy as BS
+import qualified Data.ByteString as BSI
+
+
+--- Experiments ---
+
+
+--- Import/Export ---
+
+
+-- | Runs @gnuplot@ on the given @.gpi@, passing it a @load_path@ variable to
+-- help it find Goal-generated csvs.
+runGnuplot
+    :: FilePath -- ^ Gnuplot loadpath
+    -> String -- ^ Gnuplot script
+    -> IO ()
+runGnuplot ldpth gpipth = do
+    let cmd = concat [ "gnuplot ", " -e \"load_path='", ldpth, "'\" ",gpipth,".gpi" ]
+    putStrLn $ "Running Command: " ++ cmd
+    callCommand cmd
+
+-- | Runs @gnuplot@ on the given @.gpi@, passing it a @load_path@ variable to
+-- help it find Goal-generated csvs, and a list of variables.
+runGnuplotWithVariables
+    :: FilePath -- ^ Gnuplot loadpath
+    -> String -- ^ Gnuplot script
+    -> [(String,String)] -- ^ Arguments
+    -> IO ()
+runGnuplotWithVariables ldpth gpipth args = do
+    let cmd = concat $ [ "gnuplot ", " -e \"load_path='", ldpth, "'" ]
+            ++ (mapArgs <$> args) ++ [ "\" ",gpipth,".gpi" ]
+    putStrLn $ "Running Command: " ++ cmd
+    callCommand cmd
+        where mapArgs (nm,val) = concat ["; ",nm,"='",val,"'"]
+
+
+-- | Load the given CSV file. The @.csv@ extension is automatically added.
+goalImport
+    :: FromRecord r
+    => FilePath
+    -> IO (Either String [r]) -- ^ CSVs
+goalImport flpth = do
+    bstrm <- decode NoHeader <$> BS.readFile (flpth ++ ".csv")
+    case bstrm of
+      Right as -> return . Right $ V.toList as
+      Left str -> return $ Left str
+
+-- | Load the given CSV file with headers. The @.csv@ extension is automatically added.
+goalImportNamed
+    :: FromNamedRecord r
+    => FilePath
+    -> IO (Either String [r]) -- ^ CSVs
+goalImportNamed flpth = do
+    bstrm <- decodeByName <$> BS.readFile (flpth ++ ".csv")
+    case bstrm of
+      Right as -> return . Right . V.toList $ snd as
+      Left str -> return $ Left str
+
+filePather :: FilePath -> FilePath -> IO FilePath
+filePather ldpth flnm = do
+    createDirectoryIfMissing True ldpth
+    return $ concat [ldpth,"/",flnm,".csv"]
+
+-- | Export the given CSVs to a file in the given directory. The @.csv@
+-- extension is automatically added to the file name.
+goalExport
+    :: ToRecord r
+    => FilePath -- load_path
+    -> String -- File Name
+    -> [r] -- ^ CSVs
+    -> IO ()
+goalExport ldpth flnm csvs = do
+    flpth <- filePather ldpth flnm
+    BS.writeFile flpth $ encode csvs
+
+-- | Export the given list of CSVs to a file in the given directory, seperating
+-- each set of CSVs by a single line. This causes gnuplot to the read CSV as a
+-- collection of line segments. The @.csv@ extension is automatically added to
+-- the file name.
+goalExportLines
+    :: ToRecord r
+    => FilePath
+    -> FilePath
+    -> [[r]] -- ^ CSVss
+    -> IO ()
+goalExportLines ldpth flnm csvss = do
+    flpth <- filePather ldpth flnm
+    BS.writeFile flpth . BS.tail . BS.tail . BS.concat $ BS.append "\r\n" . encode <$> csvss
+
+-- | Export the named CSVs to a file in the given directory, adding a header to
+-- the @.csv@ file.
+goalExportNamed
+    :: (ToNamedRecord r, DefaultOrdered r)
+    => FilePath
+    -> FilePath
+    -> [r] -- ^ CSVs
+    -> IO ()
+goalExportNamed ldpth flnm csvs = do
+    flpth <- filePather ldpth flnm
+    BS.writeFile flpth $ encodeDefaultOrderedByName csvs
+
+-- | Export the given list of named CSVs to a file, breaking it into a set of
+-- line segments (with headers).
+goalExportNamedLines
+    :: (ToNamedRecord r, DefaultOrdered r)
+    => FilePath
+    -> FilePath
+    -> [[r]] -- ^ CSVss
+    -> IO ()
+goalExportNamedLines ldpth flnm csvss = do
+    flpth <- filePather ldpth flnm
+    BS.writeFile flpth . BS.concat $ BS.append "\r\n" . encodeDefaultOrderedByName <$> csvss
+
+
+--- Util ---
+
+
+deCamelCaseLoop :: String -> String
+deCamelCaseLoop "" = ""
+deCamelCaseLoop (c:wrds) =
+    let (wrd,wrds') = span isLower wrds
+     in (c:wrd) ++ ' ' : deCamelCaseLoop wrds'
+
+deCamelCase :: String -> String
+deCamelCase (c:wrds) = init $ deCamelCaseLoop (toUpper c : wrds)
+deCamelCase "" = error "How is deCamelCase being run on an empty string?"
+
+deCamelCaseCSV :: Options
+deCamelCaseCSV = defaultOptions { fieldLabelModifier = deCamelCase }
+
+-- | A generic @.csv@ parser which reorganizes a header name in camel case into
+-- "human readable" text. Useful for instantiating 'FromNamedRecord'.
+goalCSVParser :: (Generic a, GFromNamedRecord (Rep a)) => NamedRecord -> Parser a
+goalCSVParser = genericParseNamedRecord deCamelCaseCSV
+
+-- | A generic @.csv@ namer which reorganizes a header name in camel case into
+-- "human readable" text. Useful for instantiating 'ToNamedRecord'.
+goalCSVNamer
+    :: (Generic a, GToRecord (Rep a) (BSI.ByteString, BSI.ByteString)) => a -> NamedRecord
+goalCSVNamer = genericToNamedRecord deCamelCaseCSV
+
+-- | A generic @.csv@ order which reorganizes a header name in camel case into
+-- "human readable" text. Useful for instantiating 'DefaultOrdered'.
+goalCSVOrder :: (Generic a, GToNamedRecordHeader (Rep a)) => a -> Header
+goalCSVOrder = genericHeaderOrder deCamelCaseCSV
+
+
diff --git a/Goal/Core/Util.hs b/Goal/Core/Util.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Core/Util.hs
@@ -0,0 +1,265 @@
+-- | A collection of generic numerical and list manipulation functions.
+module Goal.Core.Util
+    ( -- * List Manipulation
+      takeEvery
+    , breakEvery
+    , kFold
+    , kFold'
+    -- * Numeric
+    , roundSD
+    , toPi
+    , circularDistance
+    , integrate
+    , logistic
+    , logit
+    , square
+    , triangularNumber
+    -- ** List Numerics
+    , average
+    , weightedAverage
+    , circularAverage
+    , weightedCircularAverage
+    , range
+    , discretizeFunction
+    , logSumExp
+    , logIntegralExp
+    -- * Tracing
+    , traceGiven
+    -- * TypeNats
+    , finiteInt
+    , natValInt
+    , Triangular
+    -- ** Type Rationals
+    , Rat
+    , type (/)
+    , ratVal
+    ) where
+
+
+--- Imports ---
+
+
+-- Unqualified --
+
+import Numeric
+import Data.Ratio
+import Data.Proxy
+import Debug.Trace
+import Data.Finite
+import GHC.TypeNats
+
+-- Qualified --
+
+import qualified Numeric.GSL.Integration as I
+import qualified Data.List as L
+
+--- General Functions ---
+
+
+-- | Takes every nth element, starting with the head of the list.
+takeEvery :: Int -> [x] -> [x]
+{-# INLINE takeEvery #-}
+takeEvery m = map snd . filter (\(x,_) -> mod x m == 0) . zip [0..]
+
+-- | Break the list up into lists of length n.
+breakEvery :: Int -> [x] -> [[x]]
+{-# INLINE breakEvery #-}
+breakEvery _ [] = []
+breakEvery n xs = take n xs : breakEvery n (drop n xs)
+
+-- | Runs traceShow on the given element.
+traceGiven :: Show a => a -> a
+{-# INLINE traceGiven #-}
+traceGiven a = traceShow a a
+
+
+--- Numeric ---
+
+-- | Numerically integrates a 1-d function over an interval.
+integrate
+    :: Double -- ^ Error Tolerance
+    -> (Double -> Double) -- ^ Function
+    -> Double -- ^ Interval beginning
+    -> Double -- ^ Interval end
+    -> (Double,Double) -- ^ Integral
+{-# INLINE integrate #-}
+integrate errbnd = I.integrateQAGS errbnd 10000
+
+-- | Rounds the number to the specified significant digit.
+roundSD :: RealFloat x => Int -> x -> x
+{-# INLINE roundSD #-}
+roundSD n x =
+    let n' :: Int
+        n' = round $ 10^n * x
+     in fromIntegral n'/10^n
+
+-- | Value of a point on a circle, minus rotations.
+toPi :: RealFloat x => x -> x
+{-# INLINE toPi #-}
+toPi x =
+    let xpi = x / (2*pi)
+        f = xpi - fromIntegral (floor xpi :: Int)
+     in 2 * pi * f
+
+-- | Distance between two points on a circle, removing rotations.
+circularDistance :: RealFloat x => x -> x -> x
+{-# INLINE circularDistance #-}
+circularDistance x y =
+    let x' = toPi x
+        y' = toPi y
+     in min (toPi $ x' - y') (toPi $ y' - x')
+
+-- | A standard sigmoid function.
+logistic :: Floating x => x -> x
+{-# INLINE logistic #-}
+logistic x = 1 / (1 + exp (negate x))
+
+-- | The inverse of the logistic.
+logit :: Floating x => x -> x
+{-# INLINE logit #-}
+logit x = log $ x / (1 - x)
+
+-- | The square of a number (for avoiding endless default values).
+square :: Floating x => x -> x
+{-# INLINE square #-}
+square x = x^(2::Int)
+
+-- | Triangular number.
+triangularNumber :: Int -> Int
+{-# INLINE triangularNumber #-}
+triangularNumber n = flip div 2 $ n * (n+1)
+
+
+-- Lists --
+
+-- | Average value of a 'Traversable' of 'Fractional's.
+average :: (Foldable f, Fractional x) => f x -> x
+{-# INLINE average #-}
+average = uncurry (/) . L.foldl' (\(s,c) e -> (e+s,c+1)) (0,0)
+
+-- | Weighted Average given a 'Traversable' of (weight,value) pairs.
+weightedAverage :: (Foldable f, Fractional x) => f (x,x) -> x
+{-# INLINE weightedAverage #-}
+weightedAverage = uncurry (/) . L.foldl' (\(sm,nrm) (w,x) -> (sm + w*x,nrm + w)) (0,0)
+
+-- | Circular average value of a 'Traversable' of radians.
+circularAverage :: (Traversable f, RealFloat x) => f x -> x
+{-# INLINE circularAverage #-}
+circularAverage rds =
+    let snmu = average $ sin <$> rds
+        csmu = average $ cos <$> rds
+     in atan2 snmu csmu
+
+-- | Returns k (training,validation) pairs. k should be greater than or equal to 2.
+kFold :: Int -> [x] -> [([x],[x])]
+{-# INLINE kFold #-}
+kFold k xs =
+    let nvls = ceiling . (/(fromIntegral k :: Double)) . fromIntegral $ length xs
+     in L.unfoldr unfoldFun ([], breakEvery nvls xs)
+    where unfoldFun (_,[]) = Nothing
+          unfoldFun (hds,tl:tls) = Just ((concat $ hds ++ tls,tl),(tl:hds,tls))
+
+-- | Returns k (training,test,validation) pairs for early stopping algorithms. k
+-- should be greater than or equal to 3.
+kFold' :: Int -> [x] -> [([x],[x],[x])]
+{-# INLINE kFold' #-}
+kFold' k xs =
+    let nvls = ceiling . (/(fromIntegral k :: Double)) . fromIntegral $ length xs
+        brks = breakEvery nvls xs
+     in L.unfoldr unfoldFun ([], brks)
+    where unfoldFun (hds,tl:tl':tls) = Just ((concat $ hds ++ tls,tl,tl'),(tl:hds,tl':tls))
+          unfoldFun (hds,tl:tls) =
+              let (tl0:hds') = reverse hds
+               in Just ((concat $ reverse hds' ++ tls,tl,tl0),(tl:hds,tls))
+          unfoldFun (_,[]) = Nothing
+
+
+-- | Weighted Circular average value of a 'Traversable' of radians.
+weightedCircularAverage :: (Traversable f, RealFloat x) => f (x,x) -> x
+{-# INLINE weightedCircularAverage #-}
+weightedCircularAverage wxs =
+    let snmu = weightedAverage $ sinPair <$> wxs
+        csmu = weightedAverage $ cosPair <$> wxs
+     in atan2 snmu csmu
+    where sinPair (w,rd) = (w,sin rd)
+          cosPair (w,rd) = (w,cos rd)
+
+-- | Returns n numbers which uniformly partitions the interval [mn,mx].
+range
+    :: RealFloat x => x -> x -> Int -> [x]
+{-# INLINE range #-}
+range _ _ 0 = []
+range mn mx 1 = [(mn + mx) / 2]
+range mn mx n =
+    [ x * mx + (1 - x) * mn | x <- (/ (fromIntegral n - 1)) . fromIntegral <$> [0 .. n-1] ]
+
+-- | Takes range information in the form of a minimum, maximum, and sample count,
+-- a function to sample, and returns a list of pairs (x,f(x)) over the specified
+-- range.
+discretizeFunction :: Double -> Double -> Int -> (Double -> Double) -> [(Double,Double)]
+{-# INLINE discretizeFunction #-}
+discretizeFunction mn mx n f =
+    let rng = range mn mx n
+    in zip rng $ f <$> rng
+
+-- | Given a set of values, computes the "soft maximum" by way of taking the
+-- exponential of every value, summing the results, and then taking the
+-- logarithm. Incorporates some tricks to improve numerical stability.
+logSumExp :: (Ord x, Floating x, Traversable f) => f x -> x
+{-# INLINE logSumExp #-}
+logSumExp xs =
+    let mx = maximum xs
+     in (+ mx) . log1p . subtract 1 . sum $ exp . subtract mx <$> xs
+
+-- | Given a function, computes the "soft maximum" of the function by computing
+-- the integral of the exponential of the function, and taking the logarithm of
+-- the result. The maximum is first approximated on a given set of samples to
+-- improve numerical stability. Pro tip: If you want to compute the normalizer
+-- of a an exponential family probability density, provide the unnormalized
+-- log-density to this function.
+logIntegralExp
+    :: Traversable f
+    => Double -- ^ Error Tolerance
+    -> (Double -> Double) -- ^ Function
+    -> Double -- ^ Interval beginning
+    -> Double -- ^ Interval end
+    -> f Double -- ^ Samples (for approximating the max)
+    -> Double -- ^ Log-Integral-Exp
+{-# INLINE logIntegralExp #-}
+logIntegralExp err f mnbnd mxbnd xsmps =
+    let mx = maximum $ f <$> xsmps
+        expf x = exp $ f x - mx
+     in (+ mx) . log1p . subtract 1 . fst $ integrate err expf mnbnd mxbnd
+
+
+--- TypeLits ---
+
+
+-- | Type-level triangular number.
+type Triangular n = Div (n * (n + 1)) 2
+
+-- | Type level rational numbers. This implementation does not currently permit negative numbers.
+data Rat (n :: Nat) (d :: Nat)
+
+-- | Infix 'Rat'.
+type (/) n d = Rat n d
+
+-- | Recover a rational value from a 'Proxy'.
+ratVal :: (KnownNat n, KnownNat d) => Proxy (n / d) -> Rational
+{-# INLINE ratVal #-}
+ratVal = ratVal0 Proxy Proxy
+
+
+-- | 'natVal and 'fromIntegral'.
+natValInt :: KnownNat n => Proxy n -> Int
+{-# INLINE natValInt #-}
+natValInt = fromIntegral . natVal
+
+-- | 'getFinite' and 'fromIntegral'.
+finiteInt :: KnownNat n => Finite n -> Int
+{-# INLINE finiteInt #-}
+finiteInt = fromIntegral . getFinite
+
+ratVal0 :: (KnownNat n, KnownNat d) => Proxy n -> Proxy d -> Proxy (n / d) -> Rational
+{-# INLINE ratVal0 #-}
+ratVal0 prxyn prxyd _ = fromIntegral (natVal prxyn) % fromIntegral (natVal prxyd)
diff --git a/Goal/Core/Vector/Boxed.hs b/Goal/Core/Vector/Boxed.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Core/Vector/Boxed.hs
@@ -0,0 +1,245 @@
+ {-# OPTIONS_GHC -fplugin=GHC.TypeLits.KnownNat.Solver -fplugin=GHC.TypeLits.Normalise -fconstraint-solver-iterations=10 #-}
+-- | Vectors and Matrices with statically-typed dimensions based on boxed vectors.
+
+module Goal.Core.Vector.Boxed
+    ( -- * Vector
+      module Data.Vector.Sized
+      -- ** Construction
+    , doubleton
+    , range
+    , breakStream
+    , breakEvery
+    -- ** Deconstruction
+    , toPair
+    , concat
+    -- * Matrix
+    , Matrix
+    , nRows
+    , nColumns
+    -- ** Construction
+    , fromRows
+    , fromColumns
+    , matrixIdentity
+    , outerProduct
+    , diagonalConcat
+    -- ** Deconstruction
+    , toRows
+    , toColumns
+    -- ** Manipulation
+    , columnVector
+    , rowVector
+    -- ** BLAS
+    , dotProduct
+    , matrixVectorMultiply
+    , matrixMatrixMultiply
+    , inverse
+    , transpose
+    ) where
+
+
+--- Imports ---
+
+
+-- Goal --
+
+import Goal.Core.Util hiding (breakEvery,range)
+import qualified Goal.Core.Util (breakEvery)
+
+import qualified Goal.Core.Vector.Generic as G
+
+-- Unqualified --
+
+import Prelude hiding (concat,zipWith,(++),replicate)
+import qualified Data.Vector as B
+import qualified Data.Vector.Mutable as BM
+import qualified Control.Monad.ST as ST
+import qualified Data.Vector.Generic.Sized.Internal as I
+
+-- Qualified --
+
+import Data.Vector.Sized
+import GHC.TypeNats
+import Data.Proxy
+
+-- Qualified Imports --
+
+-- | Flatten a 'Vector' of 'Vector's.
+concat :: KnownNat n => Vector m (Vector n x) -> Vector (m*n) x
+{-# INLINE concat #-}
+concat = G.concat
+
+-- | Create a 'Vector' of length 2.
+doubleton :: x -> x -> Vector 2 x
+{-# INLINE doubleton #-}
+doubleton = G.doubleton
+
+-- | Partition of an interval.
+range :: (KnownNat n, Fractional x) => x -> x -> Vector n x
+{-# INLINE range #-}
+range = G.range
+
+-- | Cycles a list of elements and breaks it up into an infinite list of 'Vector's.
+breakStream :: forall n a. KnownNat n => [a] -> [Vector n a]
+{-# INLINE breakStream #-}
+breakStream as =
+    I.Vector . B.fromList <$> Goal.Core.Util.breakEvery (natValInt (Proxy :: Proxy n)) (cycle as)
+
+-- | Converts a length two 'Vector' into a pair of elements.
+toPair :: Vector 2 x -> (x,x)
+{-# INLINE toPair #-}
+toPair = G.toPair
+
+-- | Breaks a 'Vector' into a Vector of Vectors.
+breakEvery :: (KnownNat n, KnownNat k) => Vector (n*k) a -> Vector n (Vector k a)
+{-# INLINE breakEvery #-}
+breakEvery = G.breakEvery
+
+
+--- Matrices ---
+
+
+-- | Matrices with static dimensions (boxed).
+type Matrix = G.Matrix B.Vector
+
+-- | The number of rows in the 'Matrix'.
+nRows :: forall m n a . KnownNat m => Matrix m n a -> Int
+{-# INLINE nRows #-}
+nRows = G.nRows
+
+-- | The number of columns in the 'Matrix'.
+nColumns :: forall m n a . KnownNat n => Matrix m n a -> Int
+{-# INLINE nColumns #-}
+nColumns = G.nColumns
+
+-- | Convert a 'Matrix' into a 'Vector' of 'Vector's of rows.
+toRows :: (KnownNat m, KnownNat n) => Matrix m n x -> Vector m (Vector n x)
+{-# INLINE toRows #-}
+toRows = G.toRows
+
+-- | Convert a 'Matrix' into a 'Vector' of 'Vector's of columns.
+toColumns :: (KnownNat m, KnownNat n) => Matrix m n x -> Vector n (Vector m x)
+{-# INLINE toColumns #-}
+toColumns = G.toColumns
+
+
+-- | Turn a 'Vector' into a single column 'Matrix'.
+columnVector :: Vector n a -> Matrix n 1 a
+{-# INLINE columnVector #-}
+columnVector = G.columnVector
+
+-- | Turn a 'Vector' into a single row 'Matrix'.
+rowVector :: Vector n a -> Matrix 1 n a
+{-# INLINE rowVector #-}
+rowVector = G.rowVector
+
+-- | Create a 'Matrix' from a 'Vector' of row 'Vector's.
+fromRows :: KnownNat n => Vector m (Vector n x) -> Matrix m n x
+{-# INLINE fromRows #-}
+fromRows = G.fromRows
+
+-- | Create a 'Matrix' from a 'Vector' of column 'Vector's.
+fromColumns :: (KnownNat n, KnownNat m) => Vector n (Vector m x) -> Matrix m n x
+{-# INLINE fromColumns #-}
+fromColumns = G.fromColumns
+
+-- | Diagonally concatenate two matrices, padding the gaps with zeroes (pure implementation).
+diagonalConcat
+    :: (KnownNat n, KnownNat m, KnownNat o, KnownNat p, Num a)
+    => Matrix n m a -> Matrix o p a -> Matrix (n+o) (m+p) a
+{-# INLINE diagonalConcat #-}
+diagonalConcat mtx1 mtx2 =
+    let rws1 = (++ replicate 0) <$> toRows mtx1
+        rws2 = (replicate 0 ++) <$> toRows mtx2
+     in fromRows $ rws1 ++ rws2
+
+-- | Pure implementation of the dot product.
+dotProduct :: Num x => Vector n x -> Vector n x -> x
+{-# INLINE dotProduct #-}
+dotProduct = G.dotProduct
+
+-- | Pure implementation of the outer product.
+outerProduct
+    :: (KnownNat m, KnownNat n, Num x)
+    => Vector m x -> Vector n x -> Matrix m n x
+{-# INLINE outerProduct #-}
+outerProduct = G.outerProduct
+
+-- | Pure implementation of 'Matrix' transposition.
+transpose
+    :: (KnownNat m, KnownNat n, Num x)
+    => Matrix m n x -> Matrix n m x
+{-# INLINE transpose #-}
+transpose = G.transpose
+
+-- | Pure 'Matrix' x 'Vector' multiplication.
+matrixVectorMultiply
+    :: (KnownNat m, KnownNat n, Num x)
+    => Matrix m n x -> Vector n x -> Vector m x
+{-# INLINE matrixVectorMultiply #-}
+matrixVectorMultiply mtx = G.toVector . matrixMatrixMultiply mtx . columnVector
+
+-- | The identity 'Matrix'.
+matrixIdentity :: (KnownNat n, Num a) => Matrix n n a
+{-# INLINE matrixIdentity #-}
+matrixIdentity =
+    fromRows $ generate (\i -> generate (\j -> if finiteInt i == finiteInt j then 1 else 0))
+
+-- | Pure implementation of matrix inversion.
+inverse :: forall a n. (Fractional a, Ord a, KnownNat n) => Matrix n n a -> Maybe (Matrix n n a)
+{-# INLINE inverse #-}
+inverse mtx =
+    let rws = fromSized $ fromSized <$> zipWith (++) (toRows mtx) (toRows matrixIdentity)
+        n = natValInt (Proxy :: Proxy n)
+        rws' = B.foldM' eliminateRow rws $ B.generate n id
+     in G.Matrix . I.Vector . B.concatMap (B.drop n) <$> rws'
+
+-- | Pure 'Matrix' x 'Matrix' multiplication.
+matrixMatrixMultiply
+    :: forall m n o a. (KnownNat m, KnownNat n, KnownNat o, Num a)
+    => Matrix m n a -> Matrix n o a -> Matrix m o a
+{-# INLINE matrixMatrixMultiply #-}
+matrixMatrixMultiply (G.Matrix (I.Vector v)) wm =
+    let n = natValInt (Proxy :: Proxy n)
+        o = natValInt (Proxy :: Proxy o)
+        (G.Matrix (I.Vector w')) = G.transpose wm
+        f k = let (i,j) = divMod (finiteInt k) o
+                  slc1 = B.unsafeSlice (i*n) n v
+                  slc2 = B.unsafeSlice (j*n) n w'
+               in G.weakDotProduct slc1 slc2
+     in G.Matrix $ G.generate f
+
+
+--- Internal ---
+
+
+eliminateRow :: (Ord a, Fractional a) => B.Vector (B.Vector a) -> Int -> Maybe (B.Vector (B.Vector a))
+eliminateRow mtx k = do
+    mtx' <- pivotRow k mtx
+    return . nullifyRows k $ normalizePivot k mtx'
+
+pivotRow :: (Fractional a, Ord a) => Int -> B.Vector (B.Vector a) -> Maybe (B.Vector (B.Vector a))
+pivotRow k rws =
+    let l = (+k) . B.maxIndex $ abs . flip B.unsafeIndex k . B.take (B.length rws) <$> B.drop k rws
+        ak = B.unsafeIndex rws k B.! l
+     in if abs ak < 1e-10 then Nothing
+                  else ST.runST $ do
+                           mrws <- B.thaw rws
+                           BM.unsafeSwap mrws k l
+                           Just <$> B.freeze mrws
+
+normalizePivot :: Fractional a => Int -> B.Vector (B.Vector a) -> B.Vector (B.Vector a)
+normalizePivot k rws = ST.runST $ do
+    let ak = recip . flip B.unsafeIndex k $ B.unsafeIndex rws k
+    mrws <- B.thaw rws
+    BM.modify mrws ((*ak) <$>) k
+    B.freeze mrws
+
+nullifyRows :: Fractional a => Int -> B.Vector (B.Vector a) -> B.Vector (B.Vector a)
+nullifyRows k rws =
+    let rwk = B.unsafeIndex rws k
+        ak = B.unsafeIndex rwk k
+        generator i = if i == k then 0 else B.unsafeIndex (B.unsafeIndex rws i) k / ak
+        as = B.generate (B.length rws) generator
+     in B.zipWith (B.zipWith (-)) rws $ (\a -> (*a) <$> rwk) <$> as
+
+
diff --git a/Goal/Core/Vector/Generic.hs b/Goal/Core/Vector/Generic.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Core/Vector/Generic.hs
@@ -0,0 +1,233 @@
+{-# LANGUAGE StandaloneDeriving,GeneralizedNewtypeDeriving #-}
+ {-# OPTIONS_GHC -fplugin=GHC.TypeLits.KnownNat.Solver -fplugin=GHC.TypeLits.Normalise -fconstraint-solver-iterations=10 #-}
+-- | Vectors and Matrices with statically typed dimensions.
+module Goal.Core.Vector.Generic
+    ( -- * Vector
+      module Data.Vector.Generic.Sized
+    , VectorClass
+    -- * Construction
+    , doubleton
+    , range
+    , breakEvery
+    -- * Deconstruction
+    , concat
+    -- * Matrix
+    , Matrix (Matrix,toVector)
+    , nRows
+    , nColumns
+    -- ** Construction
+    , fromRows
+    , fromColumns
+    -- ** Deconstruction
+    , toPair
+    , toRows
+    , toColumns
+    -- ** Manipulation
+    , columnVector
+    , rowVector
+    -- ** BLAS
+    , transpose
+    , dotProduct
+    , weakDotProduct
+    , outerProduct
+    , matrixVectorMultiply
+    , matrixMatrixMultiply
+    ) where
+
+
+--- Imports ---
+
+
+-- Goal --
+
+import Goal.Core.Util hiding (breakEvery,range)
+
+-- Unqualified --
+
+import GHC.TypeNats
+import Data.Proxy
+import Control.DeepSeq
+import Data.Vector.Generic.Sized
+import Data.Vector.Generic.Sized.Internal
+import Foreign.Storable
+import Prelude hiding (concatMap,concat,map,sum,replicate)
+
+-- Qualified --
+
+import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Storable as S
+
+import Numeric.LinearAlgebra (Numeric)
+
+--- Vector ---
+
+
+type VectorClass = G.Vector
+
+-- | Create a 'Matrix' from a 'Vector' of 'Vector's which represent the rows.
+concat :: (KnownNat n, G.Vector v x, G.Vector v (Vector v n x)) => Vector v m (Vector v n x) -> Vector v (m*n) x
+{-# INLINE concat #-}
+concat = concatMap id
+
+-- | Collect two values into a length 2 'Vector'.
+doubleton :: G.Vector v x => x -> x -> Vector v 2 x
+{-# INLINE doubleton #-}
+doubleton x1 x2 = cons x1 $ singleton x2
+
+-- | Breaks a 'Vector' into a Vector of Vectors.
+breakEvery
+    :: forall v n k a . (G.Vector v a, G.Vector v (Vector v k a), KnownNat n, KnownNat k)
+    => Vector v (n*k) a -> Vector v n (Vector v k a)
+{-# INLINE breakEvery #-}
+breakEvery v0 =
+    let k = natValInt (Proxy :: Proxy k)
+        v = fromSized v0
+     in generate (\i -> Vector $ G.unsafeSlice (finiteInt i*k) k v)
+
+-- | Reshapes a length 2 'Vector' into a pair of values.
+toPair :: G.Vector v a => Vector v 2 a -> (a,a)
+{-# INLINE toPair #-}
+toPair v = (unsafeIndex v 0, unsafeIndex v 1)
+
+-- | Uniform partition of an interval into a 'Vector'.
+range
+    :: forall v n x. (G.Vector v x, KnownNat n, Fractional x)
+    => x -> x -> Vector v n x
+{-# INLINE range #-}
+range mn mx =
+    let n = natValInt (Proxy :: Proxy n)
+        stp = (mx - mn)/fromIntegral (n-1)
+     in enumFromStepN mn stp
+
+
+--- Matrix ---
+
+
+-- | Matrices with static dimensions.
+newtype Matrix v (m :: Nat) (n :: Nat) a = Matrix { toVector :: Vector v (m*n) a }
+    deriving (Eq,Show,NFData)
+
+deriving instance (KnownNat m, KnownNat n, Storable x) => Storable (Matrix S.Vector m n x)
+deriving instance (KnownNat m, KnownNat n, Numeric x, Num x)
+  => Num (Matrix S.Vector m n x)
+deriving instance (KnownNat m, KnownNat n, Numeric x, Fractional x)
+  => Fractional (Matrix S.Vector m n x)
+deriving instance (KnownNat m, KnownNat n, Numeric x, Floating x)
+  => Floating (Matrix S.Vector m n x)
+
+-- | Turn a 'Vector' into a single column 'Matrix'.
+columnVector :: Vector v n a -> Matrix v n 1 a
+{-# INLINE columnVector #-}
+columnVector = Matrix
+
+-- | Turn a 'Vector' into a single row 'Matrix'.
+rowVector :: Vector v n a -> Matrix v 1 n a
+{-# INLINE rowVector #-}
+rowVector = Matrix
+
+-- | Create a 'Matrix' from a 'Vector' of 'Vector's which represent the rows.
+fromRows :: (G.Vector v x, G.Vector v (Vector v n x), KnownNat n) => Vector v m (Vector v n x) -> Matrix v m n x
+{-# INLINE fromRows #-}
+fromRows = Matrix . concat
+
+-- | Create a 'Matrix' from a 'Vector' of 'Vector's which represent the columns.
+fromColumns
+    :: (G.Vector v x, G.Vector v Int, G.Vector v (Vector v n x), G.Vector v (Vector v m x), KnownNat n, KnownNat m)
+    => Vector v n (Vector v m x) -> Matrix v m n x
+{-# INLINE fromColumns #-}
+fromColumns = transpose . fromRows
+
+-- | The number of rows in the 'Matrix'.
+nRows :: forall v m n a . KnownNat m => Matrix v m n a -> Int
+{-# INLINE nRows #-}
+nRows _ = natValInt (Proxy :: Proxy m)
+
+-- | The number of columns in the 'Matrix'.
+nColumns :: forall v m n a . KnownNat n => Matrix v m n a -> Int
+{-# INLINE nColumns #-}
+nColumns _ = natValInt (Proxy :: Proxy n)
+
+-- | Convert a 'Matrix' into a 'Vector' of 'Vector's of rows.
+toRows :: (G.Vector v a, G.Vector v (Vector v n a), KnownNat n, KnownNat m)
+       => Matrix v m n a -> Vector v m (Vector v n a)
+{-# INLINE toRows #-}
+toRows (Matrix v) = breakEvery v
+
+-- | Convert a 'Matrix' into a 'Vector' of 'Vector's of columns.
+toColumns
+    :: (G.Vector v a, G.Vector v (Vector v m a), KnownNat m, KnownNat n, G.Vector v Int)
+    => Matrix v m n a -> Vector v n (Vector v m a)
+{-# INLINE toColumns #-}
+toColumns = toRows . transpose
+
+
+--- BLAS ---
+
+
+-- | Pure implementation of 'Matrix' transposition.
+transpose
+    :: forall v m n a . (KnownNat m, KnownNat n, G.Vector v Int, G.Vector v a, G.Vector v (Vector v m a))
+    => Matrix v m n a -> Matrix v n m a
+{-# INLINE transpose #-}
+transpose (Matrix v) =
+    let n = natValInt (Proxy :: Proxy n)
+     in fromRows $ generate (\j -> generate (\i -> unsafeIndex v $ finiteInt j + finiteInt i*n) :: Vector v m a)
+
+-- | Pure implementation of the dot product.
+dotProduct :: (G.Vector v x, Num x) => Vector v n x -> Vector v n x -> x
+{-# INLINE dotProduct #-}
+dotProduct v1 v2 = weakDotProduct (fromSized v1) (fromSized v2)
+
+-- | Pure implementation of the outer product.
+outerProduct
+    :: ( KnownNat m, KnownNat n, Num x
+       , G.Vector v Int, G.Vector v x, G.Vector v (Vector v n x), G.Vector v (Vector v m x), G.Vector v (Vector v 1 x) )
+     => Vector v n x -> Vector v m x -> Matrix v n m x
+{-# INLINE outerProduct #-}
+outerProduct v1 v2 = matrixMatrixMultiply (columnVector v1) (rowVector v2)
+
+-- | Pure implementation of the dot product on standard vectors.
+weakDotProduct :: (G.Vector v x, Num x) => v x -> v x -> x
+{-# INLINE weakDotProduct #-}
+weakDotProduct v1 v2 = G.foldl foldFun 0 (G.enumFromN 0 (G.length v1) :: S.Vector Int)
+    where foldFun d i = d + G.unsafeIndex v1 i * G.unsafeIndex v2 i
+
+-- | Pure 'Matrix' x 'Vector' multiplication.
+matrixVectorMultiply
+    :: (KnownNat m, KnownNat n, G.Vector v x, G.Vector v (Vector v n x), Num x)
+    => Matrix v m n x
+    -> Vector v n x
+    -> Vector v m x
+{-# INLINE matrixVectorMultiply #-}
+matrixVectorMultiply mtx v =
+    map (dotProduct v) $ toRows mtx
+
+-- | Pure 'Matrix' x 'Matrix' multiplication.
+matrixMatrixMultiply
+    :: ( KnownNat m, KnownNat n, KnownNat o, Num x
+       , G.Vector v Int, G.Vector v x, G.Vector v (Vector v m x), G.Vector v (Vector v n x), G.Vector v (Vector v o x) )
+    => Matrix v m n x
+    -> Matrix v n o x
+    -> Matrix v m o x
+{-# INLINE matrixMatrixMultiply #-}
+matrixMatrixMultiply mtx1 mtx2 =
+    fromColumns . map (matrixVectorMultiply mtx1) $ toColumns mtx2
+
+
+--- Numeric Classes ---
+
+
+--instance (Storable x, Numeric x, KnownNat n, KnownNat m)
+--  => Num (Matrix S.Vector n m x) where
+--    {-# INLINE (+) #-}
+--    (+) (Matrix (Vector v1)) (Matrix (Vector v2)) = Matrix $ Vector (H.add v1 v2)
+--    {-# INLINE (*) #-}
+--    (*) (Matrix xs) (Matrix xs') = Matrix $ xs * xs'
+--    {-# INLINE negate #-}
+--    negate (Matrix (Vector v)) = Matrix $ Vector (H.scale (-1) v)
+--    {-# INLINE abs #-}
+--    abs (Matrix xs) = Matrix $ abs xs
+--    {-# INLINE signum #-}
+--    signum (Matrix xs) = Matrix $ signum xs
+--    {-# INLINE fromInteger #-}
+--    fromInteger x = Matrix . replicate $ fromInteger x
diff --git a/Goal/Core/Vector/Generic/Internal.hs b/Goal/Core/Vector/Generic/Internal.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Core/Vector/Generic/Internal.hs
@@ -0,0 +1,11 @@
+-- | Forbidden access to the static vector constructor.
+module Goal.Core.Vector.Generic.Internal
+    ( -- * Vector
+      module Data.Vector.Generic.Sized.Internal
+    ) where
+
+
+--- Imports ---
+
+
+import Data.Vector.Generic.Sized.Internal
diff --git a/Goal/Core/Vector/Generic/Mutable.hs b/Goal/Core/Vector/Generic/Mutable.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Core/Vector/Generic/Mutable.hs
@@ -0,0 +1,11 @@
+-- | Mutable generic static vectors.
+module Goal.Core.Vector.Generic.Mutable
+    ( -- * Vector
+      module Data.Vector.Generic.Mutable.Sized
+    ) where
+
+
+--- Imports ---
+
+
+import Data.Vector.Generic.Mutable.Sized
diff --git a/Goal/Core/Vector/Storable.hs b/Goal/Core/Vector/Storable.hs
new file mode 100644
--- /dev/null
+++ b/Goal/Core/Vector/Storable.hs
@@ -0,0 +1,741 @@
+ {-# OPTIONS_GHC -fplugin=GHC.TypeLits.KnownNat.Solver -fplugin=GHC.TypeLits.Normalise -fconstraint-solver-iterations=10 #-}
+-- | Vectors and Matrices with statically typed dimensions based on storable vectors and using HMatrix where possible.
+module Goal.Core.Vector.Storable
+    ( -- * Vector
+      module Data.Vector.Storable.Sized
+      -- ** Construction
+    , doubleton
+    , range
+      -- ** Deconstruction
+    , concat
+    , breakEvery
+    , toPair
+    -- ** Computation
+    , average
+    , zipFold
+    -- * Matrix
+    , Matrix
+    , nRows
+    , nColumns
+    -- ** Construction
+    , fromRows
+    , fromColumns
+    , matrixIdentity
+    , outerProduct
+    , sumOuterProduct
+    , averageOuterProduct
+    , weightedAverageOuterProduct
+    , diagonalMatrix
+    , fromLowerTriangular
+    -- ** Deconstruction
+    , toRows
+    , toColumns
+    , lowerTriangular
+    , takeDiagonal
+    -- ** Manipulation
+    , columnVector
+    , rowVector
+    , combineTriangles
+    , diagonalConcat
+    , horizontalConcat
+    , verticalConcat
+    -- ** Computation
+    , trace
+    , withMatrix
+    -- *** BLAS
+    , scale
+    , add
+    , dotProduct
+    , dotMap
+    , matrixVectorMultiply
+    , matrixMatrixMultiply
+    , matrixMap
+    , eigens
+    , isSemiPositiveDefinite
+    , determinant
+    , inverseLogDeterminant
+    , inverse
+    , pseudoInverse
+    , matrixRoot
+    , transpose
+    -- *** Least Squares
+    , linearLeastSquares
+    , meanSquaredError
+    , rSquared
+    , l2Norm
+    , unsafeCholesky
+    -- *** Convolutions
+    , crossCorrelate2d
+    , convolve2d
+    , kernelOuterProduct
+    , kernelTranspose
+    -- ** Miscellaneous
+    , prettyPrintMatrix
+    ) where
+
+
+--- Imports ---
+
+
+-- Goal --
+
+import Goal.Core.Util hiding (average,breakEvery,range)
+
+-- Unqualified --
+
+import Data.Proxy
+import Data.Complex
+import Foreign.Storable
+import Data.Vector.Storable.Sized
+import Numeric.LinearAlgebra (Field,Numeric)
+import GHC.TypeNats
+import Prelude hiding (concat,foldr1,concatMap,replicate,(++),length,map,sum,zip,and)
+
+-- Qualified --
+
+import qualified Prelude
+import qualified Data.Vector.Storable as S
+import qualified Goal.Core.Vector.Generic as G
+import qualified Data.Vector.Generic.Sized.Internal as G
+import qualified Numeric.LinearAlgebra as H
+import qualified Data.List as L
+
+
+--- Generic ---
+
+
+-- | Matrices with static dimensions (storable).
+type Matrix = G.Matrix S.Vector
+
+-- | A fold over pairs of elements of 'Vector's of equal length.
+zipFold :: (KnownNat n, Storable x, Storable y)
+        => (z -> x -> y -> z)
+        -> z
+        -> Vector n x
+        -> Vector n y
+        -> z
+{-# INLINE zipFold #-}
+zipFold f z0 xs ys =
+    let n = length xs
+        foldfun z i = f z (unsafeIndex xs i) (unsafeIndex ys i)
+     in L.foldl' foldfun z0 [0..n-1]
+
+-- | Concatenates a vector of vectors into a single vector.
+concat :: (KnownNat n, Storable x) => Vector m (Vector n x) -> Vector (m*n) x
+{-# INLINE concat #-}
+concat = G.concat
+
+-- | Collect two values into a length 2 'Vector'.
+doubleton :: Storable x => x -> x -> Vector 2 x
+{-# INLINE doubleton #-}
+doubleton = G.doubleton
+
+-- | The number of rows in the 'Matrix'.
+nRows :: forall m n a . KnownNat m => Matrix m n a -> Int
+{-# INLINE nRows #-}
+nRows = G.nRows
+
+-- | The columns of columns in the 'Matrix'.
+nColumns :: forall m n a . KnownNat n => Matrix m n a -> Int
+{-# INLINE nColumns #-}
+nColumns = G.nColumns
+
+-- | Convert a 'Matrix' into a 'Vector' of 'Vector's of rows.
+toRows :: (KnownNat m, KnownNat n, Storable x) => Matrix m n x -> Vector m (Vector n x)
+{-# INLINE toRows #-}
+toRows = G.toRows
+
+-- | Turn a 'Vector' into a single column 'Matrix'.
+columnVector :: Vector n a -> Matrix n 1 a
+{-# INLINE columnVector #-}
+columnVector = G.columnVector
+
+-- | Turn a 'Vector' into a single row 'Matrix'.
+rowVector :: Vector n a -> Matrix 1 n a
+{-# INLINE rowVector #-}
+rowVector = G.rowVector
+
+-- | Create a 'Matrix' from a 'Vector' of 'Vector's which represent the rows.
+fromRows :: (KnownNat n, Storable x) => Vector m (Vector n x) -> Matrix m n x
+{-# INLINE fromRows #-}
+fromRows = G.fromRows
+
+-- | Uniform partition of an interval into a 'Vector'.
+range :: (KnownNat n, Fractional x, Storable x) => x -> x -> Vector n x
+{-# INLINE range #-}
+range = G.range
+
+-- | Reshapes a length 2 'Vector' into a pair of values.
+toPair :: Storable x => Vector 2 x -> (x,x)
+{-# INLINE toPair #-}
+toPair = G.toPair
+
+
+--- HMatrix ---
+
+
+-- | Converts a pure, Storable-based 'Matrix' into an HMatrix matrix.
+toHMatrix
+    :: forall m n x . (KnownNat n, KnownNat m, H.Element x, Storable x)
+    => Matrix m n x
+    -> H.Matrix x
+{-# INLINE toHMatrix #-}
+toHMatrix (G.Matrix mtx) =
+    let n = natValInt (Proxy :: Proxy n)
+        m = natValInt (Proxy :: Proxy m)
+     in if n == 0
+           then H.fromRows $ Prelude.replicate m S.empty
+           else H.reshape n $ fromSized mtx
+
+-- | Converts an HMatrix matrix into a pure, Storable-based 'Matrix'.
+fromHMatrix :: Numeric x => H.Matrix x -> Matrix m n x
+{-# INLINE fromHMatrix #-}
+fromHMatrix = G.Matrix . G.Vector . H.flatten
+
+-- | Convert a 'Matrix' into a 'Vector' of 'Vector's of columns.
+toColumns :: (KnownNat m, KnownNat n, Numeric x) => Matrix m n x -> Vector n (Vector m x)
+{-# INLINE toColumns #-}
+toColumns = toRows . transpose
+
+-- | Create a 'Matrix' from a 'Vector' of 'Vector's which represent the columns.
+fromColumns :: (KnownNat m, KnownNat n, Numeric x) => Vector n (Vector m x) -> Matrix m n x
+{-# INLINE fromColumns #-}
+fromColumns = transpose . fromRows
+
+-- | Breaks a 'Vector' into a Vector of Vectors.
+breakEvery :: forall n k a . (KnownNat n, KnownNat k, Storable a) => Vector (n*k) a -> Vector n (Vector k a)
+{-# INLINE breakEvery #-}
+breakEvery v0 =
+    let k = natValInt (Proxy :: Proxy k)
+        v = fromSized v0
+     in generate (\i -> G.Vector $ S.unsafeSlice (finiteInt i*k) k v)
+
+
+--- BLAS ---
+
+
+-- | The sum of two 'Vector's.
+add :: Numeric x => Vector n x -> Vector n x -> Vector n x
+{-# INLINE add #-}
+add (G.Vector v1) (G.Vector v2) = G.Vector (H.add v1 v2)
+
+-- | Scalar multiplication of a 'Vector'.
+scale :: Numeric x => x -> Vector n x -> Vector n x
+{-# INLINE scale #-}
+scale x (G.Vector v) = G.Vector (H.scale x v)
+
+-- | Apply a 'Vector' operation to a 'Matrix'.
+withMatrix :: (Vector (n*m) x -> Vector (n*m) x) -> Matrix n m x -> Matrix n m x
+{-# INLINE withMatrix #-}
+withMatrix f (G.Matrix v) = G.Matrix $ f v
+
+-- | Returns the lower triangular part of a square matrix.
+lowerTriangular :: forall n x . (Storable x, H.Element x, KnownNat n) => Matrix n n x -> Vector (Triangular n) x
+{-# INLINE lowerTriangular #-}
+lowerTriangular mtx =
+    let hmtx = toHMatrix mtx
+        rws = H.toRows hmtx
+        rws' = Prelude.zipWith S.take [1..] rws
+     in G.Vector $ S.concat rws'
+--    let n = natValInt (Proxy :: Proxy n)
+--        idxs = G.Vector . S.fromList
+--            $ Prelude.concat [ from2Index n <$> Prelude.zip (repeat k) [0..k] | k <- [0..n-1] ]
+--     in backpermute xs idxs
+
+-- | Constructs a `Matrix` from a lower triangular part.
+fromLowerTriangular :: forall n x . (Storable x, KnownNat n) => Vector (Triangular n) x -> Matrix n n x
+{-# INLINE fromLowerTriangular #-}
+fromLowerTriangular xs =
+    let n = natValInt (Proxy :: Proxy n)
+        idxs = generate (toTriangularIndex . to2Index n . finiteInt)
+     in G.Matrix $ backpermute xs idxs
+
+-- | Build a matrix with the given diagonal, lower triangular part given by the
+-- first matrix, and upper triangular part given by the second matrix.
+combineTriangles
+    :: (KnownNat k, Storable x)
+    => Vector k x -- ^ Diagonal
+    -> Matrix k k x -- ^ Lower triangular source
+    -> Matrix k k x -- ^ Upper triangular source
+    -> Matrix k k x
+{-# INLINE combineTriangles #-}
+combineTriangles (G.Vector diag) crs1 crs2 =
+    fromRows $ generate (generator (toRows crs1) (toRows crs2))
+        where
+            generator rws1 rws2 fnt =
+                let (G.Vector rw1) = index rws1 fnt
+                    (G.Vector rw2) = index rws2 fnt
+                    i = fromIntegral fnt
+                 in G.Vector $ S.take i rw1 S.++ S.cons (diag S.! i) (S.drop (i+1) rw2)
+
+-- | The average of a 'Vector' of elements.
+average :: (Numeric x, Fractional x) => Vector n x -> x
+{-# INLINE average #-}
+average (G.Vector v) = H.sumElements v / fromIntegral (S.length v)
+
+-- | The dot product of two numerical 'Vector's.
+dotProduct :: Numeric x => Vector n x -> Vector n x -> x
+{-# INLINE dotProduct #-}
+dotProduct v1 v2 = H.dot (fromSized v1) (fromSized v2)
+
+-- | The determinant of a 'Matrix'.
+diagonalMatrix :: forall n x . (KnownNat n, Field x) => Vector n x -> Matrix n n x
+{-# INLINE diagonalMatrix #-}
+diagonalMatrix v =
+    let n = natValInt (Proxy :: Proxy n)
+     in fromHMatrix $ H.diagRect 0 (fromSized v) n n
+
+-- | The determinant of a 'Matrix'.
+takeDiagonal :: (KnownNat n, Field x) => Matrix n n x -> Vector n x
+{-# INLINE takeDiagonal #-}
+takeDiagonal = G.Vector . H.takeDiag . toHMatrix
+
+-- | The determinant of a 'Matrix'.
+trace :: (KnownNat n, Field x) => Matrix n n x -> x
+{-# INLINE trace #-}
+trace = S.sum . H.takeDiag . toHMatrix
+
+-- | Returns the eigenvalues and eigenvectors 'Matrix'.
+eigens :: (KnownNat n, Field x) => Matrix n n x -> (Vector n (Complex Double), Vector n (Vector n (Complex Double)))
+{-# INLINE eigens #-}
+eigens mtx =
+    let (exs,evs) = H.eig $ toHMatrix mtx
+     in (G.Vector exs, G.Vector . S.fromList $ G.Vector <$> H.toColumns evs)
+
+-- | Test if the matrix is semi-positive definite.
+isSemiPositiveDefinite :: (KnownNat n, Field x) => Matrix n n x -> Bool
+{-# INLINE isSemiPositiveDefinite #-}
+isSemiPositiveDefinite =
+    and . map ((0 <=) . realPart) . fst . eigens
+
+-- | Returns the inverse, the logarithm of the absolute value of the
+-- determinant, and the sign of the determinant of a given matrix.
+inverseLogDeterminant :: (KnownNat n, Field x) => Matrix n n x -> (Matrix n n x, x, x)
+{-# INLINE inverseLogDeterminant #-}
+inverseLogDeterminant mtx =
+    let (imtx,(ldet,sgn)) = H.invlndet $ toHMatrix mtx
+     in (fromHMatrix imtx, ldet, sgn)
+
+-- | The determinant of a 'Matrix'.
+determinant :: (KnownNat n, Field x) => Matrix n n x -> x
+{-# INLINE determinant #-}
+determinant = H.det . toHMatrix
+
+-- | Transpose a 'Matrix'.
+transpose
+    :: forall m n x . (KnownNat m, KnownNat n, Numeric x)
+    => Matrix m n x
+    -> Matrix n m x
+{-# INLINE transpose #-}
+transpose (G.Matrix mtx) =
+    G.Matrix $ withVectorUnsafe (H.flatten . H.tr . H.reshape (natValInt (Proxy :: Proxy n))) mtx
+
+-- | Diagonally concatenate two matrices, padding the gaps with zeroes.
+diagonalConcat
+    :: (KnownNat n, KnownNat m, KnownNat o, KnownNat p, Numeric x)
+    => Matrix n m x -> Matrix o p x -> Matrix (n+o) (m+p) x
+{-# INLINE diagonalConcat #-}
+diagonalConcat mtx10 mtx20 =
+    let mtx1 = toHMatrix mtx10
+        mtx2 = toHMatrix mtx20
+     in fromHMatrix $ H.diagBlock [mtx1,mtx2]
+
+-- | Diagonally concatenate two matrices, padding the gaps with zeroes.
+horizontalConcat
+    :: (KnownNat n, KnownNat m, KnownNat o, Numeric x)
+    => Matrix n m x -> Matrix n o x -> Matrix n (m+o) x
+{-# INLINE horizontalConcat #-}
+horizontalConcat mtx10 mtx20 =
+    let mtx1 = toHMatrix mtx10
+        mtx2 = toHMatrix mtx20
+     in fromHMatrix $ mtx1 H.||| mtx2
+
+-- | Diagonally concatenate two matrices, padding the gaps with zeroes.
+verticalConcat
+    :: (KnownNat n, KnownNat m, KnownNat o, Numeric x)
+    => Matrix n o x -> Matrix m o x -> Matrix (n+m) o x
+{-# INLINE verticalConcat #-}
+verticalConcat mtx10 mtx20 =
+    let mtx1 = toHMatrix mtx10
+        mtx2 = toHMatrix mtx20
+     in fromHMatrix $ mtx1 H.=== mtx2
+
+-- | Invert a 'Matrix'.
+inverse :: forall n x . (KnownNat n, Field x) => Matrix n n x -> Matrix n n x
+{-# INLINE inverse #-}
+inverse (G.Matrix mtx) =
+    G.Matrix $ withVectorUnsafe (H.flatten . H.inv . H.reshape (natValInt (Proxy :: Proxy n))) mtx
+
+-- | Pseudo-Invert a 'Matrix'.
+pseudoInverse :: forall n x . (KnownNat n, Field x) => Matrix n n x -> Matrix n n x
+{-# INLINE pseudoInverse #-}
+pseudoInverse (G.Matrix mtx) =
+    G.Matrix $ withVectorUnsafe (H.flatten . H.pinv . H.reshape (natValInt (Proxy :: Proxy n))) mtx
+
+-- | Square root of a 'Matrix'.
+matrixRoot :: forall n x . (KnownNat n, Field x) => Matrix n n x -> Matrix n n x
+{-# INLINE matrixRoot #-}
+matrixRoot (G.Matrix mtx) =
+    G.Matrix $ withVectorUnsafe (H.flatten . H.sqrtm . H.reshape (natValInt (Proxy :: Proxy n))) mtx
+
+-- | The outer product of two 'Vector's.
+outerProduct :: (KnownNat m, KnownNat n, Numeric x) => Vector m x -> Vector n x -> Matrix m n x
+{-# INLINE outerProduct #-}
+outerProduct v1 v2 =
+    fromHMatrix $ H.outer (fromSized v1) (fromSized v2)
+
+-- | The summed outer product of two lists of 'Vector's.
+sumOuterProduct :: (KnownNat m, KnownNat n, Fractional x, Numeric x) => [(Vector m x,Vector n x)] -> Matrix m n x
+{-# INLINE sumOuterProduct #-}
+sumOuterProduct v12s =
+    let (v1s,v2s) = L.unzip v12s
+        mtx1 = H.fromColumns $ fromSized <$> v1s
+        mtx2 = H.fromRows $ fromSized <$> v2s
+     in fromHMatrix (mtx1 H.<> mtx2)
+
+-- | The average outer product of two lists of 'Vector's.
+averageOuterProduct :: (KnownNat m, KnownNat n, Fractional x, Numeric x) => [(Vector m x,Vector n x)] -> Matrix m n x
+{-# INLINE averageOuterProduct #-}
+averageOuterProduct v12s =
+    let (v1s,v2s) = L.unzip v12s
+        mtx1 = H.fromColumns $ fromSized <$> v1s
+        (_,n) = H.size mtx1
+        mtx2 = H.scale (1/fromIntegral n) . H.fromRows $ fromSized <$> v2s
+     in fromHMatrix (mtx1 H.<> mtx2)
+
+-- | The average outer product of two lists of 'Vector's.
+weightedAverageOuterProduct
+    :: ( KnownNat m, KnownNat n, Fractional x, Numeric x )
+    => [(x,Vector m x,Vector n x)]
+    -> Matrix m n x
+{-# INLINE weightedAverageOuterProduct #-}
+weightedAverageOuterProduct wv12s =
+    let (ws,v1s,v2s) = L.unzip3 wv12s
+        v1s' = L.zipWith H.scale ws $ fromSized <$> v1s
+        mtx1 = H.fromColumns v1s'
+        mtx2 = H.fromRows $ fromSized <$> v2s
+     in fromHMatrix (mtx1 H.<> mtx2)
+
+-- | The identity 'Matrix'.
+matrixIdentity :: forall n x . (KnownNat n, Numeric x, Num x) => Matrix n n x
+{-# INLINE matrixIdentity #-}
+matrixIdentity =
+    fromHMatrix . H.ident $ natValInt (Proxy :: Proxy n)
+
+-- | The dot products of one vector with a list of vectors.
+dotMap :: (KnownNat n, Numeric x) => Vector n x -> [Vector n x] -> [x]
+{-# INLINE dotMap #-}
+dotMap v vs =
+    let mtx' = H.fromRows $ fromSized <$> vs
+     in H.toList $ mtx' H.#> fromSized v
+--     in if S.null w
+--           then replicate 0
+--           else fmap G.Vector . H.toColumns $ toHMatrix mtx H.<> mtx'
+
+-- | Map a linear transformation over a list of 'Vector's.
+matrixMap :: (KnownNat m, KnownNat n, Numeric x)
+                     => Matrix m n x -> [Vector n x] -> [Vector m x]
+{-# INLINE matrixMap #-}
+matrixMap mtx vs =
+    let mtx' = H.fromColumns $ fromSized <$> vs
+     in fmap G.Vector . H.toColumns $ toHMatrix mtx H.<> mtx'
+--     in if S.null w
+--           then replicate 0
+--           else fmap G.Vector . H.toColumns $ toHMatrix mtx H.<> mtx'
+
+
+-- | Apply a linear transformation to a 'Vector'.
+matrixVectorMultiply :: (KnownNat m, KnownNat n, Numeric x)
+                     => Matrix m n x -> Vector n x -> Vector m x
+{-# INLINE matrixVectorMultiply #-}
+matrixVectorMultiply mtx v =
+    G.Vector $ toHMatrix mtx H.#> fromSized v
+--    let w = toHMatrix mtx H.#> fromSized v
+--     in if S.null w
+--           then replicate 0
+--           else G.Vector w
+
+-- | Multiply a 'Matrix' with a second 'Matrix'.
+matrixMatrixMultiply
+    :: (KnownNat m, KnownNat n, KnownNat o, Numeric x)
+    => Matrix m n x
+    -> Matrix n o x
+    -> Matrix m o x
+{-# INLINE matrixMatrixMultiply #-}
+matrixMatrixMultiply mtx1 mtx2 = fromHMatrix $ toHMatrix mtx1 H.<> toHMatrix mtx2
+
+-- | Pretty print the values of a 'Matrix' (for extremely simple values of pretty).
+prettyPrintMatrix :: (KnownNat m, KnownNat n, Numeric a, Show a) => Matrix m n a -> IO ()
+prettyPrintMatrix = print . toHMatrix
+
+-- | The Mean Squared difference between two vectors.
+meanSquaredError
+    :: KnownNat k
+    => Vector k Double
+    -> Vector k Double
+    -> Double
+{-# INLINE meanSquaredError #-}
+meanSquaredError ys yhts = average $ map square (ys - yhts)
+
+-- | L2 length of a vector.
+l2Norm
+    :: KnownNat k
+    => Vector k Double
+    -> Double
+{-# INLINE l2Norm #-}
+l2Norm (G.Vector xs) = H.norm_2 xs
+
+-- | Computes the coefficient of determintation for the given outputs and model
+-- predictions.
+rSquared
+    :: KnownNat k
+    => Vector k Double -- ^ Dependent variable observations
+    -> Vector k Double -- ^ Predicted Values
+    -> Double -- ^ R-squared
+{-# INLINE rSquared #-}
+rSquared ys yhts =
+    let ybr = average ys
+        ssres = sum $ map square (ys - yhts)
+        sstot = sum $ map (square . subtract ybr) ys
+     in 1 - (ssres/sstot)
+
+-- | Solves the linear least squares problem.
+linearLeastSquares
+    :: KnownNat l
+    => [Vector l Double] -- ^ Independent variable observations
+    -> [Double] -- ^ Dependent variable observations
+    -> Vector l Double -- ^ Parameter estimates
+{-# INLINE linearLeastSquares #-}
+linearLeastSquares as xs =
+    G.Vector $ H.fromRows (fromSized <$> as) H.<\> S.fromList xs
+
+
+unsafeCholesky
+    :: (KnownNat n, Field x, Storable x)
+    => Matrix n n x
+    -> Matrix n n x
+unsafeCholesky =
+    transpose . fromHMatrix . H.chol . H.trustSym . toHMatrix
+
+
+--- Convolutions ---
+
+
+-- | 2d cross-correlation of a kernel over a matrix of values.
+crossCorrelate2d
+    :: forall nk rdkr rdkc mr mc md x
+    . ( KnownNat rdkr, KnownNat rdkc, KnownNat md, KnownNat mr, KnownNat mc
+      , KnownNat nk, Numeric x, Storable x )
+      => Proxy rdkr -- ^ Number of Kernel rows
+      -> Proxy rdkc -- ^ Number of Kernel columns
+      -> Proxy mr -- ^ Number of Matrix/Image rows
+      -> Proxy mc -- ^ Number of Kernel/Image columns
+      -> Matrix nk (md*(2*rdkr+1)*(2*rdkc+1)) x -- ^ Kernels (nk is their number)
+      -> Matrix md (mr*mc) x -- ^ Image (md is the depth)
+      -> Matrix nk (mr*mc) x -- ^ Cross-correlated image
+{-# INLINE crossCorrelate2d #-}
+crossCorrelate2d prdkr prdkc pmr pmc krns (G.Matrix v) =
+    let pmd = Proxy :: Proxy md
+        mtx = im2col prdkr prdkc pmd pmr pmc v
+     in matrixMatrixMultiply krns mtx
+
+-- | The transpose of a convolutional kernel.
+kernelTranspose
+    :: (KnownNat nk, KnownNat md, KnownNat rdkr, KnownNat rdkc, Numeric x, Storable x)
+    => Proxy nk
+    -> Proxy md
+    -> Proxy rdkr
+    -> Proxy rdkc
+    -> Matrix nk (md*(2*rdkr+1)*(2*rdkc+1)) x -- ^ Kernels (nk is their number)
+    -> Matrix md (nk*(2*rdkr+1)*(2*rdkc+1)) x -- ^ Kernels (nk is their number)
+{-# INLINE kernelTranspose #-}
+kernelTranspose pnk pmd prdkr prdkc (G.Matrix kv) = G.Matrix . backpermute kv $ kernelTransposeIndices pnk pmd prdkr prdkc
+
+-- | 2d convolution of a kernel over a matrix of values. This is the adjoint of crossCorrelate2d.
+convolve2d
+    :: forall nk rdkr rdkc md mr mc x
+    . ( KnownNat rdkr, KnownNat rdkc, KnownNat mr, KnownNat mc
+      , KnownNat md, KnownNat nk, Numeric x, Storable x )
+      => Proxy rdkr -- ^ Number of Kernel rows
+      -> Proxy rdkc -- ^ Number of Kernel columns
+      -> Proxy mr -- ^ Number of Matrix/Image rows
+      -> Proxy mc -- ^ Number of Kernel/Image columns
+      -> Matrix nk (md*(2*rdkr+1)*(2*rdkc+1)) x -- ^ Kernels (nk is their number)
+      -> Matrix nk (mr*mc) x -- ^ Dual image (nk is its depth)
+      -> Matrix md (mr*mc) x -- ^ Convolved image
+{-# INLINE convolve2d #-}
+convolve2d prdkr prdkc pmr pmc krn mtxs =
+    let pnk = Proxy :: Proxy nk
+        pmd = Proxy :: Proxy md
+        krn' = kernelTranspose pnk pmd prdkr prdkc krn
+     in crossCorrelate2d prdkr prdkc pmr pmc krn' mtxs
+
+-- | The outer product of an image and a dual image to produce a convolutional kernel.
+kernelOuterProduct
+    :: forall nk rdkr rdkc md mr mc x
+    . ( KnownNat rdkr, KnownNat rdkc, KnownNat mr, KnownNat mc
+      , KnownNat md, KnownNat nk, Numeric x, Storable x )
+      => Proxy rdkr -- ^ Number of Kernel rows
+      -> Proxy rdkc -- ^ Number of Kernel columns
+      -> Proxy mr -- ^ Number of Matrix/Image rows
+      -> Proxy mc -- ^ Number of Kernel/Image columns
+      -> Matrix nk (mr*mc) x -- ^ Dual image (nk is its depth)
+      -> Matrix md (mr*mc) x -- ^ Image (md is the depth)
+      -> Matrix nk (md*(2*rdkr+1)*(2*rdkc+1)) x -- ^ Kernels
+{-# INLINE kernelOuterProduct #-}
+kernelOuterProduct prdkr prdkc pmr pmc omtx (G.Matrix v) =
+    let pmd = Proxy :: Proxy md
+        imtx = im2col prdkr prdkc pmd pmr pmc v
+     in matrixMatrixMultiply omtx $ transpose imtx
+
+
+--- Internal ---
+
+
+toTriangularIndex :: (Int,Int) -> Int
+toTriangularIndex (i,j)
+    | i >= j = triangularNumber i + j
+    | otherwise = toTriangularIndex (j,i)
+
+to2Index :: Int -> Int -> (Int,Int)
+to2Index nj ij = divMod ij nj
+
+to3Index :: Int -> Int -> Int -> (Int,Int,Int)
+{-# INLINE to3Index #-}
+to3Index nj nk ijk =
+    let nj' = nj*nk
+        (i,jk) = divMod ijk nj'
+        (j,k) = divMod jk nk
+     in (i,j,k)
+
+from3Index :: Int -> Int -> (Int,Int,Int) -> Int
+{-# INLINE from3Index #-}
+from3Index nj nk (i,j,k) =
+    let nj' = nj*nk
+     in i*nj' + j*nk + k
+
+to4Index :: Int -> Int -> Int -> Int -> (Int,Int,Int,Int)
+{-# INLINE to4Index #-}
+to4Index nj nk nl ijkl =
+    let nk' = nl*nk
+        nj' = nj*nk'
+        (i,jkl) = divMod ijkl nj'
+        (j,kl) = divMod jkl nk'
+        (k,l) = divMod kl nl
+     in (i,j,k,l)
+
+from4Index :: Int -> Int -> Int -> (Int,Int,Int,Int) -> Int
+{-# INLINE from4Index #-}
+from4Index nj nk nl (i,j,k,l) =
+    let nk' = nl*nk
+        nj' = nj*nk'
+     in i*nj' + j*nk' + k*nl + l
+
+kernelTransposeIndices
+    :: (KnownNat nk, KnownNat md, KnownNat rdkr, KnownNat rdkc)
+    => Proxy nk
+    -> Proxy md
+    -> Proxy rdkr
+    -> Proxy rdkc
+    -> Vector (nk*md*(2*rdkr+1)*(2*rdkc+1)) Int
+{-# INLINE kernelTransposeIndices #-}
+kernelTransposeIndices pnk pmd prdkr prdkc =
+    let nkrn = natValInt pnk
+        md = natValInt pmd
+        rdkr = natValInt prdkr
+        rdkc = natValInt prdkc
+        dmkr = 2*rdkr+1
+        dmkc = 2*rdkc+1
+        nl = dmkc
+        nk = dmkr
+        nj = nkrn
+        nl' = dmkc
+        nk' = dmkr
+        nj' = md
+        reIndex idx =
+            let (i,j,k,l) = to4Index nj nk nl idx
+             in from4Index nj' nk' nl' (j,i,nk-1-k,nl-1-l)
+     in generate (reIndex . fromIntegral)
+
+im2colIndices
+    :: forall rdkr rdkc mr mc md
+     . (KnownNat rdkr, KnownNat rdkc, KnownNat mr, KnownNat mc, KnownNat md)
+    => Proxy rdkr
+    -> Proxy rdkc
+    -> Proxy md
+    -> Proxy mr
+    -> Proxy mc
+    -> Vector (((2*rdkr+1)*(2*rdkc+1)*md)*(mr*mc)) Int
+{-# INLINE im2colIndices #-}
+im2colIndices prdkr prdkc _ pmr pmc =
+    let rdkr = natValInt prdkr
+        rdkc = natValInt prdkc
+        nj = (2*rdkr + 1)
+        nk = (2*rdkc + 1)
+        reWindow idx =
+            let (i,j,k) = to3Index nj nk idx
+             in windowIndices prdkr prdkc pmr pmc i j k
+          in (concatMap reWindow :: Vector ((2*rdkr+1)*(2*rdkc+1)*md) Int -> Vector (((2*rdkr+1)*(2*rdkc+1)*md)*(mr*mc)) Int) $ generate finiteInt
+
+im2col
+    :: forall rdkr rdkc md mr mc x
+    . (KnownNat rdkr, KnownNat rdkc, KnownNat mc, KnownNat md, KnownNat mr, Num x, Storable x)
+    => Proxy rdkr
+    -> Proxy rdkc
+    -> Proxy md
+    -> Proxy mr
+    -> Proxy mc
+    -> Vector (md*mr*mc) x
+    -> Matrix (md*(2*rdkr+1)*(2*rdkc+1)) (mr*mc) x
+{-# INLINE im2col #-}
+im2col prdkr prdkc pmd pmr pmc mtx =
+    let idxs = im2colIndices prdkr prdkc pmd pmr pmc
+        mtx' = padMatrix prdkr prdkc pmd pmr pmc mtx
+     in G.Matrix $ backpermute mtx' idxs
+
+windowIndices
+    :: forall rdkr rdkc mr mc . (KnownNat rdkr, KnownNat rdkc, KnownNat mr, KnownNat mc)
+    => Proxy rdkr
+    -> Proxy rdkc
+    -> Proxy mr
+    -> Proxy mc
+    -> Int
+    -> Int
+    -> Int
+    -> Vector (mr*mc) Int
+{-# INLINE windowIndices #-}
+windowIndices prdkr prdkc pmr pmc kd kr kc =
+    let rdkr = natValInt prdkr
+        rdkc = natValInt prdkc
+        mr = natValInt pmr
+        mc = natValInt pmc
+        mrc = mr*mc
+        nj' = mr + 2*rdkr
+        nk' = mc + 2*rdkc
+        reIndex idx =
+            let (j,k) = divMod idx mc
+             in from3Index nj' nk' (kd,j+kr,k+kc)
+     in G.Vector $ S.generate mrc reIndex
+
+padMatrix
+    :: forall rdkr rdkc mr mc md x
+    . (KnownNat rdkr, KnownNat rdkc, KnownNat md, KnownNat mr, KnownNat mc, Num x, Storable x)
+    => Proxy rdkr
+    -> Proxy rdkc
+    -> Proxy md
+    -> Proxy mr
+    -> Proxy mc
+    -> Vector (md*mr*mc) x
+    -> Vector (md*(mr + 2*rdkr)*(mc + 2*rdkc)) x
+{-# INLINE padMatrix #-}
+padMatrix _ _ _ _ _ v =
+    let mtxs :: Vector md (Matrix mr mc x)
+        mtxs = map G.Matrix $ breakEvery v
+        pdrs :: Vector rdkr (Vector mc x)
+        pdrs = replicate $ replicate 0
+        mtxs' = map (\mtx -> fromRows $ pdrs ++ toRows mtx ++ pdrs) mtxs
+        pdcs :: Vector rdkc (Vector (mr + 2*rdkr) x)
+        pdcs = replicate $ replicate 0
+     in concatMap G.toVector $ map (\mtx' -> G.fromColumns $ pdcs ++ G.toColumns mtx' ++ pdcs) mtxs'
+
+
diff --git a/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,12 @@
+This is the least interesting package in the Goal libraries, and serves simply
+to re-export existing libraries and provide essential utility factors in a
+manner that is compatible with Goal. Nevertheless, there are a few modules worth
+highlighting.
+
+**Goal.Core.Circuit**: Provides an implementation of monadic Mealy automata to
+facilitate simple stream-based processing.
+
+**Goal.Core.Vector.Storable**: Combines the
+[vector-sized](https://hackage.haskell.org/package/vector-sized) and
+[hmatrix](https://hackage.haskell.org/package/hmatrix) libraries to provide
+efficient linear algebra with static sizes.
diff --git a/benchmarks/convolutions.hs b/benchmarks/convolutions.hs
new file mode 100644
--- /dev/null
+++ b/benchmarks/convolutions.hs
@@ -0,0 +1,110 @@
+{-# LANGUAGE
+    TypeOperators,
+    NoStarIsType,
+    DataKinds
+    #-}
+
+import Goal.Core
+import qualified Goal.Core.Vector.Generic as G
+import qualified Goal.Core.Vector.Storable as S
+import qualified Goal.Core.Vector.Boxed as B
+
+import qualified Numeric.LinearAlgebra as H
+import qualified Criterion.Main as C
+import qualified System.Random.MWC.Probability as P
+
+
+
+--- Globals ---
+
+
+-- Sizes --
+
+type KRadius = 2
+type KDiameter = (2*KRadius + 1)
+type KNumber = 50
+type MSize = 50
+type MDepth = 50
+
+pkr :: Proxy KRadius
+pkr = Proxy
+
+pkd :: Proxy KDiameter
+pkd = Proxy
+
+pms :: Proxy MSize
+pms = Proxy
+
+kdmt,ms,krd :: Int
+krd = natValInt pkr
+kdmt = natValInt pkd
+ms = natValInt pms
+
+goalCorr
+    :: ( S.Matrix KNumber (KDiameter*KDiameter*MDepth) Double
+       , S.Matrix MDepth (MSize*MSize) Double )
+    -> S.Matrix KNumber (MSize*MSize) Double
+goalCorr (krns,mtx) = S.crossCorrelate2d pkr pkr pms pms krns mtx
+
+goalConv
+    :: ( S.Matrix KNumber (KDiameter*KDiameter*MDepth) Double
+       , S.Matrix KNumber (MSize*MSize) Double )
+    -> S.Matrix MDepth (MSize*MSize) Double
+goalConv (krns,mtx) = S.convolve2d pkr pkr pms pms krns mtx
+
+hmatrixCorr
+    :: (B.Vector KNumber (B.Vector MDepth (H.Matrix Double)), B.Vector MDepth (H.Matrix Double))
+    -> B.Vector KNumber (H.Matrix Double)
+hmatrixCorr (krnss,mtxs) = fromJust . B.fromList
+    $ [ foldr1 H.add [ H.corr2 krn mtx | (krn,mtx) <- B.toList $ B.zip krns mtxs ] | krns <- B.toList krnss ]
+
+repadHMatrix :: H.Matrix Double -> H.Matrix Double
+repadHMatrix hmtx =
+    let rw' = replicate krd . H.fromList $ replicate ms 0
+        hmtx' = H.fromRows . (rw' ++) . (++ rw') $ H.toRows hmtx
+        cl' = replicate krd . H.fromList $ replicate (ms + 2*krd) 0
+     in H.fromColumns . (cl' ++) . (++ cl') $ H.toColumns hmtx'
+
+--depadHMatrix :: H.Matrix Double -> H.Matrix Double
+--depadHMatrix hmtx =
+--    let hmtx' = H.fromRows . take mn . drop kdmt $ H.toRows hmtx
+--     in H.fromColumns . take mn . drop kdmt $ H.toColumns hmtx'
+
+
+-- Benchmark
+
+main :: IO ()
+main = do
+
+    let rnd :: P.Prob IO Double
+        rnd = P.uniformR (-1,1)
+
+    mtxv <- P.withSystemRandom . P.sample $ S.replicateM rnd
+    krnv <- P.withSystemRandom . P.sample $ S.replicateM rnd
+    mtxz <- P.withSystemRandom . P.sample $ S.replicateM rnd
+
+    let krn = G.Matrix krnv
+        mtx = G.Matrix mtxv
+        crr = goalCorr (krn,mtx)
+
+    let hmtxs = H.reshape ms . G.fromSized <$> G.convert (S.toRows mtx)
+        hkrnss = fmap (H.reshape kdmt . G.fromSized . G.convert) . B.breakEvery . G.convert <$> G.convert (S.toRows krn)
+        hcrr = hmatrixCorr (hkrnss,repadHMatrix <$> hmtxs)
+
+    putStrLn "Squared Error between HMatrix and Goal solutions:"
+    print . sum $
+        B.zipWith (\mtx1 mtx2 -> H.sumElements . H.cmap (^(2 :: Int)) . H.add mtx1 $ H.scale (-1) mtx2)
+         hcrr $ H.reshape ms . G.fromSized <$> G.convert (S.toRows crr)
+    putStrLn ""
+
+    putStrLn "Transpose Dot Product Difference:"
+    let mtx' = G.Matrix mtxz
+        cnv = goalConv (krn,mtx')
+    print $ S.dotProduct (G.toVector crr) (G.toVector mtx')
+        - S.dotProduct (G.toVector cnv) (G.toVector mtx)
+    putStrLn ""
+
+    C.defaultMain
+       [ C.bench "goal-corr" $ C.nf goalCorr (krn,mtx)
+       , C.bench "goal-conv" $ C.nf goalConv (krn,mtx')
+       , C.bench "hmatrix-corr" $ C.nf hmatrixCorr (hkrnss,hmtxs) ]
diff --git a/benchmarks/multiplications.hs b/benchmarks/multiplications.hs
new file mode 100644
--- /dev/null
+++ b/benchmarks/multiplications.hs
@@ -0,0 +1,129 @@
+{-# LANGUAGE DataKinds #-}
+
+import Goal.Core
+import qualified Goal.Core.Vector.Generic as G
+import qualified Goal.Core.Vector.Storable as S
+import qualified Goal.Core.Vector.Boxed as B
+
+import qualified Numeric.LinearAlgebra as H
+import qualified Criterion.Main as C
+import qualified System.Random.MWC.Probability as P
+
+
+--- Globals ---
+
+
+-- Sizes --
+
+type M = 1000
+type N = 10
+
+n,m :: Int
+m = 1000
+n = 10
+
+-- Matrices --
+
+goalMatrix1 :: S.Matrix M M Double
+goalMatrix1 = G.Matrix $ S.generate fromIntegral
+
+goalMatrix2 :: S.Matrix M N Double
+goalMatrix2 = G.Matrix $ S.generate fromIntegral
+
+bGoalMatrix1 :: B.Matrix M M Double
+bGoalMatrix1 = G.Matrix $ B.generate fromIntegral
+
+bGoalMatrix2 :: B.Matrix M N Double
+bGoalMatrix2 = G.Matrix $ B.generate fromIntegral
+
+goalVal :: (S.Matrix M M Double,S.Matrix M N Double) -> Double
+goalVal (m1,m2) =
+    let G.Matrix v = S.matrixMatrixMultiply m1 m2
+     in S.sum v
+
+goalVal2 :: (B.Matrix M M Double, B.Matrix M N Double) -> Double
+goalVal2 (m1,m2) =
+    let G.Matrix v = B.matrixMatrixMultiply m1 m2
+     in sum v
+
+hmatrixMatrix1 :: H.Matrix Double
+hmatrixMatrix1 = H.fromLists . take m . breakEvery m $ [0..]
+
+hmatrixMatrix2 :: H.Matrix Double
+hmatrixMatrix2 = H.fromLists . take m . breakEvery n $ [0..]
+
+hmatrixVal :: (H.Matrix Double,H.Matrix Double) -> Double
+hmatrixVal (m1,m2) =
+    let m3 = m1 H.<> m2
+     in H.sumElements m3
+
+-- Benchmark
+main :: IO ()
+main = do
+
+    let rnd :: P.Prob IO Double
+        rnd = P.uniformR (-1,1)
+
+    v1 <- P.withSystemRandom . P.sample $ S.replicateM rnd
+    v2 <- P.withSystemRandom . P.sample $ S.replicateM rnd
+
+    let m1 = G.Matrix v1
+        m2 = G.Matrix v2
+
+    let bm1 = G.Matrix $ G.convert v1
+        bm2 = G.Matrix $ G.convert v2
+
+    let m1'' = H.fromLists . take m . breakEvery m $!! S.toList v1
+        m2'' = H.fromLists . take m . breakEvery n $!! S.toList v2
+
+    C.defaultMain
+       [ C.bench "generative-goal" $ C.nf goalVal (goalMatrix1,goalMatrix2)
+       , C.bench "generative-goal2" $ C.nf goalVal2 (bGoalMatrix1,bGoalMatrix2)
+       , C.bench "generative-hmatrix" $ C.nf hmatrixVal (hmatrixMatrix1,hmatrixMatrix2)
+       , C.bench "random-goal" $ C.nf goalVal (m1,m2)
+       , C.bench "random-goal2" $ C.nf goalVal2 (bm1,bm2)
+       , C.bench "random-hmatrix" $ C.nf hmatrixVal (m1'',m2'') ]
+
+-- Sanity Check
+--sanityCheck :: IO ()
+--sanityCheck = do
+--
+--    let rnd :: P.Prob IO Double
+--        rnd = P.uniformR (-1,1)
+--
+--    putStrLn "Goal 1:"
+--    print $ S.matrixMatrixMultiply goalMatrix1 goalMatrix2
+--    putStrLn "Goal 2:"
+--    print $ B.matrixMatrixMultiply bGoalMatrix1 bGoalMatrix2
+--    putStrLn "Matrix:"
+--    print $ M.multStd2 matrixMatrix1 matrixMatrix2
+--    putStrLn "HMatrix:"
+--    print $ hmatrixMatrix1 H.<> hmatrixMatrix2
+--
+--    v1 <- P.withSystemRandom . P.sample $ S.replicateM rnd
+--    v2 <- P.withSystemRandom . P.sample $ S.replicateM rnd
+--
+--    let m1 = Matrix v1
+--        m2 = Matrix v2
+--
+--    let bm1 = Matrix $ G.convert v1
+--        bm2 = Matrix $ G.convert v2
+--
+--    let m1' = M.fromLists . take m . breakEvery m $!! S.toList v1
+--        m2' = M.fromLists . take m . breakEvery n $!! S.toList v2
+--
+--    let m1'' = H.fromLists . take m . breakEvery m $!! S.toList v1
+--        m2'' = H.fromLists . take m . breakEvery n $!! S.toList v2
+--
+--    putStrLn "Goal 1:"
+--    print $ goalVal (m1,m2)
+--    print $ S.matrixMatrixMultiply m1 m2
+--    putStrLn "Goal 2:"
+--    print $ goalVal2 (bm1,bm2)
+--    print $ B.matrixMatrixMultiply bm1 bm2
+--    putStrLn "Matrix:"
+--    print $ matrixVal (m1',m2')
+--    print $ M.multStd2 m1' m2'
+--    putStrLn "HMatrix:"
+--    print $ hmatrixVal (m1'',m2'')
+--    print $ m1'' H.<> m2''
diff --git a/benchmarks/outer-products.hs b/benchmarks/outer-products.hs
new file mode 100644
--- /dev/null
+++ b/benchmarks/outer-products.hs
@@ -0,0 +1,71 @@
+{-# LANGUAGE DataKinds,ScopedTypeVariables,FlexibleContexts #-}
+
+import Goal.Core
+import qualified Goal.Core.Vector.Generic as G
+import qualified Goal.Core.Vector.Storable as S
+import qualified Goal.Core.Vector.Boxed as B
+
+import qualified Numeric.LinearAlgebra as H
+import qualified Criterion.Main as C
+import qualified System.Random.MWC.Probability as P
+
+
+--- Globals ---
+
+
+type Rows = 1000
+type Columns = 1000
+
+n :: Int
+n = 20
+
+addMatrices :: (KnownNat n, KnownNat k) => S.Matrix n k Double -> S.Matrix n k Double -> S.Matrix n k Double
+addMatrices (G.Matrix mtx1) (G.Matrix mtx2) = G.Matrix $ S.add mtx1 mtx2
+
+
+--- Function ---
+
+
+averageOuterProduct2 v12s =
+    let ln = fromIntegral $ length v12s
+     in S.withMatrix (S.scale ln) $ foldr foldfun (G.Matrix $ S.replicate 0) v12s
+    where foldfun (v1,v2) mtx = addMatrices mtx $ S.outerProduct v1 v2
+
+averageWeightedOuterProduct2 wv12s =
+    foldr foldfun (G.Matrix $ S.replicate 0) wv12s
+        where foldfun (w,v1,v2) mtx = addMatrices mtx . S.outerProduct v1 $ S.scale w v2
+
+
+--- Main ---
+
+
+main :: IO ()
+main = do
+
+    let rnd :: P.Prob IO Double
+        rnd = P.uniformR (-1,1)
+
+    v1s :: [S.Vector Rows Double]
+        <- P.withSystemRandom . P.sample . replicateM n $ S.replicateM rnd
+    v2s :: [S.Vector Columns Double]
+        <- P.withSystemRandom . P.sample . replicateM n $ S.replicateM rnd
+
+    let ws :: [Double]
+        ws0 = [1..fromIntegral n]
+        ws = (/sum ws0) <$> ws0
+
+    let v12s = zip v1s v2s
+        wv12s = zip3 ws v1s v2s
+    let v11s = zip v1s v1s
+        wv11s = zip3 ws v1s v1s
+
+    C.defaultMain
+       [ C.bench "average-outer-product" $ C.nf averageOuterProduct2 v12s
+       , C.bench "bulk-outer-product" $ C.nf S.averageOuterProduct v12s
+       , C.bench "average-weighted-outer-product" $ C.nf averageWeightedOuterProduct2 wv12s
+       , C.bench "bulk-weighted-outer-product" $ C.nf S.weightedAverageOuterProduct wv12s
+       , C.bench "average-outer-product0" $ C.nf averageOuterProduct2 v11s
+       , C.bench "bulk-outer-product0" $ C.nf S.averageOuterProduct v11s
+       , C.bench "average-weighted-outer-product0" $ C.nf averageWeightedOuterProduct2 wv11s
+       , C.bench "bulk-weighted-outer-product0" $ C.nf S.weightedAverageOuterProduct wv11s
+       , C.bench "triangular-weighted-outer-product0" $ C.nf (S.lowerTriangular . S.weightedAverageOuterProduct) wv11s ]
diff --git a/goal-core.cabal b/goal-core.cabal
--- a/goal-core.cabal
+++ b/goal-core.cabal
@@ -1,43 +1,103 @@
+cabal-version: 3.0
 name: goal-core
-version: 0.1
-synopsis: Core imports for Geometric Optimization Libraries.
-description: Core provides a bunch of convenience functions, basic exports, as
-    well as plotting functionality, which are independent of the rest of the
-    library.
-license: BSD3
+synopsis: Common, non-geometric tools for use with Goal
+description: goal-core re-exports a number of other libraries, and provides a set of additional utility functions useful for scientific computing. In particular, implementations of Mealy Automata (Circuits), tools for working with CSV files and gnuplot, and a module which combines vector-sized vectors with hmatrix.
+license: BSD-3-Clause
 license-file: LICENSE
+extra-source-files: README.md
+version: 0.20
 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.Core,
-        Goal.Core.Plot,
-        Goal.Core.Plot.Contour
+        Goal.Core.Util,
+        Goal.Core.Project,
+        Goal.Core.Circuit,
+        Goal.Core.Vector.Storable,
+        Goal.Core.Vector.Generic,
+        Goal.Core.Vector.Generic.Internal,
+        Goal.Core.Vector.Generic.Mutable,
+        Goal.Core.Vector.Boxed
     build-depends:
-        base==4.*,
-        data-default-class==0.0.*,
-        lens==4.*,
-        containers==0.5.*,
-        colour==2.3.*,
-        gtk==0.14.*,
-        cairo==0.13.*,
-        Chart==1.5.*,
-        Chart-cairo==1.5.*,
-        Chart-gtk==1.5.*
-    default-extensions: TemplateHaskell
+        base >= 4.13 && < 4.15,
+        directory,
+        containers,
+        vector,
+        math-functions,
+        hmatrix,
+        vector-sized,
+        finite-typelits,
+        ghc-typelits-knownnat,
+        ghc-typelits-natnormalise,
+        deepseq,
+        process,
+        hmatrix-gsl,
+        primitive,
+        bytestring,
+        cassava,
+        async,
+        criterion,
+        optparse-applicative
     default-language: Haskell2010
-    ghc-options: -O2 -Wall -fno-warn-type-defaults -fno-warn-missing-signatures
+    default-extensions:
+        ScopedTypeVariables,
+        ExplicitNamespaces,
+        TypeOperators,
+        KindSignatures,
+        DataKinds,
+        RankNTypes,
+        TypeFamilies,
+        FlexibleContexts,
+        MultiParamTypeClasses,
+        ConstraintKinds,
+        GADTs,
+        NoStarIsType,
+        FlexibleInstances
+    ghc-options: -Wall -O2
 
-executable contours
-    main-is: contours.hs
-    hs-source-dirs: scripts
-    ghc-options: -Wall -O2 -threaded -rtsopts -fno-warn-type-defaults
-        -fno-warn-missing-signatures -fno-warn-unused-do-bind
+benchmark outer-products
+    type: exitcode-stdio-1.0
+    main-is: outer-products.hs
+    hs-source-dirs: benchmarks
     build-depends:
-        base==4.*,
-        goal-core==0.1
+        base,
+        hmatrix,
+        mwc-random,
+        mwc-probability,
+        criterion,
+        goal-core
     default-language: Haskell2010
+    ghc-options: -Wall -O2
+
+benchmark convolutions
+    type: exitcode-stdio-1.0
+    main-is: convolutions.hs
+    hs-source-dirs: benchmarks
+    build-depends:
+        base,
+        hmatrix,
+        mwc-random,
+        mwc-probability,
+        criterion,
+        goal-core
+    default-language: Haskell2010
+    ghc-options: -Wall -O2
+
+benchmark multiplications
+    type: exitcode-stdio-1.0
+    main-is: multiplications.hs
+    hs-source-dirs: benchmarks
+    build-depends:
+        base,
+        hmatrix,
+        mwc-random,
+        mwc-probability,
+        criterion,
+        goal-core
+    default-language: Haskell2010
+    ghc-options: -Wall -O2
diff --git a/scripts/contours.hs b/scripts/contours.hs
deleted file mode 100644
--- a/scripts/contours.hs
+++ /dev/null
@@ -1,27 +0,0 @@
-import Goal.Core
-
---- Test ---
-
-main :: IO ()
-main = do
-    let f x y = 0.2 * x + 0.2 * y + sin x + sin y
-        sz = 7.5
-        nstp = 1000
-        rng = (-sz,sz,nstp)
-        niso = 20
-        cntrs = contours rng rng niso f
-        clrs = rgbaGradient (0,0,1,1) (1,0,0,1) niso
-        lyt = toRenderable . execEC $ do
-
-            layout_title .= "Contours of a 2D Sin Curve"
-
-            sequence $ do
-
-                ((_,cntr),clr) <- zip cntrs clrs
-
-                return . plot . liftEC $ do
-
-                    plot_lines_style .= solidLine 3 clr
-                    plot_lines_values .= cntr
-
-    renderableToAspectWindow False 800 600 lyt