diff --git a/benchmark/Main.lhs b/benchmark/Main.lhs
--- a/benchmark/Main.lhs
+++ b/benchmark/Main.lhs
@@ -12,15 +12,15 @@
 
 We draw 2e3 normally distributed 2D points:
 
-> points :: [Point Int]
-> points = [Point (vector n) n | n <- [1..2000-1]]
->   where vector n = fromList [normals!!n, normals!!n*2]
->         normals  = mkNormals 0x29a
+> points :: [Point]
+> points = concat $ [0..2000-1] `forM` \n -> do
+>   return . Point . fromList $ [normals!!n, normals!!(n*2)]
+>   where normals = take 4000 $ mkNormals 0x29a
 
 Three of which form the centroids of our initial clusters:
  
 > clusters :: [Cluster]
-> clusters = take 3 . map (uncurry Cluster) $ zip [0..] (map point points)
+> clusters = map (uncurry Cluster) $ zip [0..] (take 3 points)
 
 To correctly benchmark the result of a pure function, we need to be able to
 evaluate it to normal form:
@@ -28,8 +28,8 @@
 > instance NFData Cluster where
 >   rnf (Cluster i c) = rnf i `seq` rnf c
 >
-> instance NFData a => NFData (Point a) where
->   rnf (Point a b) = rnf a `seq` rnf b
+> instance NFData Point where
+>   rnf (Point v) = rnf v
  
 Together the subject of our benchmarks:
  
diff --git a/kmeans-par.cabal b/kmeans-par.cabal
--- a/kmeans-par.cabal
+++ b/kmeans-par.cabal
@@ -1,5 +1,5 @@
 name:                kmeans-par
-version:             1.2.0
+version:             1.3.0
 synopsis:            Sequential and parallel implementations of Lloyd's algorithm.
 license:             MIT
 license-file:        LICENSE
diff --git a/src/Algorithms/Lloyd/Sequential.lhs b/src/Algorithms/Lloyd/Sequential.lhs
--- a/src/Algorithms/Lloyd/Sequential.lhs
+++ b/src/Algorithms/Lloyd/Sequential.lhs
@@ -5,29 +5,30 @@
 
 > module Algorithms.Lloyd.Sequential where
 > 
-> import Prelude hiding (zipWith, map, foldr1, replicate)
+> import Prelude hiding (zipWith, map, foldr, replicate)
+> import Control.Applicative ((<$>))
 > import Control.Monad (guard, forM_)
+> import Data.Foldable (Foldable(foldr))
 > import Data.Function (on)
 > import Data.Functor.Extras ((..:),(...:))
 > import Data.List (minimumBy)
 > import Data.Metric (Metric(..))
 > import Data.Semigroup (Semigroup(..))
-> import Data.Vector (Vector(..), toList, create, zipWith, map, foldr1, empty, replicate)
+> import Data.Vector (Vector(..), toList, create, zipWith, map, empty, replicate, cons)
+> import qualified Data.Vector as V (length)
 > import qualified Data.Vector.Mutable as MV (replicate, read, write)
  
-A data point is comprised of an n-dimensional real vector, and an
-associated value:
+A data point is represented by an n-dimensional real vector:
 
-> data Point a = Point 
+> data Point = Point 
 >   { point :: Vector Double
->   , value :: a
->   } 
+>   } deriving (Eq, Show)
 
-Two data-points are defined as equivalent if they've a point vector in
-common:
+Point isn't parametrised, so we can't define a functor (or any functor-family)
+instances. Instead, we define a specialised map, `pmap`:
 
-> instance Eq (Point a) where
->   (==) = (==) `on` point
+> pmap :: (Double -> Double) -> Point -> Point
+> pmap f = Point . map f . point
 
 To define distance between data-points, we admit arbitrary metric spaces
 over real vectors; **However:** it should be noted that convergence is
@@ -35,14 +36,19 @@
 same criterion as minimum-distance assignment (`assign` below.) For this
 reason, Euclidean-distance is popularly used for mean assignment.
 
-> useMetric :: Metric a => (Vector Double -> a) -> Vector Double -> Vector Double -> Double
-> useMetric = on distance
+> useMetric :: Metric a => (Vector Double -> a) -> Point -> Point -> Double
+> useMetric metric = distance `on` (metric . point)
+ 
+`Point`s may be combined using simple vector addition:
 
+> instance Semigroup Point where
+>   (<>) = Point ..: zipWith (+) `on` point
+
 Clusters are represented by an identifier and a centroid:
 
 > data Cluster = Cluster
 >   { identifier :: Int
->   , centroid   :: Vector Double
+>   , centroid   :: Point
 >   } deriving (Show)
 
 Clusters are treated as equivalent if they share a centroid:
