diff --git a/hmm-lapack.cabal b/hmm-lapack.cabal
--- a/hmm-lapack.cabal
+++ b/hmm-lapack.cabal
@@ -1,5 +1,5 @@
 Name:                hmm-lapack
-Version:             0.3.0.3
+Version:             0.4
 Synopsis:            Hidden Markov Models using LAPACK primitives
 Description:
   Hidden Markov Models implemented using LAPACK data types and operations.
@@ -41,7 +41,7 @@
   Changes.md
 
 Source-Repository this
-  Tag:         0.3.0.3
+  Tag:         0.4
   Type:        darcs
   Location:    http://hub.darcs.net/thielema/hmm-lapack
 
@@ -67,19 +67,18 @@
     Math.HiddenMarkovModel.Utility
     Math.HiddenMarkovModel.CSV
   Build-Depends:
-    lapack >=0.2.2 && <0.3,
+    lapack >=0.3 && <0.4,
     fixed-length >=0.2.1 && <0.3,
     tfp >=1.0 && <1.1,
     netlib-ffi >=0.1.1 && <0.2,
-    comfort-array >=0.2 && <0.4,
+    comfort-array >=0.4 && <0.5,
     QuickCheck >=2.5 && <3,
     explicit-exception >=0.1.7 && <0.2,
-    boxes >=0.1.5 && <0.2,
     lazy-csv >=0.5 && <0.6,
     random >=1.0 && <1.2,
     transformers >= 0.2 && <0.6,
-    non-empty >=0.2.1 && <0.4,
-    semigroups >=0.17 && <0.19,
+    non-empty >=0.3.2 && <0.4,
+    semigroups >=0.17 && <1.0,
     containers >=0.4.2 && <0.7,
     utility-ht >=0.0.12 && <0.1,
     deepseq >=1.3 && <1.5,
diff --git a/src/Math/HiddenMarkovModel.hs b/src/Math/HiddenMarkovModel.hs
--- a/src/Math/HiddenMarkovModel.hs
+++ b/src/Math/HiddenMarkovModel.hs
@@ -26,17 +26,18 @@
 import Math.HiddenMarkovModel.Private
           (T(..), Trained(..), mergeTrained, toCells, parseCSV)
 import Math.HiddenMarkovModel.Utility
