diff --git a/src/Data/Algorithm/TSNE.hs b/src/Data/Algorithm/TSNE.hs
--- a/src/Data/Algorithm/TSNE.hs
+++ b/src/Data/Algorithm/TSNE.hs
@@ -11,8 +11,11 @@
 import Pipes
 
 import Data.Algorithm.TSNE.Types
-import Data.Algorithm.TSNE.Internals
 import Data.Algorithm.TSNE.Utils
+import Data.Algorithm.TSNE.Preparation
+import Data.Algorithm.TSNE.Run3D
+import Data.Algorithm.TSNE.Run2D
+
 
 -- | Generates an infinite stream of 3D tSNE iterations.
 tsne3D :: TSNEOptions -> TSNEInput -> Producer TSNEOutput3D IO ()
diff --git a/src/Data/Algorithm/TSNE/Checks.hs b/src/Data/Algorithm/TSNE/Checks.hs
--- a/src/Data/Algorithm/TSNE/Checks.hs
+++ b/src/Data/Algorithm/TSNE/Checks.hs
@@ -17,3 +17,24 @@
 isRectangular :: [[a]] -> Bool
 isRectangular xss = has2DShape (shape2D xss) xss
 
+inputSize :: TSNEInput -> Int
+inputSize = length
+
+inputValueSize :: TSNEInput -> Int
+inputValueSize i = w 
+    where (w,h) = shape2D i 
+
+inputIsValid :: TSNEInput -> Either String ()
+inputIsValid [] = Left "empty input data"
+inputIsValid xss
+    | not (isRectangular xss) = Left "input data values are not all the same length"
+    | otherwise = Right () 
+
+isValidStateForInput :: Int -> TSNEInput -> TSNEState -> Either String ()
+isValidStateForInput d i st
+    | not (has2DShape (n,d) s) = Left $ "solution is wrong shape: " ++ show (shape2D s) 
+    | otherwise = Right ()
+        where
+            n = inputSize i
+            s = stSolution st  
+
diff --git a/src/Data/Algorithm/TSNE/Internals.hs b/src/Data/Algorithm/TSNE/Internals.hs
deleted file mode 100644
--- a/src/Data/Algorithm/TSNE/Internals.hs
+++ /dev/null
@@ -1,207 +0,0 @@
-module Data.Algorithm.TSNE.Internals where
-
-import Control.Applicative
-import Control.DeepSeq
-import Control.Exception (assert)
-import Data.Default (def)
-import Data.List(zipWith4)
-import Data.Random.Normal (normalsIO')
-import Pipes
---import Debug.Trace
-
-import Data.Algorithm.TSNE.Types
-import Data.Algorithm.TSNE.Utils
-import Data.Algorithm.TSNE.Checks
-
-
-inputSize :: TSNEInput -> Int
-inputSize = length
-
-inputValueSize :: TSNEInput -> Int
-inputValueSize i = w 
-    where (w,h) = shape2D i 
-
-inputIsValid :: TSNEInput -> Either String ()
-inputIsValid [] = Left "empty input data"
-inputIsValid xss
-    | not (isRectangular xss) = Left "input data values are not all the same length"
-    | otherwise = Right () 
-
-isValidStateForInput :: TSNEInput -> TSNEState -> Either String ()
-isValidStateForInput i st
-    | not (has2DShape (n,3) s) = Left $ "solution is wrong shape: " ++ show (shape2D s) 
-    | otherwise = Right ()
-        where
-            n = inputSize i
-            s = stSolution st    
-
-initState3D :: Int -> IO TSNEState
-initState3D n = do
-    s <- initSolution3D n
-    return $ TSNEState 0 s (rr 1) (rr 0)
-        where
-            rr = repeat.repeat
-
-initState2D :: Int -> IO TSNEState
-initState2D n = do
-    s <- initSolution2D n
-    return $ TSNEState 0 s (rr 1) (rr 0)
-        where
-            rr = repeat.repeat
-
-initSolution3D :: Int -> IO [[Double]]
-initSolution3D n = do
-    let ns = normalsIO' (0, 1e-4)
-    xs <- ns
-    ys <- ns
-    zs <- ns
-    return $ take n <$> [xs,ys,zs]
-
-initSolution2D :: Int -> IO [[Double]]
-initSolution2D n = do
-    let ns = normalsIO' (0, 1e-4)
-    xs <- ns
-    ys <- ns
-    return $ take n <$> [xs,ys]
-
-runTSNE3D :: TSNEOptions -> TSNEInput -> [[Probability]] -> TSNEState -> Producer TSNEOutput3D IO ()
-runTSNE3D opts vs ps st = do
-    yield $ output3D ps st
-    let st' = force $ stepTSNE opts vs ps st
-    runTSNE3D opts vs ps st'
-
-runTSNE2D :: TSNEOptions -> TSNEInput -> [[Probability]] -> TSNEState -> Producer TSNEOutput2D IO ()
-runTSNE2D opts vs ps st = do
-    yield $ output2D ps st
-    let st' = force $ stepTSNE opts vs ps st
-    runTSNE2D opts vs ps st'
-
-stepTSNE :: TSNEOptions -> TSNEInput -> [[Probability]] -> TSNEState -> TSNEState
-stepTSNE opts vs ps st = TSNEState i' s'' g' d'
-    where
-        i = stIteration st
-        s = stSolution st
-        g = stGains st
-        d = stDeltas st
-        gr = gradients ps st
-        i' = i + 1
-        s' = recenter $ z (+) s d'
-        g' = z3 newGain g d gr
-        d' = z3 (newDelta (tsneLearningRate opts) i) g' d gr
-        z = zipWith.zipWith
-        z3 = zipWith3.zipWith3
-        s'' = assert (length s' == length vs) s'
-
-newGain :: Gain -> Delta -> Gradient -> Gain
-newGain g d gr = max 0.01 g'
-    where
-        g' = if signum d == signum gr 
-                then g * 0.8
-                else g + 0.2  
-
-newDelta :: Double -> Int -> Gain -> Delta -> Gradient -> Delta
-newDelta e i g' d gr = (m * d) - (e * g' * gr)
-    where
-        m = if i < 250 then 0.5 else 0.8
-
-gradients :: [[Probability]] -> TSNEState -> [[Gradient]]
-gradients pss st = gradient <$> ss
-    where
-        gradient :: [Double] -> [Gradient]
-        gradient s = zipWith4 (f s) s pss qss qss'
-        ss = stSolution st
-        i = stIteration st
-        qss = qdist ss
-        qss' = qdist' ss 
-        f :: [Double] -> Double -> [Double] -> [Double] -> [Double] -> Gradient
-        f s x ps qs qs' = sum $ zipWith4 g s ps qs qs'
-            where
-                g y p q q' = m * (x - y)
-                    where
-                        m = 4 * (k * p - q') * q
-                        k = if i < 100 then 4 else 1
-
-solution3D :: [[Double]] -> [Position3D]
-solution3D (xs:ys:zs:_) = zip3 xs ys zs
-
-output3D :: [[Double]] -> TSNEState -> TSNEOutput3D
-output3D pss st = TSNEOutput3D i s c
-    where
-        i = stIteration st
-        s = (solution3D . stSolution) st
-        c = cost pss st
-
-solution2D :: [[Double]] -> [Position2D]
-solution2D (xs:ys:_) = zip xs ys
-
-output2D :: [[Double]] -> TSNEState -> TSNEOutput2D
-output2D pss st = TSNEOutput2D i s c
-    where
-        i = stIteration st
-        s = (solution2D . stSolution) st
-        c = cost pss st
-
-cost :: [[Double]] -> TSNEState -> Double
-cost pss st = sumsum $ (zipWith.zipWith) c pss (qdist' (stSolution st))
-    where
-        c p q = -p * log q 
-
-targetEntropy :: TSNEOptions -> Entropy
-targetEntropy = log.realToFrac.tsnePerplexity
-
-data Beta = Beta {
-    betaValue :: Double,
-    betaMin :: Double,
-    betaMax :: Double
-}
-
-neighbourProbabilities :: TSNEOptions -> TSNEInput -> [[Probability]]
-neighbourProbabilities opts vs = symmetrize $ rawNeighbourProbabilities opts vs
-
-rawNeighbourProbabilities :: TSNEOptions -> TSNEInput -> [[Probability]]
-rawNeighbourProbabilities opts vs = map np vs
-    where 
-        np a = aps (beta a) vs a
-        beta a = betaValue $ binarySearchBeta opts vs a
-
-        aps :: Double -> TSNEInput -> TSNEInputValue -> [Probability]
-        aps beta bs a = map pj' bs
-            where
-                psum = sum $ map pj bs
-                pj b 
-                    | a == b    = 0
-                    | otherwise = exp $ -(distanceSquared a b) * beta 
-                pj' b = pj b / psum
-
-binarySearchBeta :: TSNEOptions -> TSNEInput -> TSNEInputValue -> Beta
-binarySearchBeta opts vs = binarySearchBeta' opts vs 1e-4 0 (Beta 1 (-infinity) infinity)
-
-binarySearchBeta' :: TSNEOptions -> TSNEInput -> Double -> Int -> Beta -> TSNEInputValue -> Beta
-binarySearchBeta' opts bs tol i beta a
-    | i == 50            = beta
-    | abs (e - t) < tol  = beta
-    | e > t              = r $ incPrecision beta
-    | otherwise          = r $ decPrecision beta 
-        where
-            t = targetEntropy opts
-            e = entropyForInputValue (betaValue beta) bs a
-            incPrecision (Beta b _ bmax) 
-                | bmax == infinity = Beta (b * 2) b bmax
-                | otherwise        = Beta ((b + bmax) / 2) b bmax
-            decPrecision (Beta b bmin _) 
-                | bmin == -infinity = Beta (b / 2) bmin b
-                | otherwise         = Beta ((b + bmin) / 2) bmin b
-            r beta' = binarySearchBeta' opts bs tol (i+1) beta' a 
-
-entropyForInputValue :: Double -> TSNEInput -> TSNEInputValue -> Entropy
-entropyForInputValue beta bs a = sum $ map h bs
-    where
-        h b = if x > 1e-7 then -x * log x else 0
-            where x = pj' b
-        psum = sum $ map pj bs
-        pj b 
-            | a == b    = 0
-            | otherwise = exp $ -(distanceSquared a b) * beta 
-        pj' b = pj b / psum
-
-
diff --git a/src/Data/Algorithm/TSNE/Preparation.hs b/src/Data/Algorithm/TSNE/Preparation.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Algorithm/TSNE/Preparation.hs
@@ -0,0 +1,65 @@
+module Data.Algorithm.TSNE.Preparation where
+
+import Data.Algorithm.TSNE.Types
+import Data.Algorithm.TSNE.Utils
+
+
+targetEntropy :: TSNEOptions -> Entropy
+targetEntropy = log.realToFrac.tsnePerplexity
+
+data Beta = Beta {
+    betaValue :: Double,
+    betaMin :: Double,
+    betaMax :: Double
+}
+
+neighbourProbabilities :: TSNEOptions -> TSNEInput -> [[Probability]]
+neighbourProbabilities opts vs = symmetrize $ rawNeighbourProbabilities opts vs
+
+rawNeighbourProbabilities :: TSNEOptions -> TSNEInput -> [[Probability]]
+rawNeighbourProbabilities opts vs = map np vs
+    where 
+        np a = aps (beta a) vs a
+        beta a = betaValue $ binarySearchBeta opts vs a
+
+        aps :: Double -> TSNEInput -> TSNEInputValue -> [Probability]
+        aps beta bs a = map pj' bs
+            where
+                psum = sum $ map pj bs
+                pj b 
+                    | a == b    = 0
+                    | otherwise = exp $ -(distanceSquared a b) * beta 
+                pj' b = pj b / psum
+
+binarySearchBeta :: TSNEOptions -> TSNEInput -> TSNEInputValue -> Beta
+binarySearchBeta opts vs = binarySearchBeta' opts vs 1e-4 0 (Beta 1 (-infinity) infinity)
+
+binarySearchBeta' :: TSNEOptions -> TSNEInput -> Double -> Int -> Beta -> TSNEInputValue -> Beta
+binarySearchBeta' opts bs tol i beta a
+    | i == 50            = beta
+    | abs (e - t) < tol  = beta
+    | e > t              = r $ incPrecision beta
+    | otherwise          = r $ decPrecision beta 
+        where
+            t = targetEntropy opts
+            e = entropyForInputValue (betaValue beta) bs a
+            incPrecision (Beta b _ bmax) 
+                | bmax == infinity = Beta (b * 2) b bmax
+                | otherwise        = Beta ((b + bmax) / 2) b bmax
+            decPrecision (Beta b bmin _) 
+                | bmin == -infinity = Beta (b / 2) bmin b
+                | otherwise         = Beta ((b + bmin) / 2) bmin b
+            r beta' = binarySearchBeta' opts bs tol (i+1) beta' a 
+
+entropyForInputValue :: Double -> TSNEInput -> TSNEInputValue -> Entropy
+entropyForInputValue beta bs a = sum $ map h bs
+    where
+        h b = if x > 1e-7 then -x * log x else 0
+            where x = pj' b
+        psum = sum $ map pj bs
+        pj b 
+            | a == b    = 0
+            | otherwise = exp $ -(distanceSquared a b) * beta 
+        pj' b = pj b / psum
+
+
diff --git a/src/Data/Algorithm/TSNE/Run2D.hs b/src/Data/Algorithm/TSNE/Run2D.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Algorithm/TSNE/Run2D.hs
@@ -0,0 +1,42 @@
+module Data.Algorithm.TSNE.Run2D where
+
+import Control.Applicative
+import Control.DeepSeq
+import Data.Random.Normal (normalsIO')
+import Pipes
+
+import Data.Algorithm.TSNE.Types
+import Data.Algorithm.TSNE.Preparation
+import Data.Algorithm.TSNE.Stepping
+
+
+initState2D :: Int -> IO TSNEState
+initState2D n = do
+    s <- initSolution2D n
+    return $ TSNEState 0 s (rr 1) (rr 0)
+        where
+            rr = repeat.repeat
+
+initSolution2D :: Int -> IO [[Double]]
+initSolution2D n = do
+    let ns = normalsIO' (0, 1e-4)
+    xs <- ns
+    ys <- ns
+    return $ take n <$> [xs,ys]
+
+runTSNE2D :: TSNEOptions -> TSNEInput -> [[Probability]] -> TSNEState -> Producer TSNEOutput2D IO ()
+runTSNE2D opts vs ps st = do
+    yield $ output2D ps st
+    let st' = force $ stepTSNE opts vs ps st
+    runTSNE2D opts vs ps st'
+
+solution2D :: [[Double]] -> [Position2D]
+solution2D (xs:ys:_) = zip xs ys
+
+output2D :: [[Double]] -> TSNEState -> TSNEOutput2D
+output2D pss st = TSNEOutput2D i s c
+    where
+        i = stIteration st
+        s = (solution2D . stSolution) st
+        c = cost pss st
+
diff --git a/src/Data/Algorithm/TSNE/Run3D.hs b/src/Data/Algorithm/TSNE/Run3D.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Algorithm/TSNE/Run3D.hs
@@ -0,0 +1,43 @@
+module Data.Algorithm.TSNE.Run3D where
+
+import Control.Applicative
+import Control.DeepSeq
+import Data.Random.Normal (normalsIO')
+import Pipes
+
+import Data.Algorithm.TSNE.Types
+import Data.Algorithm.TSNE.Preparation
+import Data.Algorithm.TSNE.Stepping
+
+
+initState3D :: Int -> IO TSNEState
+initState3D n = do
+    s <- initSolution3D n
+    return $ TSNEState 0 s (rr 1) (rr 0)
+        where
+            rr = repeat.repeat
+
+initSolution3D :: Int -> IO [[Double]]
+initSolution3D n = do
+    let ns = normalsIO' (0, 1e-4)
+    xs <- ns
+    ys <- ns
+    zs <- ns
+    return $ take n <$> [xs,ys,zs]
+
+runTSNE3D :: TSNEOptions -> TSNEInput -> [[Probability]] -> TSNEState -> Producer TSNEOutput3D IO ()
+runTSNE3D opts vs ps st = do
+    yield $ output3D ps st
+    let st' = force $ stepTSNE opts vs ps st
+    runTSNE3D opts vs ps st'
+
+solution3D :: [[Double]] -> [Position3D]
+solution3D (xs:ys:zs:_) = zip3 xs ys zs
+
+output3D :: [[Double]] -> TSNEState -> TSNEOutput3D
+output3D pss st = TSNEOutput3D i s c
+    where
+        i = stIteration st
+        s = (solution3D . stSolution) st
+        c = cost pss st
+
diff --git a/src/Data/Algorithm/TSNE/Stepping.hs b/src/Data/Algorithm/TSNE/Stepping.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Algorithm/TSNE/Stepping.hs
@@ -0,0 +1,60 @@
+module Data.Algorithm.TSNE.Stepping where
+
+import Control.Applicative
+import Control.Exception (assert)
+import Data.List(zipWith4)
+
+import Data.Algorithm.TSNE.Types
+import Data.Algorithm.TSNE.Utils
+
+
+stepTSNE :: TSNEOptions -> TSNEInput -> [[Probability]] -> TSNEState -> TSNEState
+stepTSNE opts vs ps st = TSNEState i' s'' g' d'
+    where
+        i = stIteration st
+        s = stSolution st
+        g = stGains st
+        d = stDeltas st
+        gr = gradients ps st
+        i' = i + 1
+        s' = recenter $ z (+) s d'
+        g' = z3 newGain g d gr
+        d' = z3 (newDelta (tsneLearningRate opts) i) g' d gr
+        z = zipWith.zipWith
+        z3 = zipWith3.zipWith3
+        s'' = assert (length s' == length vs) s'
+
+newGain :: Gain -> Delta -> Gradient -> Gain
+newGain g d gr = max 0.01 g'
+    where
+        g' = if signum d == signum gr 
+                then g * 0.8
+                else g + 0.2  
+
+newDelta :: Double -> Int -> Gain -> Delta -> Gradient -> Delta
+newDelta e i g' d gr = (m * d) - (e * g' * gr)
+    where
+        m = if i < 250 then 0.5 else 0.8
+
+gradients :: [[Probability]] -> TSNEState -> [[Gradient]]
+gradients pss st = gradient <$> ss
+    where
+        gradient :: [Double] -> [Gradient]
+        gradient s = zipWith4 (f s) s pss qss qss'
+        ss = stSolution st
+        i = stIteration st
+        qss = qdist ss
+        qss' = qdist' ss 
+        f :: [Double] -> Double -> [Double] -> [Double] -> [Double] -> Gradient
+        f s x ps qs qs' = sum $ zipWith4 g s ps qs qs'
+            where
+                g y p q q' = m * (x - y)
+                    where
+                        m = 4 * (k * p - q') * q
+                        k = if i < 100 then 4 else 1
+
+cost :: [[Double]] -> TSNEState -> Double
+cost pss st = sumsum $ (zipWith.zipWith) c pss (qdist' (stSolution st))
+    where
+        c p q = -p * log q 
+
diff --git a/test/Data/Algorithm/TSNE/InternalsSpec.hs b/test/Data/Algorithm/TSNE/InternalsSpec.hs
deleted file mode 100644
--- a/test/Data/Algorithm/TSNE/InternalsSpec.hs
+++ /dev/null
@@ -1,97 +0,0 @@
-module Data.Algorithm.TSNE.InternalsSpec (main, spec) where
-
-import Data.Default (def)
-
-import Test.Hspec
-import Data.Algorithm.TSNE.Internals
-import Data.Algorithm.TSNE.Checks
-import Data.Algorithm.TSNE.Types
-import Data.Algorithm.TSNE.Utils
-
--- `main` is here so that this module can be run from GHCi on its own.  It is
--- not needed for automatic spec discovery.
-main :: IO ()
-main = hspec spec
-
-u = undefined
-
-spec :: Spec
-spec = do
-    let n = inputSize testInput
-        w = inputValueSize testInput
-
-    describe "testInput" $ do
-        it "is right shape" $ do
-            testInput `shouldSatisfy` has2DShape (64, 20)
-        it "is valid" $ do
-            inputIsValid testInput `shouldBe` Right ()
-        it "has right size" $ do
-            inputSize testInput `shouldBe` 20
-        it "has right value size" $ do
-            inputValueSize testInput `shouldBe` 64
-
-    describe "initSolution3D" $ do
-        it "is right shape" $
-            --initSolution3D 99 >>= (`shouldSatisfy` (\s -> length s == 3 && all (\xs -> length xs == 99) s))
-            initSolution3D n >>= (`shouldSatisfy` has2DShape (n,3))
-
-    describe "initState" $ do
-        it "is valid state" $
-            initState n >>= (`shouldSatisfy` isRight . (isValidStateForInput testInput))
-
-    describe "neighbourProbabilities" $ do
-        it "is right shape" $ do
-            testNeighbourProbs `shouldSatisfy` has2DShape (n,n)
-
-    describe "qdist" $ do
-        it "is right shape" $ do
-            s <- initSolution3D n
-            qdist s `shouldSatisfy` has2DShape (n,n)        
-
-    describe "qdist'" $ do
-        it "is right shape" $ do
-            s <- initSolution3D n
-            qdist' s `shouldSatisfy` has2DShape (n,n)        
-
-    describe "gradients" $ do
-        it "is right shape" $ do
-            s <- initState n
-            gradients testNeighbourProbs s `shouldSatisfy` has2DShape (n,3)
-
-    describe "stepTSNE" $ do
-        it "works" $ do 
-            s <- initState n
-            stepTSNE def testInput testNeighbourProbs s `shouldSatisfy` isRight . (isValidStateForInput testInput)
-
-
-isRight :: Either a b -> Bool
-isRight (Left _) = False
-isRight (Right _) = True
-
--- first 20 digits from the Python sklearn digits dataset
-testInput :: TSNEInput
-testInput = [
-                [0.0,0.0,5.0,13.0,9.0,1.0,0.0,0.0,0.0,0.0,13.0,15.0,10.0,15.0,5.0,0.0,0.0,3.0,15.0,2.0,0.0,11.0,8.0,0.0,0.0,4.0,12.0,0.0,0.0,8.0,8.0,0.0,0.0,5.0,8.0,0.0,0.0,9.0,8.0,0.0,0.0,4.0,11.0,0.0,1.0,12.0,7.0,0.0,0.0,2.0,14.0,5.0,10.0,12.0,0.0,0.0,0.0,0.0,6.0,13.0,10.0,0.0,0.0,0.0],
-                [0.0,0.0,0.0,12.0,13.0,5.0,0.0,0.0,0.0,0.0,0.0,11.0,16.0,9.0,0.0,0.0,0.0,0.0,3.0,15.0,16.0,6.0,0.0,0.0,0.0,7.0,15.0,16.0,16.0,2.0,0.0,0.0,0.0,0.0,1.0,16.0,16.0,3.0,0.0,0.0,0.0,0.0,1.0,16.0,16.0,6.0,0.0,0.0,0.0,0.0,1.0,16.0,16.0,6.0,0.0,0.0,0.0,0.0,0.0,11.0,16.0,10.0,0.0,0.0],
-                [0.0,0.0,0.0,4.0,15.0,12.0,0.0,0.0,0.0,0.0,3.0,16.0,15.0,14.0,0.0,0.0,0.0,0.0,8.0,13.0,8.0,16.0,0.0,0.0,0.0,0.0,1.0,6.0,15.0,11.0,0.0,0.0,0.0,1.0,8.0,13.0,15.0,1.0,0.0,0.0,0.0,9.0,16.0,16.0,5.0,0.0,0.0,0.0,0.0,3.0,13.0,16.0,16.0,11.0,5.0,0.0,0.0,0.0,0.0,3.0,11.0,16.0,9.0,0.0],
-                [0.0,0.0,7.0,15.0,13.0,1.0,0.0,0.0,0.0,8.0,13.0,6.0,15.0,4.0,0.0,0.0,0.0,2.0,1.0,13.0,13.0,0.0,0.0,0.0,0.0,0.0,2.0,15.0,11.0,1.0,0.0,0.0,0.0,0.0,0.0,1.0,12.0,12.0,1.0,0.0,0.0,0.0,0.0,0.0,1.0,10.0,8.0,0.0,0.0,0.0,8.0,4.0,5.0,14.0,9.0,0.0,0.0,0.0,7.0,13.0,13.0,9.0,0.0,0.0],
-                [0.0,0.0,0.0,1.0,11.0,0.0,0.0,0.0,0.0,0.0,0.0,7.0,8.0,0.0,0.0,0.0,0.0,0.0,1.0,13.0,6.0,2.0,2.0,0.0,0.0,0.0,7.0,15.0,0.0,9.0,8.0,0.0,0.0,5.0,16.0,10.0,0.0,16.0,6.0,0.0,0.0,4.0,15.0,16.0,13.0,16.0,1.0,0.0,0.0,0.0,0.0,3.0,15.0,10.0,0.0,0.0,0.0,0.0,0.0,2.0,16.0,4.0,0.0,0.0],
-                [0.0,0.0,12.0,10.0,0.0,0.0,0.0,0.0,0.0,0.0,14.0,16.0,16.0,14.0,0.0,0.0,0.0,0.0,13.0,16.0,15.0,10.0,1.0,0.0,0.0,0.0,11.0,16.0,16.0,7.0,0.0,0.0,0.0,0.0,0.0,4.0,7.0,16.0,7.0,0.0,0.0,0.0,0.0,0.0,4.0,16.0,9.0,0.0,0.0,0.0,5.0,4.0,12.0,16.0,4.0,0.0,0.0,0.0,9.0,16.0,16.0,10.0,0.0,0.0],
-                [0.0,0.0,0.0,12.0,13.0,0.0,0.0,0.0,0.0,0.0,5.0,16.0,8.0,0.0,0.0,0.0,0.0,0.0,13.0,16.0,3.0,0.0,0.0,0.0,0.0,0.0,14.0,13.0,0.0,0.0,0.0,0.0,0.0,0.0,15.0,12.0,7.0,2.0,0.0,0.0,0.0,0.0,13.0,16.0,13.0,16.0,3.0,0.0,0.0,0.0,7.0,16.0,11.0,15.0,8.0,0.0,0.0,0.0,1.0,9.0,15.0,11.0,3.0,0.0],
-                [0.0,0.0,7.0,8.0,13.0,16.0,15.0,1.0,0.0,0.0,7.0,7.0,4.0,11.0,12.0,0.0,0.0,0.0,0.0,0.0,8.0,13.0,1.0,0.0,0.0,4.0,8.0,8.0,15.0,15.0,6.0,0.0,0.0,2.0,11.0,15.0,15.0,4.0,0.0,0.0,0.0,0.0,0.0,16.0,5.0,0.0,0.0,0.0,0.0,0.0,9.0,15.0,1.0,0.0,0.0,0.0,0.0,0.0,13.0,5.0,0.0,0.0,0.0,0.0],
-                [0.0,0.0,9.0,14.0,8.0,1.0,0.0,0.0,0.0,0.0,12.0,14.0,14.0,12.0,0.0,0.0,0.0,0.0,9.0,10.0,0.0,15.0,4.0,0.0,0.0,0.0,3.0,16.0,12.0,14.0,2.0,0.0,0.0,0.0,4.0,16.0,16.0,2.0,0.0,0.0,0.0,3.0,16.0,8.0,10.0,13.0,2.0,0.0,0.0,1.0,15.0,1.0,3.0,16.0,8.0,0.0,0.0,0.0,11.0,16.0,15.0,11.0,1.0,0.0],
-                [0.0,0.0,11.0,12.0,0.0,0.0,0.0,0.0,0.0,2.0,16.0,16.0,16.0,13.0,0.0,0.0,0.0,3.0,16.0,12.0,10.0,14.0,0.0,0.0,0.0,1.0,16.0,1.0,12.0,15.0,0.0,0.0,0.0,0.0,13.0,16.0,9.0,15.0,2.0,0.0,0.0,0.0,0.0,3.0,0.0,9.0,11.0,0.0,0.0,0.0,0.0,0.0,9.0,15.0,4.0,0.0,0.0,0.0,9.0,12.0,13.0,3.0,0.0,0.0],
-                [0.0,0.0,1.0,9.0,15.0,11.0,0.0,0.0,0.0,0.0,11.0,16.0,8.0,14.0,6.0,0.0,0.0,2.0,16.0,10.0,0.0,9.0,9.0,0.0,0.0,1.0,16.0,4.0,0.0,8.0,8.0,0.0,0.0,4.0,16.0,4.0,0.0,8.0,8.0,0.0,0.0,1.0,16.0,5.0,1.0,11.0,3.0,0.0,0.0,0.0,12.0,12.0,10.0,10.0,0.0,0.0,0.0,0.0,1.0,10.0,13.0,3.0,0.0,0.0],
-                [0.0,0.0,0.0,0.0,14.0,13.0,1.0,0.0,0.0,0.0,0.0,5.0,16.0,16.0,2.0,0.0,0.0,0.0,0.0,14.0,16.0,12.0,0.0,0.0,0.0,1.0,10.0,16.0,16.0,12.0,0.0,0.0,0.0,3.0,12.0,14.0,16.0,9.0,0.0,0.0,0.0,0.0,0.0,5.0,16.0,15.0,0.0,0.0,0.0,0.0,0.0,4.0,16.0,14.0,0.0,0.0,0.0,0.0,0.0,1.0,13.0,16.0,1.0,0.0],
-                [0.0,0.0,5.0,12.0,1.0,0.0,0.0,0.0,0.0,0.0,15.0,14.0,7.0,0.0,0.0,0.0,0.0,0.0,13.0,1.0,12.0,0.0,0.0,0.0,0.0,2.0,10.0,0.0,14.0,0.0,0.0,0.0,0.0,0.0,2.0,0.0,16.0,1.0,0.0,0.0,0.0,0.0,0.0,6.0,15.0,0.0,0.0,0.0,0.0,0.0,9.0,16.0,15.0,9.0,8.0,2.0,0.0,0.0,3.0,11.0,8.0,13.0,12.0,4.0],
-                [0.0,2.0,9.0,15.0,14.0,9.0,3.0,0.0,0.0,4.0,13.0,8.0,9.0,16.0,8.0,0.0,0.0,0.0,0.0,6.0,14.0,15.0,3.0,0.0,0.0,0.0,0.0,11.0,14.0,2.0,0.0,0.0,0.0,0.0,0.0,2.0,15.0,11.0,0.0,0.0,0.0,0.0,0.0,0.0,2.0,15.0,4.0,0.0,0.0,1.0,5.0,6.0,13.0,16.0,6.0,0.0,0.0,2.0,12.0,12.0,13.0,11.0,0.0,0.0],
-                [0.0,0.0,0.0,8.0,15.0,1.0,0.0,0.0,0.0,0.0,1.0,14.0,13.0,1.0,1.0,0.0,0.0,0.0,10.0,15.0,3.0,15.0,11.0,0.0,0.0,7.0,16.0,7.0,1.0,16.0,8.0,0.0,0.0,9.0,16.0,13.0,14.0,16.0,5.0,0.0,0.0,1.0,10.0,15.0,16.0,14.0,0.0,0.0,0.0,0.0,0.0,1.0,16.0,10.0,0.0,0.0,0.0,0.0,0.0,10.0,15.0,4.0,0.0,0.0],
-                [0.0,5.0,12.0,13.0,16.0,16.0,2.0,0.0,0.0,11.0,16.0,15.0,8.0,4.0,0.0,0.0,0.0,8.0,14.0,11.0,1.0,0.0,0.0,0.0,0.0,8.0,16.0,16.0,14.0,0.0,0.0,0.0,0.0,1.0,6.0,6.0,16.0,0.0,0.0,0.0,0.0,0.0,0.0,5.0,16.0,3.0,0.0,0.0,0.0,1.0,5.0,15.0,13.0,0.0,0.0,0.0,0.0,4.0,15.0,16.0,2.0,0.0,0.0,0.0],
-                [0.0,0.0,0.0,8.0,15.0,1.0,0.0,0.0,0.0,0.0,0.0,12.0,14.0,0.0,0.0,0.0,0.0,0.0,3.0,16.0,7.0,0.0,0.0,0.0,0.0,0.0,6.0,16.0,2.0,0.0,0.0,0.0,0.0,0.0,7.0,16.0,16.0,13.0,5.0,0.0,0.0,0.0,15.0,16.0,9.0,9.0,14.0,0.0,0.0,0.0,3.0,14.0,9.0,2.0,16.0,2.0,0.0,0.0,0.0,7.0,15.0,16.0,11.0,0.0],
-                [0.0,0.0,1.0,8.0,15.0,10.0,0.0,0.0,0.0,3.0,13.0,15.0,14.0,14.0,0.0,0.0,0.0,5.0,10.0,0.0,10.0,12.0,0.0,0.0,0.0,0.0,3.0,5.0,15.0,10.0,2.0,0.0,0.0,0.0,16.0,16.0,16.0,16.0,12.0,0.0,0.0,1.0,8.0,12.0,14.0,8.0,3.0,0.0,0.0,0.0,0.0,10.0,13.0,0.0,0.0,0.0,0.0,0.0,0.0,11.0,9.0,0.0,0.0,0.0],
-                [0.0,0.0,10.0,7.0,13.0,9.0,0.0,0.0,0.0,0.0,9.0,10.0,12.0,15.0,2.0,0.0,0.0,0.0,4.0,11.0,10.0,11.0,0.0,0.0,0.0,0.0,1.0,16.0,10.0,1.0,0.0,0.0,0.0,0.0,12.0,13.0,4.0,0.0,0.0,0.0,0.0,0.0,12.0,1.0,12.0,0.0,0.0,0.0,0.0,1.0,10.0,2.0,14.0,0.0,0.0,0.0,0.0,0.0,11.0,14.0,5.0,0.0,0.0,0.0],
-                [0.0,0.0,6.0,14.0,4.0,0.0,0.0,0.0,0.0,0.0,11.0,16.0,10.0,0.0,0.0,0.0,0.0,0.0,8.0,14.0,16.0,2.0,0.0,0.0,0.0,0.0,1.0,12.0,12.0,11.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,11.0,3.0,0.0,0.0,0.0,0.0,0.0,0.0,5.0,11.0,0.0,0.0,0.0,1.0,4.0,4.0,7.0,16.0,2.0,0.0,0.0,7.0,16.0,16.0,13.0,11.0,1.0]
-            ]
-
-testNeighbourProbs :: [[Probability]]
-testNeighbourProbs = neighbourProbabilities def testInput
diff --git a/test/Data/Algorithm/TSNE/PreparationSpec.hs b/test/Data/Algorithm/TSNE/PreparationSpec.hs
new file mode 100644
--- /dev/null
+++ b/test/Data/Algorithm/TSNE/PreparationSpec.hs
@@ -0,0 +1,39 @@
+module Data.Algorithm.TSNE.PreparationSpec (main, spec) where
+
+import Data.Default (def)
+
+import Test.Hspec
+import Data.Algorithm.TSNE.Preparation
+import Data.Algorithm.TSNE.Checks
+import Data.Algorithm.TSNE.Types
+import Data.Algorithm.TSNE.Utils
+import Data.Algorithm.TSNE.TestInput
+import Data.Algorithm.TSNE.TestMisc
+
+-- `main` is here so that this module can be run from GHCi on its own.  It is
+-- not needed for automatic spec discovery.
+main :: IO ()
+main = hspec spec
+
+u = undefined
+
+spec :: Spec
+spec = do
+    let n = inputSize testInput
+        w = inputValueSize testInput
+
+    describe "testInput" $ do
+        it "is right shape" $ do
+            testInput `shouldSatisfy` has2DShape (64, 20)
+        it "is valid" $ do
+            inputIsValid testInput `shouldBe` Right ()
+        it "has right size" $ do
+            inputSize testInput `shouldBe` 20
+        it "has right value size" $ do
+            inputValueSize testInput `shouldBe` 64
+
+    describe "neighbourProbabilities" $ do
+        it "is right shape" $ do
+            testNeighbourProbs `shouldSatisfy` has2DShape (n,n)
+
+
diff --git a/test/Data/Algorithm/TSNE/Run2DSpec.hs b/test/Data/Algorithm/TSNE/Run2DSpec.hs
new file mode 100644
--- /dev/null
+++ b/test/Data/Algorithm/TSNE/Run2DSpec.hs
@@ -0,0 +1,32 @@
+module Data.Algorithm.TSNE.Run2DSpec (main, spec) where
+
+import Data.Default (def)
+
+import Test.Hspec
+import Data.Algorithm.TSNE.Run2D
+import Data.Algorithm.TSNE.Stepping
+import Data.Algorithm.TSNE.Checks
+import Data.Algorithm.TSNE.Types
+import Data.Algorithm.TSNE.Utils
+import Data.Algorithm.TSNE.TestInput
+import Data.Algorithm.TSNE.TestMisc
+
+-- `main` is here so that this module can be run from GHCi on its own.  It is
+-- not needed for automatic spec discovery.
+main :: IO ()
+main = hspec spec
+
+u = undefined
+
+spec :: Spec
+spec = do
+    let n = inputSize testInput
+        w = inputValueSize testInput
+
+    describe "initSolution2D" $ do
+        it "is right shape" $
+            initSolution2D n >>= (`shouldSatisfy` has2DShape (n,2))
+
+    describe "initState" $ do
+        it "is valid state" $
+            initState2D n >>= (`shouldSatisfy` isRight . (isValidStateForInput 2 testInput))
diff --git a/test/Data/Algorithm/TSNE/Run3DSpec.hs b/test/Data/Algorithm/TSNE/Run3DSpec.hs
new file mode 100644
--- /dev/null
+++ b/test/Data/Algorithm/TSNE/Run3DSpec.hs
@@ -0,0 +1,47 @@
+module Data.Algorithm.TSNE.Run3DSpec (main, spec) where
+
+import Data.Default (def)
+
+import Test.Hspec
+import Data.Algorithm.TSNE.Run3D
+import Data.Algorithm.TSNE.Stepping
+import Data.Algorithm.TSNE.Checks
+import Data.Algorithm.TSNE.Types
+import Data.Algorithm.TSNE.Utils
+import Data.Algorithm.TSNE.TestInput
+import Data.Algorithm.TSNE.TestMisc
+
+-- `main` is here so that this module can be run from GHCi on its own.  It is
+-- not needed for automatic spec discovery.
+main :: IO ()
+main = hspec spec
+
+u = undefined
+
+spec :: Spec
+spec = do
+    let n = inputSize testInput
+        w = inputValueSize testInput
+
+    describe "initSolution3D" $ do
+        it "is right shape" $
+            initSolution3D n >>= (`shouldSatisfy` has2DShape (n,3))
+
+    describe "initState" $ do
+        it "is valid state" $
+            initState3D n >>= (`shouldSatisfy` isRight . (isValidStateForInput 3 testInput))
+
+    describe "qdist" $ do
+        it "is right shape" $ do
+            s <- initSolution3D n
+            qdist s `shouldSatisfy` has2DShape (n,n)        
+
+    describe "qdist'" $ do
+        it "is right shape" $ do
+            s <- initSolution3D n
+            qdist' s `shouldSatisfy` has2DShape (n,n)        
+
+    describe "gradients" $ do
+        it "is right shape" $ do
+            s <- initState3D n
+            gradients testNeighbourProbs s `shouldSatisfy` has2DShape (n,3)
diff --git a/test/Data/Algorithm/TSNE/SteppingSpec.hs b/test/Data/Algorithm/TSNE/SteppingSpec.hs
new file mode 100644
--- /dev/null
+++ b/test/Data/Algorithm/TSNE/SteppingSpec.hs
@@ -0,0 +1,29 @@
+module Data.Algorithm.TSNE.SteppingSpec (main, spec) where
+
+import Data.Default (def)
+
+import Test.Hspec
+import Data.Algorithm.TSNE.Stepping
+import Data.Algorithm.TSNE.Run3D
+import Data.Algorithm.TSNE.Checks
+import Data.Algorithm.TSNE.Types
+import Data.Algorithm.TSNE.Utils
+import Data.Algorithm.TSNE.TestInput
+import Data.Algorithm.TSNE.TestMisc
+
+-- `main` is here so that this module can be run from GHCi on its own.  It is
+-- not needed for automatic spec discovery.
+main :: IO ()
+main = hspec spec
+
+u = undefined
+
+spec :: Spec
+spec = do
+    let n = inputSize testInput
+        w = inputValueSize testInput
+
+    describe "stepTSNE" $ do
+        it "works" $ do 
+            s <- initState3D n
+            stepTSNE def testInput testNeighbourProbs s `shouldSatisfy` isRight . (isValidStateForInput 3 testInput)
diff --git a/test/Data/Algorithm/TSNE/TestInput.hs b/test/Data/Algorithm/TSNE/TestInput.hs
new file mode 100644
--- /dev/null
+++ b/test/Data/Algorithm/TSNE/TestInput.hs
@@ -0,0 +1,35 @@
+module Data.Algorithm.TSNE.TestInput where
+
+import Data.Default (def)
+
+import Data.Algorithm.TSNE.Types
+import Data.Algorithm.TSNE.Preparation
+
+-- first 20 digits from the Python sklearn digits dataset
+testInput :: TSNEInput
+testInput = [
+                [0.0,0.0,5.0,13.0,9.0,1.0,0.0,0.0,0.0,0.0,13.0,15.0,10.0,15.0,5.0,0.0,0.0,3.0,15.0,2.0,0.0,11.0,8.0,0.0,0.0,4.0,12.0,0.0,0.0,8.0,8.0,0.0,0.0,5.0,8.0,0.0,0.0,9.0,8.0,0.0,0.0,4.0,11.0,0.0,1.0,12.0,7.0,0.0,0.0,2.0,14.0,5.0,10.0,12.0,0.0,0.0,0.0,0.0,6.0,13.0,10.0,0.0,0.0,0.0],
+                [0.0,0.0,0.0,12.0,13.0,5.0,0.0,0.0,0.0,0.0,0.0,11.0,16.0,9.0,0.0,0.0,0.0,0.0,3.0,15.0,16.0,6.0,0.0,0.0,0.0,7.0,15.0,16.0,16.0,2.0,0.0,0.0,0.0,0.0,1.0,16.0,16.0,3.0,0.0,0.0,0.0,0.0,1.0,16.0,16.0,6.0,0.0,0.0,0.0,0.0,1.0,16.0,16.0,6.0,0.0,0.0,0.0,0.0,0.0,11.0,16.0,10.0,0.0,0.0],
+                [0.0,0.0,0.0,4.0,15.0,12.0,0.0,0.0,0.0,0.0,3.0,16.0,15.0,14.0,0.0,0.0,0.0,0.0,8.0,13.0,8.0,16.0,0.0,0.0,0.0,0.0,1.0,6.0,15.0,11.0,0.0,0.0,0.0,1.0,8.0,13.0,15.0,1.0,0.0,0.0,0.0,9.0,16.0,16.0,5.0,0.0,0.0,0.0,0.0,3.0,13.0,16.0,16.0,11.0,5.0,0.0,0.0,0.0,0.0,3.0,11.0,16.0,9.0,0.0],
+                [0.0,0.0,7.0,15.0,13.0,1.0,0.0,0.0,0.0,8.0,13.0,6.0,15.0,4.0,0.0,0.0,0.0,2.0,1.0,13.0,13.0,0.0,0.0,0.0,0.0,0.0,2.0,15.0,11.0,1.0,0.0,0.0,0.0,0.0,0.0,1.0,12.0,12.0,1.0,0.0,0.0,0.0,0.0,0.0,1.0,10.0,8.0,0.0,0.0,0.0,8.0,4.0,5.0,14.0,9.0,0.0,0.0,0.0,7.0,13.0,13.0,9.0,0.0,0.0],
+                [0.0,0.0,0.0,1.0,11.0,0.0,0.0,0.0,0.0,0.0,0.0,7.0,8.0,0.0,0.0,0.0,0.0,0.0,1.0,13.0,6.0,2.0,2.0,0.0,0.0,0.0,7.0,15.0,0.0,9.0,8.0,0.0,0.0,5.0,16.0,10.0,0.0,16.0,6.0,0.0,0.0,4.0,15.0,16.0,13.0,16.0,1.0,0.0,0.0,0.0,0.0,3.0,15.0,10.0,0.0,0.0,0.0,0.0,0.0,2.0,16.0,4.0,0.0,0.0],
+                [0.0,0.0,12.0,10.0,0.0,0.0,0.0,0.0,0.0,0.0,14.0,16.0,16.0,14.0,0.0,0.0,0.0,0.0,13.0,16.0,15.0,10.0,1.0,0.0,0.0,0.0,11.0,16.0,16.0,7.0,0.0,0.0,0.0,0.0,0.0,4.0,7.0,16.0,7.0,0.0,0.0,0.0,0.0,0.0,4.0,16.0,9.0,0.0,0.0,0.0,5.0,4.0,12.0,16.0,4.0,0.0,0.0,0.0,9.0,16.0,16.0,10.0,0.0,0.0],
+                [0.0,0.0,0.0,12.0,13.0,0.0,0.0,0.0,0.0,0.0,5.0,16.0,8.0,0.0,0.0,0.0,0.0,0.0,13.0,16.0,3.0,0.0,0.0,0.0,0.0,0.0,14.0,13.0,0.0,0.0,0.0,0.0,0.0,0.0,15.0,12.0,7.0,2.0,0.0,0.0,0.0,0.0,13.0,16.0,13.0,16.0,3.0,0.0,0.0,0.0,7.0,16.0,11.0,15.0,8.0,0.0,0.0,0.0,1.0,9.0,15.0,11.0,3.0,0.0],
+                [0.0,0.0,7.0,8.0,13.0,16.0,15.0,1.0,0.0,0.0,7.0,7.0,4.0,11.0,12.0,0.0,0.0,0.0,0.0,0.0,8.0,13.0,1.0,0.0,0.0,4.0,8.0,8.0,15.0,15.0,6.0,0.0,0.0,2.0,11.0,15.0,15.0,4.0,0.0,0.0,0.0,0.0,0.0,16.0,5.0,0.0,0.0,0.0,0.0,0.0,9.0,15.0,1.0,0.0,0.0,0.0,0.0,0.0,13.0,5.0,0.0,0.0,0.0,0.0],
+                [0.0,0.0,9.0,14.0,8.0,1.0,0.0,0.0,0.0,0.0,12.0,14.0,14.0,12.0,0.0,0.0,0.0,0.0,9.0,10.0,0.0,15.0,4.0,0.0,0.0,0.0,3.0,16.0,12.0,14.0,2.0,0.0,0.0,0.0,4.0,16.0,16.0,2.0,0.0,0.0,0.0,3.0,16.0,8.0,10.0,13.0,2.0,0.0,0.0,1.0,15.0,1.0,3.0,16.0,8.0,0.0,0.0,0.0,11.0,16.0,15.0,11.0,1.0,0.0],
+                [0.0,0.0,11.0,12.0,0.0,0.0,0.0,0.0,0.0,2.0,16.0,16.0,16.0,13.0,0.0,0.0,0.0,3.0,16.0,12.0,10.0,14.0,0.0,0.0,0.0,1.0,16.0,1.0,12.0,15.0,0.0,0.0,0.0,0.0,13.0,16.0,9.0,15.0,2.0,0.0,0.0,0.0,0.0,3.0,0.0,9.0,11.0,0.0,0.0,0.0,0.0,0.0,9.0,15.0,4.0,0.0,0.0,0.0,9.0,12.0,13.0,3.0,0.0,0.0],
+                [0.0,0.0,1.0,9.0,15.0,11.0,0.0,0.0,0.0,0.0,11.0,16.0,8.0,14.0,6.0,0.0,0.0,2.0,16.0,10.0,0.0,9.0,9.0,0.0,0.0,1.0,16.0,4.0,0.0,8.0,8.0,0.0,0.0,4.0,16.0,4.0,0.0,8.0,8.0,0.0,0.0,1.0,16.0,5.0,1.0,11.0,3.0,0.0,0.0,0.0,12.0,12.0,10.0,10.0,0.0,0.0,0.0,0.0,1.0,10.0,13.0,3.0,0.0,0.0],
+                [0.0,0.0,0.0,0.0,14.0,13.0,1.0,0.0,0.0,0.0,0.0,5.0,16.0,16.0,2.0,0.0,0.0,0.0,0.0,14.0,16.0,12.0,0.0,0.0,0.0,1.0,10.0,16.0,16.0,12.0,0.0,0.0,0.0,3.0,12.0,14.0,16.0,9.0,0.0,0.0,0.0,0.0,0.0,5.0,16.0,15.0,0.0,0.0,0.0,0.0,0.0,4.0,16.0,14.0,0.0,0.0,0.0,0.0,0.0,1.0,13.0,16.0,1.0,0.0],
+                [0.0,0.0,5.0,12.0,1.0,0.0,0.0,0.0,0.0,0.0,15.0,14.0,7.0,0.0,0.0,0.0,0.0,0.0,13.0,1.0,12.0,0.0,0.0,0.0,0.0,2.0,10.0,0.0,14.0,0.0,0.0,0.0,0.0,0.0,2.0,0.0,16.0,1.0,0.0,0.0,0.0,0.0,0.0,6.0,15.0,0.0,0.0,0.0,0.0,0.0,9.0,16.0,15.0,9.0,8.0,2.0,0.0,0.0,3.0,11.0,8.0,13.0,12.0,4.0],
+                [0.0,2.0,9.0,15.0,14.0,9.0,3.0,0.0,0.0,4.0,13.0,8.0,9.0,16.0,8.0,0.0,0.0,0.0,0.0,6.0,14.0,15.0,3.0,0.0,0.0,0.0,0.0,11.0,14.0,2.0,0.0,0.0,0.0,0.0,0.0,2.0,15.0,11.0,0.0,0.0,0.0,0.0,0.0,0.0,2.0,15.0,4.0,0.0,0.0,1.0,5.0,6.0,13.0,16.0,6.0,0.0,0.0,2.0,12.0,12.0,13.0,11.0,0.0,0.0],
+                [0.0,0.0,0.0,8.0,15.0,1.0,0.0,0.0,0.0,0.0,1.0,14.0,13.0,1.0,1.0,0.0,0.0,0.0,10.0,15.0,3.0,15.0,11.0,0.0,0.0,7.0,16.0,7.0,1.0,16.0,8.0,0.0,0.0,9.0,16.0,13.0,14.0,16.0,5.0,0.0,0.0,1.0,10.0,15.0,16.0,14.0,0.0,0.0,0.0,0.0,0.0,1.0,16.0,10.0,0.0,0.0,0.0,0.0,0.0,10.0,15.0,4.0,0.0,0.0],
+                [0.0,5.0,12.0,13.0,16.0,16.0,2.0,0.0,0.0,11.0,16.0,15.0,8.0,4.0,0.0,0.0,0.0,8.0,14.0,11.0,1.0,0.0,0.0,0.0,0.0,8.0,16.0,16.0,14.0,0.0,0.0,0.0,0.0,1.0,6.0,6.0,16.0,0.0,0.0,0.0,0.0,0.0,0.0,5.0,16.0,3.0,0.0,0.0,0.0,1.0,5.0,15.0,13.0,0.0,0.0,0.0,0.0,4.0,15.0,16.0,2.0,0.0,0.0,0.0],
+                [0.0,0.0,0.0,8.0,15.0,1.0,0.0,0.0,0.0,0.0,0.0,12.0,14.0,0.0,0.0,0.0,0.0,0.0,3.0,16.0,7.0,0.0,0.0,0.0,0.0,0.0,6.0,16.0,2.0,0.0,0.0,0.0,0.0,0.0,7.0,16.0,16.0,13.0,5.0,0.0,0.0,0.0,15.0,16.0,9.0,9.0,14.0,0.0,0.0,0.0,3.0,14.0,9.0,2.0,16.0,2.0,0.0,0.0,0.0,7.0,15.0,16.0,11.0,0.0],
+                [0.0,0.0,1.0,8.0,15.0,10.0,0.0,0.0,0.0,3.0,13.0,15.0,14.0,14.0,0.0,0.0,0.0,5.0,10.0,0.0,10.0,12.0,0.0,0.0,0.0,0.0,3.0,5.0,15.0,10.0,2.0,0.0,0.0,0.0,16.0,16.0,16.0,16.0,12.0,0.0,0.0,1.0,8.0,12.0,14.0,8.0,3.0,0.0,0.0,0.0,0.0,10.0,13.0,0.0,0.0,0.0,0.0,0.0,0.0,11.0,9.0,0.0,0.0,0.0],
+                [0.0,0.0,10.0,7.0,13.0,9.0,0.0,0.0,0.0,0.0,9.0,10.0,12.0,15.0,2.0,0.0,0.0,0.0,4.0,11.0,10.0,11.0,0.0,0.0,0.0,0.0,1.0,16.0,10.0,1.0,0.0,0.0,0.0,0.0,12.0,13.0,4.0,0.0,0.0,0.0,0.0,0.0,12.0,1.0,12.0,0.0,0.0,0.0,0.0,1.0,10.0,2.0,14.0,0.0,0.0,0.0,0.0,0.0,11.0,14.0,5.0,0.0,0.0,0.0],
+                [0.0,0.0,6.0,14.0,4.0,0.0,0.0,0.0,0.0,0.0,11.0,16.0,10.0,0.0,0.0,0.0,0.0,0.0,8.0,14.0,16.0,2.0,0.0,0.0,0.0,0.0,1.0,12.0,12.0,11.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,11.0,3.0,0.0,0.0,0.0,0.0,0.0,0.0,5.0,11.0,0.0,0.0,0.0,1.0,4.0,4.0,7.0,16.0,2.0,0.0,0.0,7.0,16.0,16.0,13.0,11.0,1.0]
+            ]
+
+testNeighbourProbs :: [[Probability]]
+testNeighbourProbs = neighbourProbabilities def testInput
+
diff --git a/test/Data/Algorithm/TSNE/TestMisc.hs b/test/Data/Algorithm/TSNE/TestMisc.hs
new file mode 100644
--- /dev/null
+++ b/test/Data/Algorithm/TSNE/TestMisc.hs
@@ -0,0 +1,5 @@
+module Data.Algorithm.TSNE.TestMisc where
+
+isRight :: Either a b -> Bool
+isRight (Left _) = False
+isRight (Right _) = True
diff --git a/tsne.cabal b/tsne.cabal
--- a/tsne.cabal
+++ b/tsne.cabal
@@ -1,5 +1,5 @@
 name:                tsne
-version:             1.2.0
+version:             1.3.0
 synopsis:            t-SNE
 description:         Pure Haskell implementation of the t-SNE dimension reduction algorithm.
 homepage:            https://bitbucket.org/robagar/haskell-tsne
@@ -16,9 +16,12 @@
   hs-source-dirs:      src
   exposed-modules:     Data.Algorithm.TSNE,
                        Data.Algorithm.TSNE.Types, 
-                       Data.Algorithm.TSNE.Internals, 
                        Data.Algorithm.TSNE.Utils,
-                       Data.Algorithm.TSNE.Checks 
+                       Data.Algorithm.TSNE.Checks, 
+                       Data.Algorithm.TSNE.Preparation, 
+                       Data.Algorithm.TSNE.Stepping, 
+                       Data.Algorithm.TSNE.Run2D, 
+                       Data.Algorithm.TSNE.Run3D 
   build-depends:       base >= 4.7 && < 5,
                        data-default,
                        deepseq,
@@ -31,9 +34,14 @@
   hs-source-dirs:      test
   main-is:             Spec.hs
   other-modules:       Data.Algorithm.TSNESpec,
-                       Data.Algorithm.TSNE.InternalsSpec,
+                       Data.Algorithm.TSNE.PreparationSpec,
+                       Data.Algorithm.TSNE.SteppingSpec,
                        Data.Algorithm.TSNE.UtilsSpec,
-                       Data.Algorithm.TSNE.ChecksSpec
+                       Data.Algorithm.TSNE.ChecksSpec,
+                       Data.Algorithm.TSNE.Run3DSpec,
+                       Data.Algorithm.TSNE.Run2DSpec,
+                       Data.Algorithm.TSNE.TestInput,
+                       Data.Algorithm.TSNE.TestMisc
   build-depends:       base,
                        hspec,
                        data-default,