@@ -54,55 +60,54 @@
 of points in a set when constructing new clusters. It consists of a count of
 included points, as well as the vector sum of the constituent points:
 
-> data PointSum = PointSum Int (Vector Double) deriving (Show)
+> data PointSum = PointSum Int Point deriving (Show)
 
 Two `PointSum`s can be combined by summing corresponding values:
 
 > instance Semigroup PointSum where
->   PointSum c0 p0 <> PointSum c1 p1 = PointSum (c0+c1) (zipWith (+) p0 p1)
+>   PointSum c0 p0 <> PointSum c1 p1 = PointSum (c0+c1) (p0<>p1)
 
 Without dependent types, we can't define a valid Monoid instance for PointSum,
 as the identity varies with the dimensions of the point vector. But we can
 construct one at runtime:
 
 > emptyPointSum :: Int -> PointSum
-> emptyPointSum length = PointSum 0 $ replicate length 0
-
-A point sum is constructed by incrementally adding points likeso:
-
-> addPoint :: PointSum -> Point a -> PointSum
-> addPoint sum = (sum <>) . PointSum 1 . point
+> emptyPointSum length = PointSum 0 . Point $ replicate length 0
 
 A point sum can be turned into a cluster by computing its centroid, the average
 of all points associated with the cluster:
 
 > toCluster :: Int -> PointSum -> Cluster