-          (SquareMatrix, squareConstant, distance,
+          (squareConstant, distance, matrixDistance,
            randomItemProp, normalizeProb, attachOnes)
 
+import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix
 import qualified Numeric.LAPACK.Matrix as Matrix
 import qualified Numeric.LAPACK.Vector as Vector
+import Numeric.LAPACK.Matrix ((#!))
 
 import qualified Numeric.Netlib.Class as Class
 
 import qualified Data.Array.Comfort.Storable as StorableArray
 import qualified Data.Array.Comfort.Shape as Shape
-import qualified Data.Array.Comfort.Boxed as Array
 
 import qualified Text.CSV.Lazy.String as CSV
 
@@ -53,12 +54,12 @@
 
 
 type DiscreteTrained symbol sh prob =
-         Trained (Distr.DiscreteTrained symbol sh prob) sh prob
-type Discrete symbol sh prob = T (Distr.Discrete symbol sh prob) sh prob
+         Trained (Distr.Discrete symbol) sh prob
+type Discrete symbol sh prob = T (Distr.Discrete symbol) sh prob
 
 type GaussianTrained emiSh stateSh a =
-         Trained (Distr.GaussianTrained emiSh stateSh a) stateSh a
-type Gaussian emiSh stateSh a = T (Distr.Gaussian emiSh stateSh a) stateSh a
+         Trained (Distr.Gaussian emiSh) stateSh a
+type Gaussian emiSh stateSh a = T (Distr.Gaussian emiSh) stateSh a
 
 
 {- |
@@ -68,9 +69,8 @@
 You can use this as a starting point for 'Normalized.trainUnsupervised'.
 -}
 uniform ::
-   (Distr.Info distr, Distr.StateShape distr ~ sh, Shape.C sh,
-    Distr.Probability distr ~ prob) =>
-   distr -> T distr sh prob
+   (Distr.Info typ, Shape.C sh, Class.Real prob) =>
+   Distr.T typ sh prob -> T typ sh prob
 uniform distr =
    let sh = Distr.statesShape distr
        c = recip $ fromIntegral $ Shape.size sh
@@ -82,31 +82,29 @@
 
 
 probabilitySequence ::
-   (Traversable f, Distr.EmissionProb distr,
-    Distr.StateShape distr ~ sh, Shape.Indexed sh, Shape.Index sh ~ state,
-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>
-   T distr sh prob -> f (state, emission) -> f prob
+   (Distr.EmissionProb typ, Shape.Indexed sh, Shape.Index sh ~ state,
+    Class.Real prob, Distr.Emission typ prob ~ emission, Traversable f) =>
+   T typ sh prob -> f (state, emission) -> f prob
 probabilitySequence hmm =
    snd
    .
    mapAccumL
       (\index (s, e) ->
-         ((transition hmm StorableArray.!) . flip (,) s,
+         ((transition hmm #!) . flip (,) s,
           index s * Distr.emissionStateProb (distribution hmm) e s))
       (initial hmm StorableArray.!)
 
 generate ::
-   (Rnd.RandomGen g, Ord prob, Rnd.Random prob, Distr.Generate distr,
-    Distr.StateShape distr ~ sh, Shape.Indexed sh, Shape.Index sh ~ state,
-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>
-   T distr sh prob -> g -> [emission]
+   (Distr.Generate typ, Shape.Indexed sh, Class.Real prob,
+    Rnd.RandomGen g, Rnd.Random prob, Distr.Emission typ prob ~ emission) =>
+   T typ sh prob -> g -> [emission]
 generate hmm = map snd . generateLabeled hmm
 
 generateLabeled ::
-   (Rnd.RandomGen g, Ord prob, Rnd.Random prob, Distr.Generate distr,
-    Distr.StateShape distr ~ sh, Shape.Indexed sh, Shape.Index sh ~ state,
-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>
-   T distr sh prob -> g -> [(state, emission)]
+   (Distr.Generate typ, Shape.Indexed sh, Shape.Index sh ~ state,
+    Rnd.RandomGen g, Rnd.Random prob,
+    Class.Real prob, Distr.Emission typ prob ~ emission) =>
+   T typ sh prob -> g -> [(state, emission)]
 generateLabeled hmm =
    MS.evalState $
    flip MS.evalStateT (initial hmm) $
@@ -123,31 +121,24 @@
 Contribute a manually labeled emission sequence to a HMM training.
 -}
 trainSupervised ::
-   (Distr.StateShape distr ~ sh, Shape.Index sh ~ state,
-    Distr.Estimate tdistr distr,
-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>
-   sh -> NonEmpty.T [] (state, emission) -> Trained tdistr sh prob
+   (Distr.Estimate typ, Shape.Indexed sh, Shape.Index sh ~ state,
+    Class.Real prob, Distr.Emission typ prob ~ emission) =>
+   sh -> NonEmpty.T [] (state, emission) -> Trained typ sh prob
 trainSupervised sh xs =
    let getState (s, _x) = s
    in  Trained {
-          trainedInitial =
-             StorableArray.fromAssociations sh 0
-                [(getState (NonEmpty.head xs), 1)],
+          trainedInitial = Vector.unit sh $ getState $ NonEmpty.head xs,
           trainedTransition =
-             Matrix.transpose $
-             StorableArray.accumulate (+) (squareConstant sh 0) $
+             Matrix.transpose $ ArrMatrix.fromVector $
+             StorableArray.accumulate (+)
+                (ArrMatrix.toVector $ squareConstant sh 0) $
              attachOnes $ NonEmpty.mapAdjacent (,) $ fmap getState xs,
-          trainedDistribution =
-             Distr.accumulateEmissions $ Array.map attachOnes $
-             Array.accumulate (flip (:))
-                (Array.fromList sh $ replicate (Shape.size sh) [])
-                (NonEmpty.flatten xs)
+          trainedDistribution = Distr.accumulateEmissions sh xs
        }
 
 finishTraining ::
-   (Shape.C sh, Eq sh,
-    Distr.Estimate tdistr distr, Distr.Probability distr ~ prob) =>
-   Trained tdistr sh prob -> T distr sh prob
+   (Distr.Estimate typ, Shape.C sh, Eq sh, Class.Real prob) =>
+   Trained typ sh prob -> T typ sh prob
 finishTraining hmm =
    Cons {
       initial = normalizeProb $ trainedInitial hmm,
@@ -156,18 +147,15 @@
    }
 
 normalizeProbColumns ::
-   (Shape.C sh, Eq sh, Class.Real a) => SquareMatrix sh a -> SquareMatrix sh a
+   (Shape.C sh, Eq sh, Class.Real a) => Matrix.Square sh a -> Matrix.Square sh a
 normalizeProbColumns m =
    Matrix.scaleColumns (StorableArray.map recip (Matrix.columnSums m)) m
 
 trainMany ::
-   (Shape.C sh, Eq sh,
-    Distr.Estimate tdistr distr, Distr.Probability distr ~ prob,
-    Foldable f) =>
-   (trainingData -> Trained tdistr sh prob) ->
-   NonEmpty.T f trainingData -> T distr sh prob
-trainMany train =
-   finishTraining . NonEmpty.foldl1Map mergeTrained train
+   (Distr.Estimate typ, Shape.C sh, Eq sh, Class.Real prob, Foldable f) =>
+   (trainingData -> Trained typ sh prob) ->
+   NonEmpty.T f trainingData -> T typ sh prob
+trainMany train = finishTraining . NonEmpty.foldl1Map mergeTrained train
 
 
 
@@ -181,25 +169,22 @@
 should suffice for defining an abort criterion.
 -}
 deviation ::
-   (Shape.C sh, Eq sh, Class.Real prob, Ord prob) =>
-   T distr sh prob -> T distr sh prob -> prob
+   (Shape.C sh, Eq sh, Class.Real prob) =>
+   T typ sh prob -> T typ sh prob -> prob
 deviation hmm0 hmm1 =
    distance (initial hmm0) (initial hmm1)
    `max`
-   distance (transition hmm0) (transition hmm1)
+   matrixDistance (transition hmm0) (transition hmm1)
 
 
 toCSV ::
-   (Distr.ToCSV distr, Shape.Indexed sh, Class.Real prob, Show prob) =>
-   T distr sh prob -> String
+   (Distr.ToCSV typ, Shape.Indexed sh, Class.Real prob, Show prob) =>
+   T typ sh prob -> String
 toCSV hmm =
-   CSV.ppCSVTable $ snd $ CSV.toCSVTable $ HMMCSV.padTable "" $
-   toCells hmm
+   CSV.ppCSVTable $ snd $ CSV.toCSVTable $ HMMCSV.padTable "" $ toCells hmm
 
 fromCSV ::
-   (Distr.FromCSV distr, Distr.StateShape distr ~ stateSh,
-    Shape.Indexed stateSh, Shape.Index stateSh ~ state,
-    Class.Real prob, Read prob) =>
-   (Int -> stateSh) -> String -> ME.Exceptional String (T distr stateSh prob)
+   (Distr.FromCSV typ, Shape.Indexed sh, Eq sh, Class.Real prob, Read prob) =>
+   (Int -> sh) -> String -> ME.Exceptional String (T typ sh prob)
 fromCSV makeShape =
    MS.evalStateT (parseCSV makeShape) . map HMMCSV.fixShortRow . CSV.parseCSV
diff --git a/src/Math/HiddenMarkovModel/CSV.hs b/src/Math/HiddenMarkovModel/CSV.hs
--- a/src/Math/HiddenMarkovModel/CSV.hs
+++ b/src/Math/HiddenMarkovModel/CSV.hs
@@ -1,16 +1,15 @@
 module Math.HiddenMarkovModel.CSV where
 
-import Math.HiddenMarkovModel.Utility (SquareMatrix, vectorDim)
+import Math.HiddenMarkovModel.Utility (vectorDim)
 
 import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape
 import qualified Numeric.LAPACK.Matrix as Matrix
 import qualified Numeric.LAPACK.Vector as Vector
-import Numeric.LAPACK.Matrix (ZeroInt)
+import Numeric.LAPACK.Matrix (ShapeInt)
 import Numeric.LAPACK.Vector (Vector)
 
 import qualified Numeric.Netlib.Class as Class
 
-import qualified Data.Array.Comfort.Storable as ComfortArray
 import qualified Data.Array.Comfort.Shape as Shape
 
 import qualified Text.CSV.Lazy.String as CSV
@@ -32,7 +31,7 @@
 cellsFromVector = map show . Vector.toList
 
 cellsFromSquare ::
-   (Shape.Indexed sh, Show a, Class.Real a) => SquareMatrix sh a -> [[String]]
+   (Shape.Indexed sh, Show a, Class.Real a) => Matrix.Square sh a -> [[String]]
 cellsFromSquare = map (map show . Vector.toList) . Matrix.toRows
 
 padTable :: a -> [[a]] -> [[a]]
@@ -110,20 +109,21 @@
 
 parseVectorCells ::
    (Read a, Class.Real a) =>
-   CSVParser (Vector ZeroInt a)
+   CSVParser (Vector ShapeInt a)
 parseVectorCells =
    parseVectorFields =<< getRow
 
+-- ToDo: Maybe check row consistency already here?
 parseVectorFields ::
    (Read a, Class.Real a) =>
-   CSV.CSVRow -> CSVParser (Vector ZeroInt a)
+   CSV.CSVRow -> CSVParser (Vector ShapeInt a)
 parseVectorFields =
    MT.lift . fmap Vector.autoFromList . mapM parseNumberCell .
    Rev.dropWhile (null . CSV.csvFieldContent)
 
 parseNonEmptyVectorCells ::
    (Read a, Class.Real a) =>
-   CSVParser (Vector ZeroInt a)
+   CSVParser (Vector ShapeInt a)
 parseNonEmptyVectorCells = do
    v <- parseVectorCells
    assert (vectorDim v > 0) "no data for vector"
@@ -143,14 +143,14 @@
 
 parseSquareMatrixCells ::
    (Shape.C sh, Read a, Class.Real a) =>
-   sh -> CSVParser (SquareMatrix sh a)
+   sh -> CSVParser (Matrix.Square sh a)
 parseSquareMatrixCells sh = do
    let n = Shape.size sh
    rows <- replicateM n parseVectorCells
    assert (not $ null rows) "no rows"
    assert (all ((n==) . vectorDim) rows) "inconsistent matrix dimensions"
    return $
-      ComfortArray.reshape (MatrixShape.square MatrixShape.RowMajor sh) $
+      Matrix.reshape (MatrixShape.square MatrixShape.RowMajor sh) $
       Matrix.fromRows (Shape.ZeroBased n) rows
 
 parseStringList :: CSV.CSVRow -> CSVParser [String]
diff --git a/src/Math/HiddenMarkovModel/Distribution.hs b/src/Math/HiddenMarkovModel/Distribution.hs
--- a/src/Math/HiddenMarkovModel/Distribution.hs
+++ b/src/Math/HiddenMarkovModel/Distribution.hs
@@ -1,28 +1,34 @@
 {-# LANGUAGE TypeFamilies #-}
-{-# LANGUAGE FlexibleContexts #-}
-{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE EmptyDataDecls #-}
 module Math.HiddenMarkovModel.Distribution (
-   Emission, Probability, StateShape,
-   Info(..), Generate(..), EmissionProb(..), Estimate(..),
+   T(..), Trained(..), Emission,
+   Show(..), NFData(..), Format(..),
+   Info(..), Generate(..), EmissionProb(..),
+   Estimate(..), accumulateEmissionVectors,
 
-   Discrete(..), DiscreteTrained(..),
-   Gaussian(..), GaussianTrained(..), gaussian,
+   Discrete, discreteFromList,
+   Gaussian, gaussian, gaussianTrained,
 
    ToCSV(..), FromCSV(..), HMMCSV.CSVParser, CSVSymbol(..),
    ) where
 
 import qualified Math.HiddenMarkovModel.CSV as HMMCSV
-import Math.HiddenMarkovModel.Utility (SquareMatrix, randomItemProp, vectorDim)
+import Math.HiddenMarkovModel.Utility (randomItemProp, vectorDim)
 
 import qualified Numeric.LAPACK.Matrix.HermitianPositiveDefinite as HermitianPD
 import qualified Numeric.LAPACK.Matrix.Hermitian as Hermitian
 import qualified Numeric.LAPACK.Matrix.Triangular as Triangular
 import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape
+import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix
 import qualified Numeric.LAPACK.Matrix as Matrix
 import qualified Numeric.LAPACK.Vector as Vector
-import Numeric.LAPACK.Matrix ((<#))
+import qualified Numeric.LAPACK.Format as Format
+import qualified Numeric.LAPACK.Output as Output
+import Numeric.LAPACK.Matrix ((-*#), (-/#), (#/\), (|*-), (#!))
 import Numeric.LAPACK.Vector (Vector)
-import Numeric.LAPACK.Format (FormatArray, Format(format))
+import Numeric.LAPACK.Format (FormatArray)
+import Numeric.LAPACK.Output (Output)
 
 import qualified Numeric.Netlib.Class as Class
 import Foreign.Storable (Storable)
@@ -31,334 +37,397 @@
 import qualified Data.Array.Comfort.Shape as Shape
 import qualified Data.Array.Comfort.Boxed as Array
 import Data.Array.Comfort.Boxed (Array, (!))
+import Data.Array.Comfort.Shape ((:+:)((:+:)))
 
 import qualified System.Random as Rnd
 
 import qualified Text.CSV.Lazy.String as CSV
-import qualified Text.PrettyPrint.Boxes as TextBox
-import Text.PrettyPrint.Boxes ((<>), (<+>))
 import Text.Read.HT (maybeRead)
 import Text.Printf (printf)
 
 import qualified Control.Monad.Exception.Synchronous as ME
 import qualified Control.Monad.Trans.Class as MT
 import qualified Control.Monad.Trans.State as MS
-import Control.DeepSeq (NFData, rnf)
+import qualified Control.DeepSeq as DeepSeq
 import Control.Monad (liftM2)
-import Control.Applicative (liftA2, (<|>))
+import Control.Applicative (liftA2)
 
+import qualified Data.NonEmpty.Map as NonEmptyMap
 import qualified Data.NonEmpty as NonEmpty
-import qualified Data.Foldable as Fold
+import qualified Data.Semigroup as Sg
 import qualified Data.Map as Map
-import qualified Data.Set as Set
 import qualified Data.List.HT as ListHT
 import qualified Data.List as List
 import Data.Functor.Identity (Identity(Identity), runIdentity)
-import Data.Foldable (Foldable, foldMap)
-import Data.Tuple.HT (mapFst, fst3, swap)
-import Data.Monoid (Endo(Endo, appEndo))
-import Data.Map (Map)
-import Data.Maybe (fromMaybe, listToMaybe)
+import Data.Tuple.HT (snd3)
+import Data.Set (Set)
+import Data.Maybe (listToMaybe)
 
-import Prelude ()
-import Prelude2010
+import qualified Prelude as P
+import Prelude2010 hiding (Show, showsPrec)
 
 
-type HermitianMatrix sh = Hermitian.Hermitian sh
-type UpperTriangular sh = Triangular.Upper sh
 
+data family T typ sh prob
+data family Trained typ sh prob
 
-type family Probability distr
-type family Emission distr
-type family StateShape distr
+type family Emission typ prob
 
 
-class (Class.Real (Probability distr)) => Info distr where
-   statesShape :: distr -> StateShape distr
+class Show typ where
+   showsPrec ::
+      (Shape.C sh, P.Show sh, P.Show prob, Storable prob) =>
+      Int -> T typ sh prob -> ShowS
+   showsPrecTrained ::
+      (Shape.C sh, P.Show sh, P.Show prob, Storable prob) =>
+      Int -> Trained typ sh prob -> ShowS
 
-class (Class.Real (Probability distr)) => Generate distr where
+instance
+   (Show typ, Shape.C sh, P.Show sh, P.Show prob, Storable prob) =>
+      P.Show (T typ sh prob) where
+   showsPrec = showsPrec
+
+instance
+   (Show typ, Shape.C sh, P.Show sh, P.Show prob, Storable prob) =>
+      P.Show (Trained typ sh prob) where
+   showsPrec = showsPrecTrained
+
+
+class NFData typ where
+   rnf ::
+      (DeepSeq.NFData sh, DeepSeq.NFData prob, Shape.C sh) =>
+      T typ sh prob -> ()
+   rnfTrained ::
+      (DeepSeq.NFData sh, DeepSeq.NFData prob, Shape.C sh) =>
+      Trained typ sh prob -> ()
+
+instance
+   (NFData typ, DeepSeq.NFData sh, DeepSeq.NFData prob, Shape.C sh) =>
+      DeepSeq.NFData (T typ sh prob) where
+   rnf = rnf
+
+instance
+   (NFData typ, DeepSeq.NFData sh, DeepSeq.NFData prob, Shape.C sh) =>
+      DeepSeq.NFData (Trained typ sh prob) where
+   rnf = rnfTrained
+
+
+class Format typ where
+   format ::
+      (Shape.C sh, Output out, Class.Real prob) =>
+      String -> T typ sh prob -> out
+
+instance
+   (Format typ, Shape.C sh, Class.Real prob) =>
+      Format.Format (T typ sh prob) where
+   format = format
+
+
+
+class Info typ where
+   statesShape :: (Shape.C sh) => T typ sh prob -> sh
+   statesShapeTrained :: (Shape.C sh) => Trained typ sh prob -> sh
+
+class Generate typ where
    generate ::
-      (Rnd.RandomGen g, Emission distr ~ emission, StateShape distr ~ sh) =>
-      distr -> Shape.Index sh -> MS.State g emission
+      (Shape.Indexed sh, Class.Real prob, Rnd.Random prob, Rnd.RandomGen g) =>
+      T typ sh prob -> Shape.Index sh -> MS.State g (Emission typ prob)
 
-class
-   (Shape.Indexed (StateShape distr), Class.Real (Probability distr)) =>
-      EmissionProb distr where
+class EmissionProb typ where
+   mapStatesShape ::
+      (Shape.C sh0, Shape.C sh1) =>
+      (sh0 -> sh1) -> T typ sh0 prob -> T typ sh1 prob
    {-
    This function could be implemented generically in terms of emissionStateProb
    but that would require an Info constraint.
    -}
    emissionProb ::
-      distr -> Emission distr -> Vector (StateShape distr) (Probability distr)
+      (Shape.C sh, Class.Real prob) =>
+      T typ sh prob -> Emission typ prob -> Vector sh prob
    emissionStateProb ::
-      distr -> Emission distr -> Shape.Index (StateShape distr) -> Probability distr
+      (Shape.Indexed sh, Class.Real prob) =>
+      T typ sh prob -> Emission typ prob -> Shape.Index sh -> prob
    emissionStateProb distr e s = emissionProb distr e StorableArray.! s
 
-class
-   (Distribution tdistr ~ distr, Trained distr ~ tdistr, EmissionProb distr) =>
-      Estimate tdistr distr where
-   type Distribution tdistr
-   type Trained distr
+class (EmissionProb typ) => Estimate typ where
    accumulateEmissions ::
-      (Probability distr ~ prob, StateShape distr ~ sh) =>
-      Array sh [(Emission distr, prob)] -> tdistr
-   -- could as well be in Semigroup class
-   combine :: tdistr -> tdistr -> tdistr
-   normalize :: tdistr -> distr
+      (Shape.Indexed sh, Class.Real prob, Shape.Index sh ~ state) =>
+      sh -> NonEmpty.T [] (state, Emission typ prob) -> Trained typ sh prob
+   trainVector ::
+      (Shape.C sh, Eq sh, Class.Real prob) =>
+      Emission typ prob -> Vector sh prob -> Trained typ sh prob
+   combine ::
+      (Shape.C sh, Eq sh, Class.Real prob) =>
+      Trained typ sh prob -> Trained typ sh prob -> Trained typ sh prob
+   normalize ::
+      (Shape.C sh, Eq sh, Class.Real prob) =>
+      Trained typ sh prob -> T typ sh prob
 
+accumulateEmissionVectors ::
+   (Estimate typ, Shape.C sh, Eq sh, Class.Real prob) =>
+   NonEmpty.T [] (Emission typ prob, Vector sh prob) -> Trained typ sh prob
+accumulateEmissionVectors = NonEmpty.foldl1Map combine (uncurry trainVector)
 
-newtype Discrete symbol sh prob = Discrete (Map symbol (Vector sh prob))
-   deriving (Show)
+instance
+   (Estimate typ, Shape.C sh, Eq sh, Class.Real prob) =>
+      Sg.Semigroup (Trained typ sh prob) where
+   (<>) = combine
 
-newtype
-   DiscreteTrained symbol sh prob =
-      DiscreteTrained (Map symbol (Vector sh prob))
-   deriving (Show)
 
-type instance Probability (Discrete symbol sh prob) = prob
-type instance Emission (Discrete symbol sh prob) = symbol
-type instance StateShape (Discrete symbol sh prob) = sh
+data Discrete symbol
 
+newtype instance T (Discrete symbol) sh prob =
+      Discrete (Matrix.General (Set symbol) sh prob)
 
-instance
-   (NFData sh, NFData prob, NFData symbol) =>
-      NFData (Discrete symbol sh prob) where
-   rnf (Discrete m) = rnf m
+newtype instance Trained (Discrete symbol) sh prob =
+      DiscreteTrained (NonEmptyMap.T symbol (Vector sh prob))
 
-instance
-   (NFData sh, NFData prob, NFData symbol) =>
-      NFData (DiscreteTrained symbol sh prob) where
-   rnf (DiscreteTrained m) = rnf m
+type instance Emission (Discrete symbol) prob = symbol
 
-instance
-   (FormatArray sh, Class.Real prob, Format symbol) =>
-      Format (Discrete symbol sh prob) where
+
+instance (P.Show symbol, Ord symbol) => Show (Discrete symbol) where
+   showsPrec prec (Discrete m) = P.showsPrec prec m
+   showsPrecTrained prec (DiscreteTrained m) = P.showsPrec prec m
+
+instance (DeepSeq.NFData symbol) => NFData (Discrete symbol) where
+   rnf (Discrete m) = DeepSeq.rnf m
+   rnfTrained (DiscreteTrained m) = DeepSeq.rnf m
+
+instance (P.Show symbol, Ord symbol) => Format (Discrete symbol) where
    format fmt (Discrete m) =
-      TextBox.vsep 1 TextBox.left $
-      map (\(sym,v) -> format fmt sym <> TextBox.char ':' <+> format fmt v) $
-      Map.toAscList m
+      Output.formatAligned $
+      map (\(sym,v) ->
+            map (Identity . Output.text) $
+            (show sym ++ ":") : map (printFmt fmt) (Vector.toList v)) $
+      Array.toAssociations $ Matrix.toRowArray m
 
-instance
-   (Shape.C sh, Class.Real prob, Ord symbol) =>
-      Info (Discrete symbol sh prob) where
-   statesShape (Discrete m) = StorableArray.shape $ snd $ Map.findMin m
+-- cf. Data.Bifunctor.Flip
+newtype Flip f b a = Flip {getFlip :: f a b}
 
-instance
-   (Shape.Indexed sh, Class.Real prob, Ord symbol, Ord prob, Rnd.Random prob) =>
-      Generate (Discrete symbol sh prob) where
-   generate (Discrete m) state =
-      randomItemProp $ Map.toAscList $ fmap (StorableArray.! state) m
+printFmt :: (Class.Real a) => String -> a -> String
+printFmt fmt =
+   getFlip $ Class.switchReal (Flip $ printf fmt) (Flip $ printf fmt)
 
-instance
-   (Shape.Indexed sh, Class.Real prob, Ord symbol) =>
-      EmissionProb (Discrete symbol sh prob) where
-   emissionProb (Discrete m) =
-      mapLookup "emitDiscrete: unknown emission symbol" m
+instance (Ord symbol) => Info (Discrete symbol) where
+   statesShape (Discrete m) = Matrix.width m
+   statesShapeTrained (DiscreteTrained m) = discreteStateShape m
 
-instance
-   (Shape.Indexed sh, Eq sh, Class.Real prob, Ord symbol) =>
-      Estimate (DiscreteTrained symbol sh prob) (Discrete symbol sh prob) where
-   type Distribution (DiscreteTrained symbol sh prob) = Discrete symbol sh prob
-   type Trained (Discrete symbol sh prob) = DiscreteTrained symbol sh prob
-   accumulateEmissions grouped =
-      let set = Set.toAscList $ foldMap (Set.fromList . map fst) grouped
-          emi = Map.fromAscList $ zip set [0..]
-      in  DiscreteTrained $ Map.fromAscList $ zip set $
-          transposeVectorList $
-          Array.map
-             (StorableArray.accumulate (+)
-                 (Vector.constant (Shape.ZeroBased $ length set) 0) .
-              map (mapFst
-                 (mapLookup "estimateDiscrete: unknown emission symbol" emi)))
-             grouped
+instance (Ord symbol) => Generate (Discrete symbol) where
+   generate (Discrete m) =
+      randomItemProp . StorableArray.toAssociations . Matrix.takeColumn m
+
+instance (Ord symbol) => EmissionProb (Discrete symbol) where
+   mapStatesShape f (Discrete m) = Discrete $ Matrix.mapWidth f m
+   emissionProb (Discrete m) = Matrix.takeRow m
+   emissionStateProb (Discrete m) x s = m #! (x,s)
+
+instance (Ord symbol) => Estimate (Discrete symbol) where
+   accumulateEmissions sh =
+      DiscreteTrained .
+      NonEmptyMap.map
+         (StorableArray.reshape sh .
+          StorableArray.fromAssociations 0 (Shape.Deferred sh) .
+          Map.toList) .
+      NonEmptyMap.fromListWith (Map.unionWith (+)) .
+      fmap (\(state,sym) -> (sym, Map.singleton (Shape.deferIndex sh state) 1))
+   trainVector sym = DiscreteTrained . NonEmptyMap.singleton sym
    combine (DiscreteTrained distr0) (DiscreteTrained distr1) =
-      DiscreteTrained $ Map.unionWith Vector.add distr0 distr1
+      DiscreteTrained $ NonEmptyMap.unionWith Vector.add distr0 distr1
    normalize (DiscreteTrained distr) =
-      Discrete $ if Map.null distr then distr else normalizeProbVecs distr
+      Discrete $ normalizeProbColumns $ discreteFromMap distr
 
-transposeVectorList ::
-   (Shape.C sh, Eq sh, Class.Real a) =>
-   Array sh (Vector Matrix.ZeroInt a) -> [Vector sh a]
-transposeVectorList xs =
-   case Array.toList xs of
-      [] -> []
-      x:_ -> Matrix.toRows $ Matrix.fromColumnArray (StorableArray.shape x) xs
+normalizeProbColumns ::
+   (Shape.C height, Shape.C width, Eq width, Class.Real a) =>
+   Matrix.General height width a -> Matrix.General height width a
+normalizeProbColumns m = m #/\ Matrix.columnSums m
 
-normalizeProbVecs ::
-   (Shape.C sh, Eq sh, Foldable f, Functor f, Class.Real a) =>
-   f (Vector sh a) -> f (Vector sh a)
-normalizeProbVecs vs =
-   let factors =
-         StorableArray.map recip $ List.foldl1' Vector.add $ Fold.toList vs
-   in fmap (Vector.mul factors) vs
+discreteStateShape ::
+   (Shape.C sh) => NonEmptyMap.T symbol (Vector sh prob) -> sh
+discreteStateShape =
+   StorableArray.shape . snd . fst . NonEmptyMap.minViewWithKey
 
-mapLookup :: (Ord k) => String -> Map.Map k a -> k -> a
-mapLookup msg dict x = Map.findWithDefault (error msg) x dict
+discreteFromMap ::
+   (Ord symbol, Shape.C sh, Eq sh, Class.Real prob) =>
+   NonEmptyMap.T symbol (Vector sh prob) -> Matrix.General (Set symbol) sh prob
+discreteFromMap m =
+   Matrix.fromRowArray (discreteStateShape m) $
+   Array.fromMap $ NonEmptyMap.flatten m
 
+discreteFromList ::
+   (Ord symbol, Shape.C sh, Eq sh, Class.Real prob) =>
+   NonEmpty.T [] (symbol, Vector sh prob) -> T (Discrete symbol) sh prob
+discreteFromList = Discrete . discreteFromMap . NonEmptyMap.fromList
 
-newtype Gaussian emiSh stateSh a =
-      Gaussian (Array stateSh (Vector emiSh a, UpperTriangular emiSh a, a))
-   deriving (Show)
 
-newtype GaussianTrained emiSh stateSh a =
+
+data Gaussian emiSh
+
+newtype instance T (Gaussian emiSh) stateSh a =
+   Gaussian (Array stateSh (a, Vector emiSh a, Triangular.Upper emiSh a))
+
+newtype instance Trained (Gaussian emiSh) stateSh a =
    GaussianTrained
-      (Array stateSh
-         (Maybe (Vector emiSh a, HermitianMatrix emiSh a, a)))
-   deriving (Show)
+      (StorableArray.Array (stateSh, MatrixShape.Hermitian (():+:emiSh)) a)
 
-type instance Probability (Gaussian emiSh stateSh a) = a
-type instance Emission (Gaussian emiSh stateSh a) = Vector emiSh a
-type instance StateShape (Gaussian emiSh stateSh a) = stateSh
+type instance Emission (Gaussian emiSh) a = Vector emiSh a
 
 
-instance
-   (NFData emiSh, NFData stateSh, Shape.C stateSh, NFData a, Storable a) =>
-      NFData (Gaussian emiSh stateSh a) where
-   rnf (Gaussian params) = rnf params
+instance (Shape.C emiSh, P.Show emiSh) => Show (Gaussian emiSh) where
+   showsPrec prec (Gaussian m) = P.showsPrec prec m
+   showsPrecTrained prec (GaussianTrained m) = P.showsPrec prec m
 
-instance
-   (NFData emiSh, NFData stateSh, Shape.C stateSh, NFData a, Storable a) =>
-      NFData (GaussianTrained emiSh stateSh a) where
-   rnf (GaussianTrained params) = rnf params
+instance (DeepSeq.NFData emiSh) => NFData (Gaussian emiSh) where
+   rnf (Gaussian params) = DeepSeq.rnf params
+   rnfTrained (GaussianTrained params) = DeepSeq.rnf params
 
 
-instance
-   (FormatArray emiSh, Shape.C stateSh, Class.Real a) =>
-      Format (Gaussian emiSh stateSh a) where
+instance (FormatArray emiSh) => Format (Gaussian emiSh) where
    format = runFormatGaussian $ Class.switchReal formatGaussian formatGaussian
 
-newtype FormatGaussian emiSh stateSh a =
+newtype FormatGaussian out emiSh stateSh a =
    FormatGaussian
-      {runFormatGaussian :: String -> Gaussian emiSh stateSh a -> TextBox.Box}
+      {runFormatGaussian :: String -> T (Gaussian emiSh) stateSh a -> out}
 
 formatGaussian ::
-   (FormatArray emiSh, Shape.C stateSh, Class.Real a, Format a) =>
-   FormatGaussian emiSh stateSh a
+   (FormatArray emiSh, Shape.C stateSh,
+    Class.Real a, Format.Format a, Output out) =>
+   FormatGaussian out emiSh stateSh a
 formatGaussian =
-   FormatGaussian $ \fmt (Gaussian params) -> format fmt $ Array.toList params
+   FormatGaussian $ \fmt (Gaussian params) ->
+      Format.format fmt $ Array.toList params
 
 
-instance
-   (Shape.Indexed stateSh, Eq stateSh, Class.Real a) =>
-      Info (Gaussian emiSh stateSh a) where
+instance Info (Gaussian emiSh) where
    statesShape (Gaussian params) = Array.shape params
+   statesShapeTrained (GaussianTrained params) =
+      fst $ StorableArray.shape params
 
-instance
-   (Shape.C emiSh, Eq emiSh, Shape.Indexed stateSh, Eq stateSh, Class.Real a) =>
-      Generate (Gaussian emiSh stateSh a) where
+instance (Shape.C emiSh, Eq emiSh) => Generate (Gaussian emiSh) where
    generate (Gaussian allParams) state = do
-      let (center, covarianceChol, _c) = allParams ! state
+      let (_c, center, covarianceChol) = allParams ! state
       seed <- MS.state Rnd.random
       return $
          Vector.add center $
          Vector.random Vector.Normal (StorableArray.shape center) seed
-            <# covarianceChol
+            -*# covarianceChol
 
-instance
-   (Shape.C emiSh, Eq emiSh, Shape.Indexed stateSh, Eq stateSh, Class.Real a) =>
-      EmissionProb (Gaussian emiSh stateSh a) where
+instance (Shape.C emiSh, Eq emiSh) => EmissionProb (Gaussian emiSh) where
+   mapStatesShape f (Gaussian m) = Gaussian $ Array.mapShape f m
    emissionProb (Gaussian allParams) x =
-      Vector.fromList (Array.shape allParams) $
-      map (emissionProbGen x) $ Array.toList allParams
+      StorableArray.fromBoxed $ fmap (gaussianEmissionProb x) allParams
    emissionStateProb (Gaussian allParams) x s =
-      emissionProbGen x $ allParams ! s
+      gaussianEmissionProb x $ allParams ! s
 
-emissionProbGen ::
+gaussianEmissionProb ::
    (Shape.C emiSh, Eq emiSh, Class.Real a) =>
-   Vector emiSh a -> (Vector emiSh a, UpperTriangular emiSh a, a) -> a
-emissionProbGen x (center, covarianceChol, c) =
-   let x0 =
-         Matrix.solveVector (Triangular.transpose covarianceChol) $
-         Vector.sub x center
-   in  c * cexp ((-1/2) * Vector.inner x0 x0)
+   Vector emiSh a -> (a, Vector emiSh a, Triangular.Upper emiSh a) -> a
+gaussianEmissionProb x (c, center, covarianceChol) =
+   c * expSquared (Vector.sub x center -/# covarianceChol)
 
+expSquared :: (Shape.C sh, Class.Real a) => Vector sh a -> a
+expSquared =
+   getNorm $ Class.switchReal (Norm expSquaredAux) (Norm expSquaredAux)
 
-instance
-   (Shape.C emiSh, Eq emiSh, Shape.Indexed stateSh, Eq stateSh, Class.Real a) =>
-      Estimate
-         (GaussianTrained emiSh stateSh a)
-         (Gaussian emiSh stateSh a) where
-   type Distribution (GaussianTrained emiSh stateSh a) =
-            Gaussian emiSh stateSh a
-   type Trained (Gaussian emiSh stateSh a) = GaussianTrained emiSh stateSh a
-   accumulateEmissions =
-      let params xs =
-            (NonEmpty.foldl1Map Vector.add (uncurry $ flip Vector.scale) xs,
-             covarianceReal $ fmap swap xs,
-             Fold.sum $ fmap snd xs)
-      in  GaussianTrained . fmap (fmap params . NonEmpty.fetch)
-   combine (GaussianTrained distr0) (GaussianTrained distr1) =
-      let comb (center0, covariance0, weight0)
-               (center1, covariance1, weight1) =
-             (Vector.add center0 center1,
-              Vector.add covariance0 covariance1,
-              weight0 + weight1)
-      in  GaussianTrained $ Array.zipWith (maybePlus comb) distr0 distr1
+newtype Norm f a = Norm {getNorm :: f a -> a}
+
+expSquaredAux ::
+   (Shape.C sh, Class.Floating a, Vector.RealOf a ~ ar, Class.Real ar) =>
+   Vector sh a -> ar
+expSquaredAux x = exp ((-1/2) * Vector.norm2Squared x)
+
+
+instance (Shape.C emiSh, Eq emiSh) => Estimate (Gaussian emiSh) where
+   accumulateEmissions sh xs =
+      let emiSh = StorableArray.shape $ snd $ NonEmpty.head xs
+          hermSh = MatrixShape.hermitian MatrixShape.RowMajor (():+:emiSh)
+      in GaussianTrained $
+         Matrix.toRowMajor . Matrix.fromRowArray hermSh . Array.reshape sh .
+         Array.accumulate Vector.add
+            (Array.replicate (Shape.Deferred sh) (Vector.zero hermSh)) .
+         map (\(state,v) -> (Shape.deferIndex sh state, extendedHermitian v)) .
+         NonEmpty.flatten
+            $ xs
+   trainVector xs probs =
+      GaussianTrained $ Matrix.toRowMajor $ probs |*- extendedHermitian xs
+   combine (GaussianTrained m0) (GaussianTrained m1) =
+      GaussianTrained $ Vector.add m0 m1
    {-
      Sum_i (xi-m) * (xi-m)^T
    = Sum_i xi*xi^T + Sum_i m*m^T - Sum_i xi*m^T - Sum_i m*xi^T
    = Sum_i xi*xi^T - Sum_i m*m^T
    = Sum_i xi*xi^T - n * m*m^T
    -}
-   normalize (GaussianTrained distr) =
-      let params (centerSum, covarianceSum, weight) =
-             let c = recip weight
-                 center = Vector.scale c centerSum
+   normalize (GaussianTrained m) =
+      let params (weight, centerSum, covarianceSum) =
+             let c = recip (weight#!((),()))
+                 center = Vector.scale c $ Matrix.flattenRow centerSum
              in  (center,
-                  Vector.sub (Vector.scale c covarianceSum)
+                  Matrix.sub
+                     (Matrix.scaleRealReal c covarianceSum)
                      (Hermitian.outer MatrixShape.RowMajor center))
-      in  Gaussian $
-          fmap
-             (gaussianParameters . params .
-              fromMaybe
-                (error "Distribution.normalize: undefined array element")) $
-          distr
-
-maybePlus :: (a -> a -> a) -> Maybe a -> Maybe a -> Maybe a
-maybePlus f mx my = liftA2 f mx my <|> mx <|> my
+      in Gaussian $
+         fmap (gaussianParameters . params .
+               Hermitian.split . ArrMatrix.fromVector) $
+         Matrix.toRowArray $ Matrix.fromRowMajor m
 
+extendedHermitian ::
+   (Shape.C emiSh, Class.Floating a) =>
+   StorableArray.Array emiSh a ->
+   StorableArray.Array (MatrixShape.Hermitian (():+:emiSh)) a
+extendedHermitian =
+   ArrMatrix.toVector .
+   Hermitian.outer MatrixShape.RowMajor . Vector.append (Vector.one ())
 
-newtype CovarianceReal f emiSh a =
-   CovarianceReal
-      {getCovarianceReal :: f (a, Vector emiSh a) -> HermitianMatrix emiSh a}
+{- |
+input array must be non-empty
+-}
+gaussianTrained ::
+   (Shape.C emiSh, Eq emiSh, Shape.C stateSh, Class.Real prob) =>
+   Array stateSh (prob, Vector emiSh prob, Matrix.Hermitian emiSh prob) ->
+   Trained (Gaussian emiSh) stateSh prob
+gaussianTrained =
+   GaussianTrained . Matrix.toRowMajor .
+   matrixFromRowArray "HMM.Distribution.gaussianTrained" .
+   fmap
+      (\(weight, center, covariance) ->
+         ArrMatrix.toVector $
+         Hermitian.stack
+            (Hermitian.fromList MatrixShape.RowMajor () [weight])
+            (Matrix.singleRow MatrixShape.RowMajor center)
+            covariance)
 
-covarianceReal ::
-   (Shape.C emiSh, Eq emiSh, Class.Real a) =>
-   NonEmpty.T [] (a, Vector emiSh a) -> HermitianMatrix emiSh a
-covarianceReal =
-   getCovarianceReal $
-   Class.switchReal
-      (CovarianceReal $ Hermitian.sumRank1NonEmpty MatrixShape.RowMajor)
-      (CovarianceReal $ Hermitian.sumRank1NonEmpty MatrixShape.RowMajor)
+matrixFromRowArray ::
+   (Shape.C width, Eq width, Shape.C height, Class.Real a) =>
+   String ->
+   Array height (StorableArray.Array width a) ->
+   Matrix.General height width a
+matrixFromRowArray name xs =
+   case Array.toList xs of
+      [] -> error $ name ++ ": empty array"
+      x:_ -> Matrix.fromRowArray (StorableArray.shape x) xs
 
 gaussian ::
    (Shape.C emiSh, Shape.C stateSh, Class.Real prob) =>
-   Array stateSh (Vector emiSh prob, HermitianMatrix emiSh prob) ->
-   Gaussian emiSh stateSh prob
-gaussian = consGaussian . fmap gaussianParameters
+   Array stateSh (Vector emiSh prob, Matrix.Hermitian emiSh prob) ->
+   T (Gaussian emiSh) stateSh prob
+gaussian = Gaussian . fmap gaussianParameters
 
 gaussianParameters ::
    (Shape.C emiSh, Class.Real prob) =>
-   (Vector emiSh prob, HermitianMatrix emiSh prob) ->
-   (Vector emiSh prob, UpperTriangular emiSh prob, prob)
+   (Vector emiSh prob, Matrix.Hermitian emiSh prob) ->
+   (prob, Vector emiSh prob, Triangular.Upper emiSh prob)
 gaussianParameters (center, covariance) =
    gaussianFromCholesky center $ HermitianPD.decompose covariance
 
-consGaussian ::
-   (Shape.C stateSh) =>
-   Array stateSh (Vector emiSh a, UpperTriangular emiSh a, a) ->
-   Gaussian emiSh stateSh a
-consGaussian = Gaussian
-
 gaussianFromCholesky ::
    (Shape.C emiSh, Class.Real prob) =>
-   Vector emiSh prob -> UpperTriangular emiSh prob ->
-   (Vector emiSh prob, UpperTriangular emiSh prob, prob)
+   Vector emiSh prob -> Triangular.Upper emiSh prob ->
+   (prob, Vector emiSh prob, Triangular.Upper emiSh prob)
 gaussianFromCholesky center covarianceChol =
    let covarianceSqrtDet =
          Vector.product $ Triangular.takeDiagonal covarianceChol
-   in  (center, covarianceChol,
-        recip (sqrt2pi ^ vectorDim center * covarianceSqrtDet))
+   in  (recip (sqrt2pi ^ vectorDim center * covarianceSqrtDet),
+        center, covarianceChol)
 
 sqrt2pi :: (Class.Real a) => a
 sqrt2pi = runIdentity $ Class.switchReal sqrt2piAux sqrt2piAux
@@ -366,16 +435,16 @@
 sqrt2piAux :: (Floating a) => Identity a
 sqrt2piAux = Identity $ sqrt (2*pi)
 
-cexp :: (Class.Real a) => a -> a
-cexp = appEndo $ Class.switchReal (Endo exp) (Endo exp)
 
-
-
-class ToCSV distr where
-   toCells :: distr -> [[String]]
+class ToCSV typ where
+   toCells ::
+      (Shape.C sh, Class.Real prob, P.Show prob) =>
+      T typ sh prob -> [[String]]
 
-class FromCSV distr where
-   parseCells :: StateShape distr -> HMMCSV.CSVParser distr
+class FromCSV typ where
+   parseCells ::
+      (Shape.C sh, Eq sh, Class.Real prob, Read prob) =>
+      sh -> HMMCSV.CSVParser (T typ sh prob)
 
 class (Ord symbol) => CSVSymbol symbol where
    cellFromSymbol :: symbol -> String
@@ -390,21 +459,18 @@
    symbolFromCell = maybeRead
 
 
-instance
-   (Shape.C sh, Class.Real prob, Show prob, Read prob, CSVSymbol symbol) =>
-      ToCSV (Discrete symbol sh prob) where
+instance (CSVSymbol symbol) => ToCSV (Discrete symbol) where
    toCells (Discrete m) =
       map
          (\(symbol, probs) ->
             cellFromSymbol symbol : HMMCSV.cellsFromVector probs) $
-      Map.toAscList m
+      Array.toAssociations $ Matrix.toRowArray m
 
-instance
-   (Shape.C sh, Class.Real prob, Show prob, Read prob, CSVSymbol symbol) =>
-      FromCSV (Discrete symbol sh prob) where
+instance (CSVSymbol symbol) => FromCSV (Discrete symbol) where
    parseCells n =
-      fmap (Discrete . Map.fromList) $
-      HMMCSV.manyRowsUntilEnd $ parseSymbolProb n
+      let p = parseSymbolProb n
+      in fmap discreteFromList $
+         liftA2 NonEmpty.Cons (HMMCSV.getRow >>= p) (HMMCSV.manyRowsUntilEnd p)
 
 parseSymbolProb ::
    (Shape.C sh, Class.Real prob, Read prob, CSVSymbol symbol) =>
@@ -426,36 +492,30 @@
                 return $ StorableArray.reshape sh v)
 
 
-instance
-   (Shape.Indexed emiSh, Shape.Indexed stateSh,
-    Class.Real a, Eq a, Show a, Read a) =>
-      ToCSV (Gaussian emiSh stateSh a) where
+instance (Shape.Indexed emiSh) => ToCSV (Gaussian emiSh) where
    toCells (Gaussian params) =
       List.intercalate [[]] $
       map
-         (\(center, covarianceChol, _) ->
+         (\(_, center, covarianceChol) ->
             HMMCSV.cellsFromVector center :
             HMMCSV.cellsFromSquare (Triangular.toSquare covarianceChol)) $
       Array.toList params
 
-instance
-   (emiSh ~ Matrix.ZeroInt, Shape.Indexed stateSh,
-    Class.Real a, Eq a, Show a, Read a) =>
-      FromCSV (Gaussian emiSh stateSh a) where
+instance (emiSh ~ Matrix.ShapeInt) => FromCSV (Gaussian emiSh) where
    parseCells sh = do
       let n = Shape.size sh
       gs <- HMMCSV.manySepUntilEnd parseSingleGaussian
       HMMCSV.assert (length gs == n) $
          printf "number of states (%d) and number of Gaussians (%d) mismatch"
             n (length gs)
-      let sizes = map (vectorDim . fst3) gs
+      let sizes = map (vectorDim . snd3) gs
       HMMCSV.assert (ListHT.allEqual sizes) $
          printf "dimensions of emissions mismatch: %s" (show sizes)
-      return $ consGaussian $ Array.fromList sh gs
+      return $ Gaussian $ Array.fromList sh gs
 
 parseSingleGaussian ::
-   (emiSh ~ Matrix.ZeroInt, Class.Real prob, Eq prob, Read prob) =>
-   HMMCSV.CSVParser (Vector emiSh prob, UpperTriangular emiSh prob, prob)
+   (emiSh ~ Matrix.ShapeInt, Class.Real prob, Eq prob, Read prob) =>
+   HMMCSV.CSVParser (prob, Vector emiSh prob, Triangular.Upper emiSh prob)
 parseSingleGaussian = do
    center <- HMMCSV.parseNonEmptyVectorCells
    covarianceCholSquare <-
@@ -475,6 +535,8 @@
 -}
 isUpperTriang ::
    (Shape.C sh, Class.Real a, Eq a) =>
-   SquareMatrix sh a -> UpperTriangular sh a -> Bool
+   Matrix.Square sh a -> Triangular.Upper sh a -> Bool
 isUpperTriang m mt =
-   Vector.toList m == Vector.toList (Triangular.toSquare mt)
+   Vector.toList (ArrMatrix.toVector m)
+   ==
+   Vector.toList (ArrMatrix.toVector (Triangular.toSquare mt))
diff --git a/src/Math/HiddenMarkovModel/Example/CirclePrivate.hs b/src/Math/HiddenMarkovModel/Example/CirclePrivate.hs
--- a/src/Math/HiddenMarkovModel/Example/CirclePrivate.hs
+++ b/src/Math/HiddenMarkovModel/Example/CirclePrivate.hs
@@ -46,7 +46,7 @@
 hmm :: HMM
 hmm =
    HMM.Cons {
-      HMM.initial = normalizeProb $ Vector.constant stateSet 1,
+      HMM.initial = normalizeProb $ Vector.one stateSet,
       HMM.transition =
          squareFromLists stateSet $
             stateVector 0.9 0.0 0.0 0.1 :
diff --git a/src/Math/HiddenMarkovModel/Example/SineWave.hs b/src/Math/HiddenMarkovModel/Example/SineWave.hs
--- a/src/Math/HiddenMarkovModel/Example/SineWave.hs
+++ b/src/Math/HiddenMarkovModel/Example/SineWave.hs
@@ -40,7 +40,7 @@
 hmm :: HMM
 hmm =
    HMM.Cons {
-      HMM.initial = normalizeProb $ Vector.constant stateSet 1,
+      HMM.initial = normalizeProb $ Vector.one stateSet,
       HMM.transition =
          squareFromLists stateSet $
             stateVector 0.9 0.0 0.0 0.1 :
diff --git a/src/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs b/src/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs
--- a/src/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs
+++ b/src/Math/HiddenMarkovModel/Example/TrafficLightPrivate.hs
@@ -7,17 +7,14 @@
 
 import qualified Numeric.LAPACK.Vector as Vector
 import Numeric.LAPACK.Vector (Vector)
-import Numeric.LAPACK.Format (Format(format))
 
 import qualified Data.Array.Comfort.Shape as Shape
 
-import qualified Text.PrettyPrint.Boxes as TextBox
 import Text.Read.HT (maybeRead)
 
 import Control.DeepSeq (NFData(rnf))
 import Control.Monad (liftM2)
 
-import qualified Data.Map as Map
 import qualified Data.NonEmpty as NonEmpty
 import qualified Data.List.HT as ListHT
 import Data.NonEmpty ((!:))
@@ -31,9 +28,6 @@
    rnf Red = ()
    rnf _ = ()
 
-instance Format Color where
-   format _fmt = TextBox.text . show
-
 {- |
 Using 'show' and 'read' is not always a good choice
 since they must format and parse Haskell expressions
@@ -67,9 +61,9 @@
             stateVector 0.0 0.0 0.2 0.8 :
             [],
       HMM.distribution =
-         Distr.Discrete $ Map.fromList $
-            (Red,    stateVector 1 0 0 0) :
-            (Yellow, stateVector 0 1 0 1):
+         Distr.discreteFromList $
+            (Red,    stateVector 1 0 0 0) !:
+            (Yellow, stateVector 0 1 0 1) :
             (Green,  stateVector 0 0 1 0) :
             []
    }
@@ -86,8 +80,8 @@
             stateVector 0.2 0.2 0.3 0.3 :
             [],
       HMM.distribution =
-         Distr.Discrete $ Map.fromList $
-            (Red,    stateVector 0.6 0.2 0.2 0.2) :
+         Distr.discreteFromList $
+            (Red,    stateVector 0.6 0.2 0.2 0.2) !:
             (Yellow, stateVector 0.2 0.6 0.2 0.6) :
             (Green,  stateVector 0.2 0.2 0.6 0.2) :
             []
diff --git a/src/Math/HiddenMarkovModel/Named.hs b/src/Math/HiddenMarkovModel/Named.hs
--- a/src/Math/HiddenMarkovModel/Named.hs
+++ b/src/Math/HiddenMarkovModel/Named.hs
@@ -1,6 +1,4 @@
 {-# LANGUAGE TypeFamilies #-}
-{-# LANGUAGE FlexibleContexts #-}
-{-# LANGUAGE UndecidableInstances #-}
 module Math.HiddenMarkovModel.Named (
    T(..),
    Discrete,
@@ -42,30 +40,31 @@
 Although 'nameFromStateMap' and 'stateFromNameMap' are exported
 you must be careful to keep them consistent when you alter them.
 -}
-data T distr sh ix prob =
+data T typ sh ix prob =
    Cons {
-      model :: HMM.T distr sh prob,
+      model :: HMM.T typ sh prob,
       nameFromStateMap :: Array sh String,
       stateFromNameMap :: Map String ix
    }
    deriving (Show)
 
+type Simple typ sh prob = T typ sh (Shape.Index sh) prob
 type Discrete symbol stateSh prob =
-      T (Distr.Discrete symbol stateSh prob) stateSh (Shape.Index stateSh) prob
+      Simple (Distr.Discrete symbol) stateSh prob
 type Gaussian emiSh stateSh a =
-      T (Distr.Gaussian emiSh stateSh a) stateSh (Shape.Index stateSh) a
+      Simple (Distr.Gaussian emiSh) stateSh a
 
 
 instance
-   (NFData distr, NFData sh, NFData ix, NFData prob,
+   (Distr.NFData typ, NFData sh, NFData ix, NFData prob,
     Shape.C sh, Storable prob) =>
-      NFData (T distr sh ix prob) where
+      NFData (T typ sh ix prob) where
    rnf hmm = rnf (model hmm, nameFromStateMap hmm, stateFromNameMap hmm)
 
 
 fromModelAndNames ::
-   (Shape.Indexed sh, Shape.Index sh ~ state) =>
-   HMM.T distr sh prob -> [String] -> T distr sh state prob
+   (Shape.Indexed sh) =>
+   HMM.T typ sh prob -> [String] -> Simple typ sh prob
 fromModelAndNames md names =
    let m = Array.fromList (StorableArray.shape $ HMM.initial md) names
    in  Cons {
@@ -82,26 +81,24 @@
 
 
 toCSV ::
-   (Distr.ToCSV distr, Shape.Indexed sh, Class.Real prob, Show prob) =>
-   T distr sh ix prob -> String
+   (Distr.ToCSV typ, Shape.Indexed sh, Class.Real prob, Show prob) =>
+   Simple typ sh prob -> String
 toCSV hmm =
    CSV.ppCSVTable $ snd $ CSV.toCSVTable $ HMMCSV.padTable "" $
       Array.toList (nameFromStateMap hmm) : HMM.toCells (model hmm)
 
 fromCSV ::
-   (Distr.FromCSV distr, Distr.StateShape distr ~ stateSh,
-    Shape.Indexed stateSh, Shape.Index stateSh ~ state,
+   (Distr.FromCSV typ, Shape.Indexed stateSh, Eq stateSh,
     Class.Real prob, Read prob) =>
    (Int -> stateSh) ->
-   String -> ME.Exceptional String (T distr stateSh state prob)
+   String -> ME.Exceptional String (Simple typ stateSh prob)
 fromCSV makeShape =
    MS.evalStateT (parseCSV makeShape) . map HMMCSV.fixShortRow . CSV.parseCSV
 
 parseCSV ::
-   (Distr.FromCSV distr, Distr.StateShape distr ~ stateSh,
-    Shape.Indexed stateSh, Shape.Index stateSh ~ state,
+   (Distr.FromCSV typ, Shape.Indexed stateSh, Eq stateSh,
     Class.Real prob, Read prob) =>
-   (Int -> stateSh) -> HMMCSV.CSVParser (T distr stateSh state prob)
+   (Int -> stateSh) -> HMMCSV.CSVParser (Simple typ stateSh prob)
 parseCSV makeShape = do
    names <- HMMCSV.parseStringList =<< HMMCSV.getRow
    let duplicateNames =
diff --git a/src/Math/HiddenMarkovModel/Normalized.hs b/src/Math/HiddenMarkovModel/Normalized.hs
--- a/src/Math/HiddenMarkovModel/Normalized.hs
+++ b/src/Math/HiddenMarkovModel/Normalized.hs
@@ -10,12 +10,12 @@
 import qualified Math.HiddenMarkovModel.Distribution as Distr
 import Math.HiddenMarkovModel.Private
           (T(..), Trained(..), emission,
-           biscaleTransition, matrixMaxMul, sumTransitions)
-import Math.HiddenMarkovModel.Utility
-         (SquareMatrix, normalizeFactor, normalizeProb)
+           biscaleTransition, revealGen, sumTransitions)
+import Math.HiddenMarkovModel.Utility (normalizeFactor, normalizeProb)
 
+import qualified Numeric.LAPACK.Matrix as Matrix
 import qualified Numeric.LAPACK.Vector as Vector
-import Numeric.LAPACK.Matrix ((<#), (#>))
+import Numeric.LAPACK.Matrix ((-*#), (#*|))
 import Numeric.LAPACK.Vector (Vector)
 
 import qualified Numeric.Netlib.Class as Class
@@ -23,15 +23,12 @@
 import qualified Control.Functor.HT as Functor
 
 import qualified Data.Array.Comfort.Storable as StorableArray
-import qualified Data.Array.Comfort.Boxed as Array
 import qualified Data.Array.Comfort.Shape as Shape
 
 import qualified Data.NonEmpty.Class as NonEmptyC
 import qualified Data.NonEmpty as NonEmpty
 import qualified Data.Foldable as Fold
-import qualified Data.List as List
-import Data.Traversable (Traversable, mapAccumL)
-import Data.Tuple.HT (mapFst, mapSnd, swap)
+import Data.Traversable (Traversable)
 
 
 {- |
@@ -40,43 +37,43 @@
 that it may be rounded to zero in the choosen number type.
 -}
 logLikelihood ::
-   (Distr.EmissionProb distr, Distr.StateShape distr ~ sh, Eq sh, Floating prob,
-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,
+   (Distr.EmissionProb typ, Shape.C sh, Eq sh, Floating prob,
+    Class.Real prob, Distr.Emission typ prob ~ emission,
     Traversable f) =>
-   T distr sh prob -> NonEmpty.T f emission -> prob
+   T typ sh prob -> NonEmpty.T f emission -> prob
 logLikelihood hmm = Fold.sum . fmap (log . fst) . alpha hmm
 
 alpha ::
-   (Distr.EmissionProb distr, Distr.StateShape distr ~ sh, Eq sh,
-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,
+   (Distr.EmissionProb typ, Shape.C sh, Eq sh,
+    Class.Real prob, Distr.Emission typ prob ~ emission,
     Traversable f) =>
-   T distr sh prob ->
+   T typ sh prob ->
    NonEmpty.T f emission -> NonEmpty.T f (prob, Vector sh prob)
 alpha hmm (NonEmpty.Cons x xs) =
    let normMulEmiss y = normalizeFactor . Vector.mul (emission hmm y)
    in  NonEmpty.scanl
-          (\(_,alphai) xi -> normMulEmiss xi (transition hmm #> alphai))
+          (\(_,alphai) xi -> normMulEmiss xi (transition hmm #*| alphai))
           (normMulEmiss x (initial hmm))
           xs
 
 beta ::
-   (Distr.EmissionProb distr, Distr.StateShape distr ~ sh, Eq sh,
-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,
+   (Distr.EmissionProb typ, Shape.C sh, Eq sh,
+    Class.Real prob, Distr.Emission typ prob ~ emission,
     Traversable f, NonEmptyC.Reverse f) =>
-   T distr sh prob ->
+   T typ sh prob ->
    f (prob, emission) -> NonEmpty.T f (Vector sh prob)
 beta hmm =
    nonEmptyScanr
       (\(ci,xi) betai ->
          Vector.scale (recip ci) $
-         Vector.mul (emission hmm xi) betai <# transition hmm)
-      (Vector.constant (StorableArray.shape $ initial hmm) 1)
+         Vector.mul (emission hmm xi) betai -*# transition hmm)
+      (Vector.one $ StorableArray.shape $ initial hmm)
 
 alphaBeta ::
-   (Distr.EmissionProb distr, Distr.StateShape distr ~ sh, Eq sh,
-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,
+   (Distr.EmissionProb typ, Shape.C sh, Eq sh,
+    Class.Real prob, Distr.Emission typ prob ~ emission,
     Traversable f, NonEmptyC.Zip f, NonEmptyC.Reverse f) =>
-   T distr sh prob ->
+   T typ sh prob ->
    NonEmpty.T f emission ->
    (NonEmpty.T f (prob, Vector sh prob), NonEmpty.T f (Vector sh prob))
 alphaBeta hmm xs =
@@ -86,19 +83,19 @@
 
 
 xiFromAlphaBeta ::
-   (Distr.EmissionProb distr, Distr.StateShape distr ~ sh, Eq sh,
-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,
+   (Distr.EmissionProb typ, Shape.C sh, Eq sh,
+    Class.Real prob, Distr.Emission typ prob ~ emission,
     Traversable f, NonEmptyC.Zip f) =>
-   T distr sh prob ->
+   T typ sh prob ->
    NonEmpty.T f emission ->
    NonEmpty.T f (prob, Vector sh prob) ->
    NonEmpty.T f (Vector sh prob) ->
-   f (SquareMatrix sh prob)
+   f (Matrix.Square sh prob)
 xiFromAlphaBeta hmm xs calphas betas =
    let (cs,alphas) = Functor.unzip calphas
    in  NonEmptyC.zipWith4
           (\x alpha0 c1 beta1 ->
-             Vector.scale (recip c1) $ biscaleTransition hmm x alpha0 beta1)
+             Matrix.scale (recip c1) $ biscaleTransition hmm x alpha0 beta1)
           (NonEmpty.tail xs)
           (NonEmpty.init alphas)
           (NonEmpty.tail cs)
@@ -119,22 +116,10 @@
 It is found using the Viterbi algorithm.
 -}
 reveal ::
-   (Distr.EmissionProb distr, Distr.StateShape distr ~ sh,
-    Shape.InvIndexed sh, Eq sh, Shape.Index sh ~ state,
-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,
-    Traversable f, NonEmptyC.Reverse f) =>
-   T distr sh prob -> NonEmpty.T f emission -> NonEmpty.T f state
-reveal hmm (NonEmpty.Cons x xs) =
-   fmap (Shape.revealIndex (StorableArray.shape $ initial hmm)) $
-   uncurry (NonEmpty.scanr (StorableArray.!)) $
-   mapFst
-      (fst . Vector.argAbsMaximum .
-       StorableArray.mapShape Shape.Deferred) $
-   mapAccumL
-      (\alphai xi ->
-         swap $ mapSnd (Vector.mul (emission hmm xi)) $
-         matrixMaxMul (transition hmm) $ normalizeProb alphai)
-      (Vector.mul (emission hmm x) (initial hmm)) xs
+   (Distr.EmissionProb typ, Shape.InvIndexed sh, Eq sh, Shape.Index sh ~ state,
+    Distr.Emission typ prob ~ emission, Class.Real prob, Traversable f) =>
+   T typ sh prob -> NonEmpty.T f emission -> NonEmpty.T f state
+reveal = revealGen normalizeProb
 
 
 {- |
@@ -155,9 +140,9 @@
 This is done by the Baum-Welch algorithm.
 -}
 trainUnsupervised ::
-   (Distr.Estimate tdistr distr, Distr.StateShape distr ~ sh, Eq sh,
-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>
-   T distr sh prob -> NonEmpty.T [] emission -> Trained tdistr sh prob
+   (Distr.Estimate typ, Shape.C sh, Eq sh,
+    Class.Real prob, Distr.Emission typ prob ~ emission) =>
+   T typ sh prob -> NonEmpty.T [] emission -> Trained typ sh prob
 trainUnsupervised hmm xs =
    let (alphas, betas) = alphaBeta hmm xs
        zetas = zetaFromAlphaBeta alphas betas
@@ -168,8 +153,5 @@
           trainedTransition =
              sumTransitions hmm $ xiFromAlphaBeta hmm xs alphas betas,
           trainedDistribution =
-             Distr.accumulateEmissions $
-             Array.fromList (StorableArray.shape zeta0) $
-             map (zip (NonEmpty.flatten xs)) $
-             List.transpose $ map Vector.toList $ NonEmpty.flatten zetas
+             Distr.accumulateEmissionVectors $ NonEmptyC.zip xs zetas
        }
diff --git a/src/Math/HiddenMarkovModel/Pattern.hs b/src/Math/HiddenMarkovModel/Pattern.hs
--- a/src/Math/HiddenMarkovModel/Pattern.hs
+++ b/src/Math/HiddenMarkovModel/Pattern.hs
@@ -30,8 +30,10 @@
 import qualified Math.HiddenMarkovModel.Distribution as Distr
 import qualified Math.HiddenMarkovModel as HMM
 import Math.HiddenMarkovModel.Private (Trained(..))
-import Math.HiddenMarkovModel.Utility (SquareMatrix, squareConstant)
+import Math.HiddenMarkovModel.Utility (squareConstant)
 
+import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix
+import qualified Numeric.LAPACK.Matrix as Matrix
 import qualified Numeric.LAPACK.Vector as Vector
 import qualified Numeric.LAPACK.ShapeStatic as ShapeStatic
 
@@ -44,16 +46,16 @@
 import Data.FixedLength ((!:))
 
 import qualified Type.Data.Num.Unary.Literal as TypeNum
-import Type.Base.Proxy (Proxy(Proxy))
 
-import qualified Data.Map as Map
+import qualified Data.NonEmpty.Map as NonEmptyMap
+import qualified Data.NonEmpty as NonEmpty
 import Data.Semigroup (Semigroup, (<>), stimes)
 
 import Prelude hiding (replicate)
 
 
 newtype T sh prob =
-   Cons (sh -> (Shape.Index sh, SquareMatrix sh prob, Shape.Index sh))
+   Cons (sh -> (Shape.Index sh, Matrix.Square sh prob, Shape.Index sh))
 
 atom ::
    (Shape.Indexed sh, Shape.Index sh ~ state, Class.Real prob) =>
@@ -77,7 +79,7 @@
    Cons $ \n ->
       case (f n, g n) of
          ((sai, ma, sao), (sbi, mb, sbo)) ->
-            (sai, increment (sbi,sao) 1 $ Vector.add ma mb, sbo)
+            (sai, increment (sbi,sao) 1 $ Matrix.add ma mb, sbo)
 
 replicate ::
    (Shape.Indexed sh, Class.Real prob) => Int -> T sh prob -> T sh prob
@@ -86,25 +88,26 @@
       case f sh of
          (si, m, so) ->
             let k = fromIntegral ki
-            in  (si, increment (si,so) (k-1) $ Vector.scale k m, so)
+            in  (si, increment (si,so) (k-1) $ Matrix.scale k m, so)
 
 increment ::
    (Shape.Indexed sh, Shape.Index sh ~ state, Class.Real a) =>
-   (state, state) -> a -> SquareMatrix sh a -> SquareMatrix sh a
-increment (i,j) x m  =  StorableArray.accumulate (+) m [((i,j), x)]
+   (state, state) -> a -> Matrix.Square sh a -> Matrix.Square sh a
+increment (i,j) x =
+   ArrMatrix.lift1 $ flip (StorableArray.accumulate (+)) [((i,j), x)]
 
 
 finish ::
-   (Shape.Indexed sh, Class.Real prob) =>
-   sh -> tdistr -> T sh prob -> Trained tdistr sh prob
-finish sh tdistr (Cons f) =
-   case f sh of
-      (si, m, _so) ->
-         Trained {
-            trainedInitial = StorableArray.fromAssociations sh 0 [(si,1)],
-            trainedTransition = m,
-            trainedDistribution = tdistr
-         }
+   (Distr.Info typ, Shape.Indexed sh, Class.Real prob) =>
+   Distr.Trained typ sh prob -> T sh prob -> Trained typ sh prob
+finish tdistr (Cons f) =
+   let sh = Distr.statesShapeTrained tdistr
+       (si, m, _so) = f sh
+   in Trained {
+         trainedInitial = Vector.unit sh si,
+         trainedTransition = m,
+         trainedDistribution = tdistr
+      }
 
 
 _example :: HMM.DiscreteTrained Char (ShapeStatic.ZeroBased TypeNum.U2) Double
@@ -112,10 +115,9 @@
    let a = atom FL.i0
        b = atom FL.i1
        distr =
-          Distr.DiscreteTrained $ Map.fromList $
-          ('a', ShapeStatic.vector $ 1!:2!:FL.end) :
+          Distr.DiscreteTrained $ NonEmptyMap.fromList $
+          ('a', ShapeStatic.vector $ 1!:2!:FL.end) NonEmpty.!:
           ('b', ShapeStatic.vector $ 4!:3!:FL.end) :
           ('c', ShapeStatic.vector $ 0!:1!:FL.end) :
           []
-   in finish (ShapeStatic.ZeroBased Proxy) distr $
-      replicate 5 $ replicate 10 a <> replicate 20 b
+   in finish distr $ replicate 5 $ replicate 10 a <> replicate 20 b
diff --git a/src/Math/HiddenMarkovModel/Private.hs b/src/Math/HiddenMarkovModel/Private.hs
--- a/src/Math/HiddenMarkovModel/Private.hs
+++ b/src/Math/HiddenMarkovModel/Private.hs
@@ -1,15 +1,16 @@
 {-# LANGUAGE TypeFamilies #-}
-{-# LANGUAGE UndecidableInstances #-}
 module Math.HiddenMarkovModel.Private where
 
 import qualified Math.HiddenMarkovModel.Distribution as Distr
 import qualified Math.HiddenMarkovModel.CSV as HMMCSV
-import Math.HiddenMarkovModel.Utility (SquareMatrix, diagonal)
+import Math.HiddenMarkovModel.Utility (diagonal)
 
+import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix
+import qualified Numeric.LAPACK.Matrix.Square as Square
 import qualified Numeric.LAPACK.Matrix as Matrix
 import qualified Numeric.LAPACK.Vector as Vector
 import qualified Numeric.LAPACK.Format as Format
-import Numeric.LAPACK.Matrix ((<#), (<#>), (#>))
+import Numeric.LAPACK.Matrix ((-*#), (##*#), (#*##), (#*|))
 import Numeric.LAPACK.Vector (Vector)
 
 import qualified Numeric.Netlib.Class as Class
@@ -20,15 +21,15 @@
 import Foreign.Storable (Storable)
 
 import qualified Data.Array.Comfort.Storable as StorableArray
-import qualified Data.Array.Comfort.Boxed as Array
 import qualified Data.Array.Comfort.Shape as Shape
 
 import qualified Data.NonEmpty.Class as NonEmptyC
 import qualified Data.NonEmpty as NonEmpty
 import qualified Data.Semigroup as Sg
 import qualified Data.List as List
+import Data.Semigroup ((<>))
 import Data.Traversable (Traversable, mapAccumL)
-import Data.Tuple.HT (mapPair, mapFst, mapSnd, swap)
+import Data.Tuple.HT (mapFst, mapSnd, swap)
 
 
 {- |
@@ -46,88 +47,84 @@
 with respect to popular HMM descriptions.
 But I think this is the natural orientation, because this way
 you can write \"transition matrix times probability column vector\".
-
-The type has two type parameters,
-although the one for the distribution would be enough.
-However, replacing @prob@ by @Distr.Probability distr@
-would prohibit the derived Show and Read instances.
 -}
-data T distr sh prob =
+data T typ sh prob =
    Cons {
       initial :: Vector sh prob,
-      transition :: SquareMatrix sh prob,
-      distribution :: distr
+      transition :: Matrix.Square sh prob,
+      distribution :: Distr.T typ sh prob
    }
    deriving (Show)
 
 instance
-   (NFData distr, NFData sh, NFData prob, Storable prob) =>
-      NFData (T distr sh prob) where
+   (Distr.NFData typ, NFData sh, Shape.C sh, NFData prob, Storable prob) =>
+      NFData (T typ sh prob) where
    rnf (Cons initial_ transition_ distribution_) =
       rnf (initial_, transition_, distribution_)
 
 instance
-   (Class.Real prob, Format.FormatArray sh, Format.Format distr) =>
-      Format.Format (T distr sh prob) where
+   (Distr.Format typ, Format.FormatArray sh, Class.Real prob) =>
+      Format.Format (T typ sh prob) where
    format fmt (Cons initial_ transition_ distribution_) =
       Format.format fmt (initial_, transition_, distribution_)
 
+mapStatesShape ::
+   (Distr.EmissionProb typ, Shape.C sh0, Shape.C sh1) =>
+   (sh0 -> sh1) -> T typ sh0 prob -> T typ sh1 prob
+mapStatesShape f hmm =
+   Cons {
+      initial = StorableArray.mapShape f $ initial hmm,
+      transition = Square.mapSize f $ transition hmm,
+      distribution = Distr.mapStatesShape f $ distribution hmm
+   }
 
+
 emission ::
-   (Shape.C sh, Eq sh, sh ~ Distr.StateShape distr, Distr.EmissionProb distr,
-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>
-   T distr sh prob -> emission -> Vector sh prob
+   (Distr.EmissionProb typ, Shape.C sh, Eq sh, Class.Real prob) =>
+   T typ sh prob -> Distr.Emission typ prob -> Vector sh prob
 emission  =  Distr.emissionProb . distribution
 
 
 forward ::
-   (Shape.C sh, Eq sh, sh ~ Distr.StateShape distr, Distr.EmissionProb distr,
-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,
-    Traversable f) =>
-   T distr sh prob -> NonEmpty.T f emission -> prob
+   (Distr.EmissionProb typ, Shape.C sh, Eq sh, Class.Real prob,
+    Distr.Emission typ prob ~ emission, Traversable f) =>
+   T typ sh prob -> NonEmpty.T f emission -> prob
 forward hmm = Vector.sum . NonEmpty.last . alpha hmm
 
 alpha ::
-   (Shape.C sh, Eq sh, sh ~ Distr.StateShape distr, Distr.EmissionProb distr,
-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,
-    Traversable f) =>
-   T distr sh prob ->
-   NonEmpty.T f emission -> NonEmpty.T f (Vector sh prob)
+   (Distr.EmissionProb typ, Shape.C sh, Eq sh, Class.Real prob,
+    Distr.Emission typ prob ~ emission, Traversable f) =>
+   T typ sh prob -> NonEmpty.T f emission -> NonEmpty.T f (Vector sh prob)
 alpha hmm (NonEmpty.Cons x xs) =
    NonEmpty.scanl
-      (\alphai xi -> Vector.mul (emission hmm xi) (transition hmm #> alphai))
+      (\alphai xi -> Vector.mul (emission hmm xi) (transition hmm #*| alphai))
       (Vector.mul (emission hmm x) (initial hmm))
       xs
 
 
 backward ::
-   (Shape.C sh, Eq sh, sh ~ Distr.StateShape distr, Distr.EmissionProb distr,
-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,
-    Traversable f) =>
-   T distr sh prob -> NonEmpty.T f emission -> prob
+   (Distr.EmissionProb typ, Shape.C sh, Eq sh, Class.Real prob,
+    Distr.Emission typ prob ~ emission, Traversable f) =>
+   T typ sh prob -> NonEmpty.T f emission -> prob
 backward hmm (NonEmpty.Cons x xs) =
-   Vector.sum $
-   Vector.mul (initial hmm) $
+   Vector.dot (initial hmm) $
    Vector.mul (emission hmm x) $
    NonEmpty.head $ beta hmm xs
 
 beta ::
-   (Shape.C sh, Eq sh, sh ~ Distr.StateShape distr, Distr.EmissionProb distr,
-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,
-    Traversable f) =>
-   T distr sh prob ->
-   f emission -> NonEmpty.T f (Vector sh prob)
+   (Distr.EmissionProb typ, Shape.C sh, Eq sh, Class.Real prob,
+    Distr.Emission typ prob ~ emission, Traversable f) =>
+   T typ sh prob -> f emission -> NonEmpty.T f (Vector sh prob)
 beta hmm =
    NonEmpty.scanr
-      (\xi betai -> Vector.mul (emission hmm xi) betai <# transition hmm)
-      (Vector.constant (StorableArray.shape $ initial hmm) 1)
+      (\xi betai -> Vector.mul (emission hmm xi) betai -*# transition hmm)
+      (Vector.one $ StorableArray.shape $ initial hmm)
 
 
 alphaBeta ::
-   (Shape.C sh, Eq sh, sh ~ Distr.StateShape distr, Distr.EmissionProb distr,
-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,
-    Traversable f) =>
-   T distr sh prob ->
+   (Distr.EmissionProb typ, Shape.C sh, Eq sh, Class.Real prob,
+    Distr.Emission typ prob ~ emission, Traversable f) =>
+   T typ sh prob ->
    NonEmpty.T f emission ->
    (prob, NonEmpty.T f (Vector sh prob), NonEmpty.T f (Vector sh prob))
 alphaBeta hmm xs =
@@ -139,29 +136,28 @@
 
 
 biscaleTransition ::
-   (Shape.C sh, Eq sh, sh ~ Distr.StateShape distr,
-    Distr.EmissionProb distr, Distr.Probability distr ~ prob) =>
-   T distr sh prob -> Distr.Emission distr ->
-   Vector sh prob -> Vector sh prob -> SquareMatrix sh prob
+   (Distr.EmissionProb typ, Shape.C sh, Eq sh, Class.Real prob) =>
+   T typ sh prob -> Distr.Emission typ prob ->
+   Vector sh prob -> Vector sh prob -> Matrix.Square sh prob
 biscaleTransition hmm x alpha0 beta1 =
-   diagonal (Vector.mul (emission hmm x) beta1)
-   <#>
-   transition hmm
-   <#>
+   (diagonal (Vector.mul (emission hmm x) beta1)
+    #*##
+    transition hmm)
+   ##*#
    diagonal alpha0
 
 xiFromAlphaBeta ::
-   (Shape.C sh, Eq sh, sh ~ Distr.StateShape distr, Distr.EmissionProb distr,
-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>
-   T distr sh prob -> prob ->
+   (Distr.EmissionProb typ, Shape.C sh, Eq sh, Class.Real prob,
+    Distr.Emission typ prob ~ emission) =>
+   T typ sh prob -> prob ->
    NonEmpty.T [] emission ->
    NonEmpty.T [] (Vector sh prob) ->
    NonEmpty.T [] (Vector sh prob) ->
-   [SquareMatrix sh prob]
+   [Matrix.Square sh prob]
 xiFromAlphaBeta hmm recipLikelihood xs alphas betas =
    zipWith3
       (\x alpha0 beta1 ->
-         Vector.scale recipLikelihood $
+         Matrix.scale recipLikelihood $
          biscaleTransition hmm x alpha0 beta1)
       (NonEmpty.tail xs)
       (NonEmpty.init alphas)
@@ -169,7 +165,7 @@
 
 zetaFromXi ::
    (Shape.C sh, Eq sh, Class.Real prob) =>
-   [SquareMatrix sh prob] -> [Vector sh prob]
+   [Matrix.Square sh prob] -> [Vector sh prob]
 zetaFromXi = map Matrix.columnSums
 
 zetaFromAlphaBeta ::
@@ -191,35 +187,41 @@
 than the smallest representable number.
 -}
 reveal ::
-   (Shape.InvIndexed sh, Shape.Index sh ~ state,
-    Eq sh, sh ~ Distr.StateShape distr, Distr.EmissionProb distr,
-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission,
-    Traversable f) =>
-   T distr sh prob -> NonEmpty.T f emission -> NonEmpty.T f state
-reveal hmm (NonEmpty.Cons x xs) =
-   fmap (Shape.revealIndex (StorableArray.shape $ initial hmm)) $
+   (Distr.EmissionProb typ, Shape.InvIndexed sh, Eq sh, Shape.Index sh ~ state,
+    Distr.Emission typ prob ~ emission, Class.Real prob, Traversable f) =>
+   T typ sh prob -> NonEmpty.T f emission -> NonEmpty.T f state
+reveal = revealGen id
+
+revealGen ::
+   (Distr.EmissionProb typ, Shape.InvIndexed sh, Eq sh, Shape.Index sh ~ state,
+    Distr.Emission typ prob ~ emission, Class.Real prob, Traversable f) =>
+   (Vector (Shape.Deferred sh) prob -> Vector (Shape.Deferred sh) prob) ->
+   T typ sh prob -> NonEmpty.T f emission -> NonEmpty.T f state
+revealGen normalize hmm =
+   fmap (Shape.revealIndex (StorableArray.shape $ initial hmm)) .
+   revealStorable normalize (mapStatesShape Shape.Deferred hmm)
+
+revealStorable ::
+   (Distr.EmissionProb typ, Shape.InvIndexed sh, Eq sh,
+    Shape.Index sh ~ state, Storable state,
+    Distr.Emission typ prob ~ emission, Class.Real prob, Traversable f) =>
+   (Vector sh prob -> Vector sh prob) ->
+   T typ sh prob -> NonEmpty.T f emission -> NonEmpty.T f state
+revealStorable normalize hmm (NonEmpty.Cons x xs) =
    uncurry (NonEmpty.scanr (StorableArray.!)) $
-   mapFst
-      (fst . Vector.argAbsMaximum .
-       StorableArray.mapShape Shape.Deferred) $
+   mapFst (fst . Vector.argAbsMaximum) $
    mapAccumL
       (\alphai xi ->
          swap $ mapSnd (Vector.mul (emission hmm xi)) $
-         matrixMaxMul (transition hmm) alphai)
+         matrixMaxMul (transition hmm) $ normalize alphai)
       (Vector.mul (emission hmm x) (initial hmm)) xs
 
 matrixMaxMul ::
-   (Shape.Indexed sh, Eq sh, Shape.Index sh ~ ix, Class.Real a) =>
-   SquareMatrix sh a -> Vector sh a ->
-   (Vector (Shape.Deferred sh) (Shape.DeferredIndex ix), Vector sh a)
-matrixMaxMul m v =
-   let sh = StorableArray.shape v
-   in mapPair (Vector.fromList (Shape.Deferred sh), Vector.fromList sh) $
-      unzip $
-      map (Vector.argAbsMaximum .
-           StorableArray.mapShape Shape.Deferred .
-           Vector.mul v) $
-      Matrix.toRows m
+   (Shape.InvIndexed sh, Eq sh, Shape.Index sh ~ ix, Storable ix,
+    Class.Real a) =>
+   Matrix.Square sh a -> Vector sh a ->
+   (Vector sh ix, Vector sh a)
+matrixMaxMul m v = Matrix.rowArgAbsMaximums $ Matrix.scaleColumns v m
 
 
 
@@ -239,35 +241,35 @@
 
 * derive it from state sequence patterns, cf. "Math.HiddenMarkovModel.Pattern".
 -}
-data Trained distr sh prob =
+data Trained typ sh prob =
    Trained {
       trainedInitial :: Vector sh prob,
-      trainedTransition :: SquareMatrix sh prob,
-      trainedDistribution :: distr
+      trainedTransition :: Matrix.Square sh prob,
+      trainedDistribution :: Distr.Trained typ sh prob
    }
    deriving (Show)
 
 instance
-   (NFData distr, NFData sh, NFData prob, Storable prob) =>
-      NFData (Trained distr sh prob) where
+   (Distr.NFData typ, NFData sh, Shape.C sh, NFData prob, Storable prob) =>
+      NFData (Trained typ sh prob) where
    rnf hmm =
       rnf (trainedInitial hmm, trainedTransition hmm, trainedDistribution hmm)
 
 
 sumTransitions ::
    (Shape.C sh, Eq sh, Class.Real e) =>
-   T distr sh e -> [SquareMatrix sh e] -> SquareMatrix sh e
+   T typ sh e -> [Matrix.Square sh e] -> Matrix.Square sh e
 sumTransitions hmm =
-   List.foldl' Vector.add
-      (Vector.constant (StorableArray.shape $ transition hmm) 0)
+   List.foldl' Matrix.add $
+   Matrix.zero $ ArrMatrix.shape $ transition hmm
 
 {- |
 Baum-Welch algorithm
 -}
 trainUnsupervised ::
-   (Distr.Estimate tdistr distr, Distr.StateShape distr ~ sh, Eq sh,
-    Distr.Probability distr ~ prob, Distr.Emission distr ~ emission) =>
-   T distr sh prob -> NonEmpty.T [] emission -> Trained tdistr sh prob
+   (Distr.Estimate typ, Shape.C sh, Eq sh, Class.Real prob,
+    Distr.Emission typ prob ~ emission) =>
+   T typ sh prob -> NonEmpty.T [] emission -> Trained typ sh prob
 trainUnsupervised hmm xs =
    let (recipLikelihood, alphas, betas) = alphaBeta hmm xs
        zetas = zetaFromAlphaBeta recipLikelihood alphas betas
@@ -279,37 +281,31 @@
              sumTransitions hmm $
              xiFromAlphaBeta hmm recipLikelihood xs alphas betas,
           trainedDistribution =
-             Distr.accumulateEmissions $
-             Array.fromList (StorableArray.shape zeta0) $
-             map (zip (NonEmpty.flatten xs)) $
-             List.transpose $ map Vector.toList $ NonEmpty.flatten zetas
+             Distr.accumulateEmissionVectors $ NonEmptyC.zip xs zetas
        }
 
 
 mergeTrained ::
-   (Shape.C sh, Eq sh,
-    Distr.Estimate tdistr distr, Distr.Probability distr ~ prob) =>
-   Trained tdistr sh prob -> Trained tdistr sh prob -> Trained tdistr sh prob
+   (Distr.Estimate typ, Shape.C sh, Eq sh, Class.Real prob) =>
+   Trained typ sh prob -> Trained typ sh prob -> Trained typ sh prob
 mergeTrained hmm0 hmm1 =
    Trained {
       trainedInitial = Vector.add (trainedInitial hmm0) (trainedInitial hmm1),
       trainedTransition =
-         Vector.add (trainedTransition hmm0) (trainedTransition hmm1),
+         Matrix.add (trainedTransition hmm0) (trainedTransition hmm1),
       trainedDistribution =
-         Distr.combine
-            (trainedDistribution hmm0) (trainedDistribution hmm1)
+         trainedDistribution hmm0 <> trainedDistribution hmm1
    }
 
 instance
-   (Shape.C sh, Eq sh,
-    Distr.Estimate tdistr distr, Distr.Probability distr ~ prob) =>
-      Sg.Semigroup (Trained tdistr sh prob) where
+   (Distr.Estimate typ, Shape.C sh, Eq sh, Class.Real prob) =>
+      Sg.Semigroup (Trained typ sh prob) where
    (<>) = mergeTrained
 
 
 toCells ::
-   (Distr.ToCSV distr, Shape.Indexed sh, Class.Real prob, Show prob) =>
-   T distr sh prob -> [[String]]
+   (Distr.ToCSV typ, Shape.Indexed sh, Class.Real prob, Show prob) =>
+   T typ sh prob -> [[String]]
 toCells hmm =
    (HMMCSV.cellsFromVector $ initial hmm) :
    (HMMCSV.cellsFromSquare $ transition hmm) ++
@@ -317,9 +313,9 @@
    (Distr.toCells $ distribution hmm)
 
 parseCSV ::
-   (Distr.FromCSV distr, Distr.StateShape distr ~ stateSh, Shape.C stateSh,
+   (Distr.FromCSV typ, Shape.C stateSh, Eq stateSh,
     Class.Real prob, Read prob) =>
-   (Int -> stateSh) -> HMMCSV.CSVParser (T distr stateSh prob)
+   (Int -> stateSh) -> HMMCSV.CSVParser (T typ stateSh prob)
 parseCSV makeShape = do
    v <-
       StorableArray.mapShape (makeShape . Shape.zeroBasedSize) <$>
diff --git a/src/Math/HiddenMarkovModel/Test.hs b/src/Math/HiddenMarkovModel/Test.hs
--- a/src/Math/HiddenMarkovModel/Test.hs
+++ b/src/Math/HiddenMarkovModel/Test.hs
@@ -11,8 +11,10 @@
 import qualified Math.HiddenMarkovModel.Normalized as Normalized
 import qualified Math.HiddenMarkovModel.Private as Priv
 import qualified Math.HiddenMarkovModel.Distribution as Distr
-import Math.HiddenMarkovModel.Utility (SquareMatrix, squareFromLists, distance)
+import Math.HiddenMarkovModel.Utility
+         (squareFromLists, distance, matrixDistance)
 
+import qualified Numeric.LAPACK.Matrix as Matrix
 import qualified Numeric.LAPACK.Vector as Vector
 import qualified Numeric.LAPACK.ShapeStatic as ShapeStatic
 import Numeric.LAPACK.Vector (Vector)
@@ -30,10 +32,12 @@
 import Control.DeepSeq (deepseq)
 
 import qualified Data.NonEmpty.Class as NonEmptyC
+import qualified Data.NonEmpty.Map as NonEmptyMap
 import qualified Data.NonEmpty as NonEmpty
 import qualified Data.Traversable as Trav
 import qualified Data.Foldable as Fold
 import qualified Data.Map as Map
+import Data.NonEmpty ((!:))
 import Data.Tuple.HT (mapSnd)
 
 import Text.Printf (printf)
@@ -53,8 +57,8 @@
             stateVector 0.1 0.1 0.2 0.8 :
             [],
       HMM.distribution =
-         Distr.Discrete $ Map.fromList $
-            ('a', stateVector 1 0 0 0) :
+         Distr.discreteFromList $
+            ('a', stateVector 1 0 0 0) !:
             ('b', stateVector 0 1 0 1) :
             ('c', stateVector 0 0 1 0) :
             []
@@ -83,8 +87,7 @@
             1 -> True
             _ -> error "invalid emission probability (must be 0 or 1)") $
    Vector.toList $
-   Map.findWithDefault (error "invalid character") c $
-   case HMM.distribution hmm of Distr.Discrete m -> m
+   case HMM.distribution hmm of Distr.Discrete m -> Matrix.takeRow m c
 
 {- |
 Should all be equal.
@@ -140,7 +143,7 @@
 {- |
 Lists should be equal
 -}
-xis :: ([SquareMatrix StateSet Double], [SquareMatrix StateSet Double])
+xis :: ([Matrix.Square StateSet Double], [Matrix.Square StateSet Double])
 xis =
    let (recipLikelihood, alphas, betas) = Priv.alphaBeta hmm sequ
    in  (Priv.xiFromAlphaBeta hmm recipLikelihood sequ alphas betas,
@@ -153,7 +156,8 @@
 xisDiff :: (Bool, Double)
 xisDiff =
    case xis of
-      (x0,x1) -> (length x0 == length x1, maximum $ zipWith distance x0 x1)
+      (x0,x1) ->
+         (length x0 == length x1, maximum $ zipWith matrixDistance x0 x1)
 
 
 reveal :: Bool
@@ -172,13 +176,16 @@
 trainUnsupervisedDiff =
    case trainUnsupervised of
       (hmm0,hmm1) ->
-         (distance (Priv.trainedTransition hmm0) (Priv.trainedTransition hmm1),
+         (matrixDistance
+             (Priv.trainedTransition hmm0) (Priv.trainedTransition hmm1),
           distance
              (Priv.trainedInitial hmm0) (Priv.trainedInitial hmm1),
           case (Priv.trainedDistribution hmm0, Priv.trainedDistribution hmm1) of
              (Distr.DiscreteTrained m0, Distr.DiscreteTrained m1) ->
-                (Map.size m0 == Map.size m1,
-                 Fold.maximum $ Map.intersectionWith distance m0 m1))
+                (NonEmptyMap.size m0 == NonEmptyMap.size m1,
+                 Fold.maximum $
+                 Map.intersectionWith distance
+                    (NonEmptyMap.flatten m0) (NonEmptyMap.flatten m1)))
 
 
 nonEmptyScanr :: Int -> [Int] -> Bool
diff --git a/src/Math/HiddenMarkovModel/Utility.hs b/src/Math/HiddenMarkovModel/Utility.hs
--- a/src/Math/HiddenMarkovModel/Utility.hs
+++ b/src/Math/HiddenMarkovModel/Utility.hs
@@ -4,10 +4,12 @@
 import qualified Numeric.LAPACK.Matrix.Hermitian as Hermitian
 import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape
 import qualified Numeric.LAPACK.Matrix.Square as Square
+import qualified Numeric.LAPACK.Matrix.Array as ArrMatrix
 import qualified Numeric.LAPACK.Matrix as Matrix
 import qualified Numeric.LAPACK.Vector as Vector
 import Numeric.LAPACK.Matrix.Triangular (Diagonal)
-import Numeric.LAPACK.Vector (Vector)
+import Numeric.LAPACK.Matrix.Array (ArrayMatrix)
+import Numeric.LAPACK.Vector (Vector, (.*|))
 
 import qualified Numeric.Netlib.Class as Class
 
@@ -22,15 +24,13 @@
 import qualified Control.Monad.Trans.State as MS
 
 
-type SquareMatrix sh = Square.Square sh
-
 normalizeProb :: (Shape.C sh, Class.Real a) => Vector sh a -> Vector sh a
 normalizeProb = snd . normalizeFactor
 
 normalizeFactor :: (Shape.C sh, Class.Real a) => Vector sh a -> (a, Vector sh a)
 normalizeFactor xs =
    let c = Vector.sum xs
-   in  (c, Vector.scale (recip c) xs)
+   in  (c, recip c .*| xs)
 
 -- see htam:Stochastic
 randomItemProp ::
@@ -57,11 +57,13 @@
 
 
 squareConstant ::
-   (Shape.C sh, Class.Real a) => sh -> a -> SquareMatrix sh a
-squareConstant = Vector.constant . MatrixShape.square MatrixShape.RowMajor
+   (Shape.C sh, Class.Real a) => sh -> a -> Matrix.Square sh a
+squareConstant =
+   (ArrMatrix.fromVector .) .
+   Vector.constant . MatrixShape.square MatrixShape.RowMajor
 
 squareFromLists ::
-   (Shape.C sh, Eq sh, Storable a) => sh -> [Vector sh a] -> SquareMatrix sh a
+   (Shape.C sh, Eq sh, Storable a) => sh -> [Vector sh a] -> Matrix.Square sh a
 squareFromLists sh =
    Square.fromGeneral . Matrix.fromRowArray sh . Array.fromList sh
 
@@ -79,3 +81,8 @@
    Class.switchReal
       (Distance $ (Vector.normInf .) . Vector.sub)
       (Distance $ (Vector.normInf .) . Vector.sub)
+
+matrixDistance ::
+   (Shape.C sh, Eq sh, Class.Real a) =>
+   ArrayMatrix sh a -> ArrayMatrix sh a -> a
+matrixDistance a b = distance (ArrMatrix.toVector a) (ArrMatrix.toVector b)