-> toCluster cid (PointSum count vector) = Cluster
+> toCluster cid (PointSum count point) = Cluster
 >   { identifier = cid
->   , centroid   = map (//count) vector
+>   , centroid   = pmap (// count) point
 >   }
 >
 > (//) :: Double -> Int -> Double
 > x // y = x / fromIntegral y
-           
+         
 After each iteration, points are associated with the cluster having the nearest
 centroid:
 
-> closestCluster :: Metric a => (Vector Double -> a) -> [Cluster] -> Point b -> Cluster
-> closestCluster (useMetric -> d) clusters (point->point) = fst . minimumBy (compare `on` snd) $ do
+> closestCluster :: Metric a => (Vector Double -> a) -> [Cluster] -> Point -> Cluster
+> closestCluster (useMetric -> d) clusters point = fst . minimumBy (compare `on` snd) $ do
 >   cluster <- clusters
 >   return (cluster, point `d` centroid cluster)
 >
-> assign :: Metric a => (Vector Double -> a) -> [Cluster] -> [Point b] -> Vector PointSum
+> assign :: Metric a => (Vector Double -> a) -> [Cluster] -> [Point] -> Vector (Vector Point)
 > assign metric clusters points = let nc = length clusters in create $ do
->   vector <- MV.replicate nc $ emptyPointSum nc
+>   vector <- MV.replicate nc empty
 >   points `forM_` \point -> do
->     let cluster  = closestCluster metric clusters point 
->         position = identifier cluster
->     sum <- MV.read vector position
->     MV.write vector position $! addPoint sum point
+>     let position = identifier $ closestCluster metric clusters point
+>     points' <- MV.read vector position
+>     MV.write vector position $! point `cons` points'
 >   return vector
 >
+> assignPS :: Metric a => (Vector Double -> a) -> [Cluster] -> [Point] -> Vector PointSum
+> assignPS metric clusters points = reduce <$> assign metric clusters points
+>   where reduce  = foldr (<>) (emptyPointSum length') . map (PointSum 1)
+>         length' = V.length . point $ head points
+>
 > makeNewClusters :: Vector PointSum -> [Cluster]
 > makeNewClusters vector = do
 >   (pointSum@(PointSum count _), index) <- zip (toList vector) [0..]
@@ -111,23 +116,28 @@
 >                     -- computing the centroid.
 >   return $ toCluster index pointSum
 >
-> step :: Metric a => (Vector Double -> a) -> [Cluster] -> [Point b] -> [Cluster]
-> step = makeNewClusters ...: assign
+> step :: Metric a => (Vector Double -> a) -> [Cluster] -> [Point] -> [Cluster]
+> step = makeNewClusters ...: assignPS
+>
 
 The algorithm consists of iteratively finding the centroid of each existing
 cluster, then reallocating points according to which centroid is closest, until
 convergence. As the algorithm isn't guaranteed to converge, we cut execution if
 convergence hasn't been observed after eighty iterations:
 
-> kmeans :: Metric a => Int -> (Vector Double -> a) -> [Point b] -> [Cluster] -> [Cluster]
-> kmeans expectDivergent metric = kmeans' expectDivergent metric 0 
+> computeClusters :: Metric a => Int -> (Vector Double -> a) -> [Point] -> [Cluster] -> [Cluster]
+> computeClusters expectDivergent metric = computeClusters' expectDivergent metric 0 
 >
-> kmeans' :: Metric a => Int -> (Vector Double -> a) -> Int -> [Point b] -> [Cluster] -> [Cluster]
-> kmeans' expectDivergent metric iterations points clusters 
+> computeClusters' :: Metric a => Int -> (Vector Double -> a) -> Int -> [Point] -> [Cluster] -> [Cluster]
+> computeClusters' expectDivergent metric iterations points clusters 
 >   | iterations >= expectDivergent = clusters
 >   | clusters' == clusters         = clusters 
->   | otherwise                     = kmeans' expectDivergent metric (succ iterations) points clusters'
+>   | otherwise                     = computeClusters' expectDivergent metric (succ iterations) points clusters'
 >   where clusters' = step metric clusters points
+>
+> kmeans :: Metric a => Int -> (Vector Double -> a) -> [Point] -> [Cluster] -> Vector (Vector Point)
+> kmeans expectDivergent metric points initial = assign metric clusters points
+>   where clusters = computeClusters 80 metric points initial
 
 A note on initialisation: typically clusters are randomly assigned to data
 points prior to the first update step.
diff --git a/src/Algorithms/Lloyd/Strategies.lhs b/src/Algorithms/Lloyd/Strategies.lhs
--- a/src/Algorithms/Lloyd/Strategies.lhs
+++ b/src/Algorithms/Lloyd/Strategies.lhs
@@ -12,7 +12,7 @@
 > import Data.Metric (Metric(..))
 > import Data.Semigroup (Semigroup(..))
 > import Data.Vector (Vector(..), zipWith)
-> import Algorithms.Lloyd.Sequential (Cluster(..), Point(..), PointSum(..), makeNewClusters, assign)
+> import Algorithms.Lloyd.Sequential (Cluster(..), Point(..), PointSum(..), makeNewClusters, assignPS, assign)
 
 We can combine two vectors of some same type $t$ provided we know how to
 combine two $t$s:
@@ -23,8 +23,8 @@
 Step is modified to, given a partitioned list of points, perform
 classification in parallel:
 
-> step :: Metric a => (Vector Double -> a) -> [Cluster] -> [[Point b]] -> [Cluster]
-> step = makeNewClusters . foldr1 (<>) . with (parTraversable rseq) ..: map ..: assign 
+> step :: Metric a => (Vector Double -> a) -> [Cluster] -> [[Point]] -> [Cluster]
+> step = makeNewClusters . foldr1 (<>) . with (parTraversable rseq) ..: map ..: assignPS
 >
 > with :: Strategy a -> a -> a
 > with = flip using
@@ -38,12 +38,16 @@
 are too few items, and those items vary in cost, some of our cores may
 be unused for part of the computation.
 
-> kmeans :: Metric a => Int -> (Vector Double -> a) -> Int -> [Point b] -> [Cluster] -> [Cluster]
-> kmeans expectDivergent metric = kmeans' expectDivergent metric 0  ..: chunksOf
+> computeClusters :: Metric a => Int -> (Vector Double -> a) -> Int -> [Point] -> [Cluster] -> [Cluster]
+> computeClusters expectDivergent metric = computeClusters' expectDivergent metric 0  ..: chunksOf
 >
-> kmeans' :: Metric a => Int -> (Vector Double -> a) -> Int -> [[Point b]] -> [Cluster] -> [Cluster]
-> kmeans' expectDivergent metric iterations points clusters 
+> computeClusters' :: Metric a => Int -> (Vector Double -> a) -> Int -> [[Point]] -> [Cluster] -> [Cluster]
+> computeClusters' expectDivergent metric iterations points clusters 
 >   | iterations >= expectDivergent = clusters
 >   | clusters' == clusters         = clusters 
->   | otherwise                     = kmeans' expectDivergent metric (succ iterations) points clusters'
+>   | otherwise                     = computeClusters' expectDivergent metric (succ iterations) points clusters'
 >   where clusters' = step metric clusters points
+>
+> kmeans :: Metric a => Int -> (Vector Double -> a) -> Int -> [Point] -> [Cluster] -> Vector (Vector Point)
+> kmeans expectDivergent metric iterations points initial = assign metric clusters points
+>   where clusters = computeClusters 80 metric iterations points initial