diff --git a/ChangeLog.md b/ChangeLog.md
--- a/ChangeLog.md
+++ b/ChangeLog.md
@@ -2,7 +2,14 @@
 # Changelog for ELynx
 
 
-## Version 0.3.1; in development
+## Version 0.3.2; in development
+
+
+## Version 0.3.1
+
+-   Use Attoparsec.
+-   Use ByteString consistently.
+-   Remove elynx-tools dependency from libaries.
 
 
 ## Version 0.3.0
diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -2,7 +2,7 @@
 
 # The ELynx Suite
 
-Version: 0.3.0.
+Version: 0.3.1.
 Reproducible evolution made easy.
 
 <p align="center"><img src="https://travis-ci.org/dschrempf/elynx.svg?branch=master"/></p>
@@ -90,9 +90,9 @@
 
     slynx --help | head -n -16
 
-    ELynx Suite version 0.3.0.
+    ELynx Suite version 0.3.1.
     Developed by Dominik Schrempf.
-    Compiled on July 30, 2020, at 12:40 pm, UTC.
+    Compiled on August 3, 2020, at 12:12 pm, UTC.
     
     Usage: slynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME] 
                  [-f|--force] COMMAND
@@ -136,9 +136,9 @@
 
     slynx concatenate --help
 
-    ELynx Suite version 0.3.0.
+    ELynx Suite version 0.3.1.
     Developed by Dominik Schrempf.
-    Compiled on July 30, 2020, at 12:40 pm, UTC.
+    Compiled on August 3, 2020, at 12:12 pm, UTC.
     
     Usage: slynx concatenate (-a|--alphabet NAME) INPUT-FILE
       Concatenate sequences found in input files.
@@ -157,9 +157,9 @@
 
     slynx examine --help
 
-    ELynx Suite version 0.3.0.
+    ELynx Suite version 0.3.1.
     Developed by Dominik Schrempf.
-    Compiled on July 30, 2020, at 12:40 pm, UTC.
+    Compiled on August 3, 2020, at 12:12 pm, UTC.
     
     Usage: slynx examine (-a|--alphabet NAME) INPUT-FILE [--per-site]
       Examine sequences. If data is a multi sequence alignment, additionally analyze columns.
@@ -179,9 +179,9 @@
 
     slynx filter-rows --help
 
-    ELynx Suite version 0.3.0.
+    ELynx Suite version 0.3.1.
     Developed by Dominik Schrempf.
-    Compiled on July 30, 2020, at 12:40 pm, UTC.
+    Compiled on August 3, 2020, at 12:12 pm, UTC.
     
     Usage: slynx filter-rows (-a|--alphabet NAME) INPUT-FILE [--longer-than LENGTH] 
                              [--shorter-than LENGTH] [--standard-characters]
@@ -202,9 +202,9 @@
 
     slynx filter-columns --help
 
-    ELynx Suite version 0.3.0.
+    ELynx Suite version 0.3.1.
     Developed by Dominik Schrempf.
-    Compiled on July 30, 2020, at 12:40 pm, UTC.
+    Compiled on August 3, 2020, at 12:12 pm, UTC.
     
     Usage: slynx filter-columns (-a|--alphabet NAME) INPUT-FILE 
                                 [--standard-chars DOUBLE]
@@ -226,9 +226,9 @@
 
     slynx simulate --help
 
-    ELynx Suite version 0.3.0.
+    ELynx Suite version 0.3.1.
     Developed by Dominik Schrempf.
-    Compiled on July 30, 2020, at 12:40 pm, UTC.
+    Compiled on August 3, 2020, at 12:12 pm, UTC.
     
     Usage: slynx simulate (-t|--tree-file Name) [-s|--substitution-model MODEL] 
                           [-m|--mixture-model MODEL] [-e|--edm-file NAME] 
@@ -305,9 +305,9 @@
 
     slynx sub-sample --help
 
-    ELynx Suite version 0.3.0.
+    ELynx Suite version 0.3.1.
     Developed by Dominik Schrempf.
-    Compiled on July 30, 2020, at 12:40 pm, UTC.
+    Compiled on August 3, 2020, at 12:12 pm, UTC.
     
     Usage: slynx sub-sample (-a|--alphabet NAME) INPUT-FILE
                             (-n|--number-of-sites INT)
@@ -335,9 +335,9 @@
 
     slynx translate --help
 
-    ELynx Suite version 0.3.0.
+    ELynx Suite version 0.3.1.
     Developed by Dominik Schrempf.
-    Compiled on July 30, 2020, at 12:40 pm, UTC.
+    Compiled on August 3, 2020, at 12:12 pm, UTC.
     
     Usage: slynx translate (-a|--alphabet NAME) INPUT-FILE (-r|--reading-frame INT)
                            (-u|--universal-code CODE)
@@ -360,9 +360,9 @@
 
     tlynx --help | head -n -16
 
-    ELynx Suite version 0.3.0.
+    ELynx Suite version 0.3.1.
     Developed by Dominik Schrempf.
-    Compiled on July 30, 2020, at 12:40 pm, UTC.
+    Compiled on August 3, 2020, at 12:12 pm, UTC.
     
     Usage: tlynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME] 
                  [-f|--force] COMMAND
@@ -399,9 +399,9 @@
 
     tlynx compare --help
 
-    ELynx Suite version 0.3.0.
+    ELynx Suite version 0.3.1.
     Developed by Dominik Schrempf.
-    Compiled on July 30, 2020, at 12:40 pm, UTC.
+    Compiled on August 3, 2020, at 12:12 pm, UTC.
     
     Usage: tlynx compare [-n|--normalize] [-b|--bipartitions] [-t|--intersect] 
                          [-f|--newick-format FORMAT] NAMES
@@ -428,9 +428,9 @@
 
     tlynx examine --help
 
-    ELynx Suite version 0.3.0.
+    ELynx Suite version 0.3.1.
     Developed by Dominik Schrempf.
-    Compiled on July 30, 2020, at 12:40 pm, UTC.
+    Compiled on August 3, 2020, at 12:12 pm, UTC.
     
     Usage: tlynx examine INPUT-FILE [-f|--newick-format FORMAT]
       Compute summary statistics of phylogenetic trees.
@@ -452,9 +452,9 @@
 
     tlynx simulate --help
 
-    ELynx Suite version 0.3.0.
+    ELynx Suite version 0.3.1.
     Developed by Dominik Schrempf.
-    Compiled on July 30, 2020, at 12:40 pm, UTC.
+    Compiled on August 3, 2020, at 12:12 pm, UTC.
     
     Usage: tlynx simulate [-t|--nTrees INT] [-n|--nLeaves INT] [-H|--height DOUBLE] 
                           [-M|--condition-on-mrca] [-l|--lambda DOUBLE] 
@@ -491,9 +491,9 @@
 
     elynx --help | head -n -16
 
-    ELynx Suite version 0.3.0.
+    ELynx Suite version 0.3.1.
     Developed by Dominik Schrempf.
-    Compiled on July 30, 2020, at 12:40 pm, UTC.
+    Compiled on August 3, 2020, at 12:12 pm, UTC.
     
     Usage: elynx COMMAND
       Validate and redo past ELynx analyses
diff --git a/elynx-markov.cabal b/elynx-markov.cabal
--- a/elynx-markov.cabal
+++ b/elynx-markov.cabal
@@ -4,10 +4,10 @@
 --
 -- see: https://github.com/sol/hpack
 --
--- hash: d6c9194049e2cfd023ea390aefa6b761f0297f5484539c466a7f5b646e156877
+-- hash: c05f62a676b8e0fa2f4741dbd0fc99542f112c2421694e8e1b3a5e59ab5f235b
 
 name:           elynx-markov
-version:        0.3.0
+version:        0.3.1
 synopsis:       Simulate molecular sequences along trees
 description:    Examine, modify, and simulate molecular sequences in a reproducible way. Please see the README on GitHub at <https://github.com/dschrempf/elynx>.
 category:       Bioinformatics
@@ -50,16 +50,15 @@
       src
   ghc-options: -Wall -fllvm
   build-depends:
-      base >=4.7 && <5
+      async
+    , attoparsec
+    , base >=4.7 && <5
     , bytestring
     , containers
     , elynx-seq
-    , elynx-tools
-    , elynx-tree
     , hmatrix
     , integration
     , math-functions
-    , megaparsec
     , mwc-random
     , parallel
     , primitive
@@ -86,10 +85,8 @@
     , containers
     , elynx-markov
     , elynx-tools
-    , elynx-tree
     , hmatrix
     , hspec
-    , hspec-megaparsec
     , mwc-random
     , vector
   default-language: Haskell2010
diff --git a/src/ELynx/Data/MarkovProcess/AminoAcid.hs b/src/ELynx/Data/MarkovProcess/AminoAcid.hs
--- a/src/ELynx/Data/MarkovProcess/AminoAcid.hs
+++ b/src/ELynx/Data/MarkovProcess/AminoAcid.hs
@@ -22,13 +22,14 @@
   )
 where
 
+import Data.ByteString.Internal (c2w)
 import Data.List (elemIndex)
 import Data.Maybe (fromMaybe)
+import qualified Data.Vector.Storable as V
 import Data.Word (Word8)
 import ELynx.Data.Alphabet.Alphabet
 import ELynx.Data.MarkovProcess.RateMatrix
 import ELynx.Data.MarkovProcess.SubstitutionModel
-import ELynx.Tools
 import Numeric.LinearAlgebra
 
 n :: Int
@@ -321,10 +322,16 @@
 lgExch :: ExchangeabilityMatrix
 lgExch = pamlToAlphaMat $ exchFromListLower n lgExchRawPaml
 
+normalizeSumVec :: Vector Double -> Vector Double
+normalizeSumVec v = V.map (/ s) v
+  where
+    s = V.sum v
+{-# INLINE normalizeSumVec #-}
+
 -- Stationary distribution in PAML order.
 lgStatDistPaml :: StationaryDistribution
 lgStatDistPaml =
-  normalizeSumVec 1.0 $
+  normalizeSumVec $
     fromList
       [ 0.079066,
         0.055941,
@@ -564,7 +571,7 @@
 -- WAG stationary distribution in PAML order.
 wagStatDistPaml :: StationaryDistribution
 wagStatDistPaml =
-  normalizeSumVec 1.0 $
+  normalizeSumVec $
     fromList
       [ 0.0866279,
         0.043972,
@@ -602,15 +609,22 @@
   where
     nm = fromMaybe "WAG-Custom" mnm
 
+matrixSetDiagToZero :: Matrix R -> Matrix R
+matrixSetDiagToZero m = m - diag (takeDiag m)
+{-# INLINE matrixSetDiagToZero #-}
+
 uniformExch :: ExchangeabilityMatrix
 uniformExch = matrixSetDiagToZero $ matrix n $ replicate (n * n) 1.0
 
 poissonExch :: ExchangeabilityMatrix
 poissonExch = uniformExch
 
+uniformVec :: Vector Double
+uniformVec = V.replicate n (1 / fromIntegral n)
+
 -- | Poisson substitution model.
 poisson :: SubstitutionModel
-poisson = substitutionModel Protein "Poisson" [] (uniformVec n) poissonExch
+poisson = substitutionModel Protein "Poisson" [] uniformVec poissonExch
 
 -- | Poisson substitution model with maybe a name and a custom stationary distribution.
 poissonCustom :: Maybe String -> StationaryDistribution -> SubstitutionModel
diff --git a/src/ELynx/Data/MarkovProcess/CXXModels.hs b/src/ELynx/Data/MarkovProcess/CXXModels.hs
--- a/src/ELynx/Data/MarkovProcess/CXXModels.hs
+++ b/src/ELynx/Data/MarkovProcess/CXXModels.hs
@@ -10,7 +10,7 @@
 --
 -- Creation date: Tue Feb 26 16:44:33 2019.
 --
--- Quang, L. S., Gascuel, O. & Lartillot, N. Empirical profile mixture models for
+-- Quang, BL. S., Gascuel, O. & Lartillot, N. Empirical profile mixture models for
 -- phylogenetic reconstruction. Bioinformatics 24, 2317–2323 (2008).
 --
 -- XXX: For now, I only provide Poisson exchangeabilities.
diff --git a/src/ELynx/Data/MarkovProcess/CXXModelsData.hs b/src/ELynx/Data/MarkovProcess/CXXModelsData.hs
--- a/src/ELynx/Data/MarkovProcess/CXXModelsData.hs
+++ b/src/ELynx/Data/MarkovProcess/CXXModelsData.hs
@@ -10,7 +10,7 @@
 --
 -- Creation date: Tue Feb 26 17:17:35 2019.
 --
--- Quang, L. S., Gascuel, O. & Lartillot, N. Empirical profile mixture models for
+-- Quang, BL. S., Gascuel, O. & Lartillot, N. Empirical profile mixture models for
 -- phylogenetic reconstruction. Bioinformatics 24, 2317–2323 (2008).
 module ELynx.Data.MarkovProcess.CXXModelsData
   ( c10StatDists,
diff --git a/src/ELynx/Data/MarkovProcess/GammaRateHeterogeneity.hs b/src/ELynx/Data/MarkovProcess/GammaRateHeterogeneity.hs
--- a/src/ELynx/Data/MarkovProcess/GammaRateHeterogeneity.hs
+++ b/src/ELynx/Data/MarkovProcess/GammaRateHeterogeneity.hs
@@ -20,7 +20,7 @@
   )
 where
 
-import qualified Data.ByteString.Lazy.Char8 as L
+import qualified Data.ByteString.Lazy.Char8 as BL
 import Data.List.NonEmpty (NonEmpty)
 import qualified Data.List.NonEmpty as N
 import qualified ELynx.Data.MarkovProcess.MixtureModel as M
@@ -32,10 +32,10 @@
 import Prelude hiding (repeat)
 
 -- | Short summary of gamma rate heterogeneity parameters.
-summarizeGammaRateHeterogeneity :: Int -> Double -> [L.ByteString]
+summarizeGammaRateHeterogeneity :: Int -> Double -> [BL.ByteString]
 summarizeGammaRateHeterogeneity n alpha =
   map
-    L.pack
+    BL.pack
     [ "Discrete gamma rate heterogeneity.",
       "Number of categories: " ++ show n,
       "Shape parameter of gamma distribution: " ++ show alpha,
@@ -106,9 +106,8 @@
     -- The last category has to be calculated with an improper integration.
     lastMean = integralAToInf meanFunc (quantiles !! (n - 1))
 
--- The error of integration.
 eps :: Double
-eps = 1e-6
+eps = 1e-8
 
 -- The integration method to use
 method :: (Double -> Double) -> Double -> Double -> [Result]
diff --git a/src/ELynx/Data/MarkovProcess/MixtureModel.hs b/src/ELynx/Data/MarkovProcess/MixtureModel.hs
--- a/src/ELynx/Data/MarkovProcess/MixtureModel.hs
+++ b/src/ELynx/Data/MarkovProcess/MixtureModel.hs
@@ -14,11 +14,10 @@
 module ELynx.Data.MarkovProcess.MixtureModel
   ( -- * Types
     Weight,
-    Component,
-    MixtureModel (name),
+    Component (weight, substModel),
+    MixtureModel (name, alphabet, components),
 
     -- * Getters
-    getAlphabet,
     getWeights,
     getSubstitutionModels,
 
@@ -30,26 +29,14 @@
     scale,
     normalize,
     appendNameComponents,
-
-    -- * Output
-    summarizeComponent,
-    summarize,
   )
 where
 
-import qualified Data.ByteString.Builder as L
-import qualified Data.ByteString.Lazy.Char8 as L
-import Data.List.NonEmpty hiding (zip)
+import qualified Data.List.NonEmpty as N
 import Data.Semigroup
-import ELynx.Data.Alphabet.Alphabet
+import ELynx.Data.Alphabet.Alphabet hiding (all)
 import qualified ELynx.Data.MarkovProcess.SubstitutionModel as S
-import ELynx.Tools
-import Prelude hiding
-  ( head,
-    length,
-    map,
-    zipWith,
-  )
+import Prelude
 
 -- | Mixture model component weight.
 type Weight = Double
@@ -66,43 +53,40 @@
   { -- | Name
     name :: S.Name,
     alphabet :: Alphabet,
-    components :: NonEmpty Component
+    components :: N.NonEmpty Component
   }
   deriving (Show, Read)
 
--- | Get alphabet used with mixture model. Throws error if components use
--- different 'Alphabet's.
-getAlphabet :: MixtureModel -> Alphabet
-getAlphabet = alphabet
-
 -- | Get weights.
-getWeights :: MixtureModel -> NonEmpty Weight
-getWeights = map weight . components
+getWeights :: MixtureModel -> N.NonEmpty Weight
+getWeights = N.map weight . components
 
 -- | Get substitution models.
-getSubstitutionModels :: MixtureModel -> NonEmpty S.SubstitutionModel
-getSubstitutionModels = map substModel . components
+getSubstitutionModels :: MixtureModel -> N.NonEmpty S.SubstitutionModel
+getSubstitutionModels = N.map substModel . components
 
 -- | Create a mixture model from a list of substitution models.
 fromSubstitutionModels ::
-  S.Name -> NonEmpty Weight -> NonEmpty S.SubstitutionModel -> MixtureModel
+  S.Name -> N.NonEmpty Weight -> N.NonEmpty S.SubstitutionModel -> MixtureModel
 fromSubstitutionModels n ws sms =
-  if allEqual $ toList alphs
-    then MixtureModel n (head alphs) comps
+  if allEqual $ N.toList alphs
+    then MixtureModel n (N.head alphs) comps
     else
       error
         "fromSubstitutionModels: alphabets of substitution models are not equal."
   where
-    comps = zipWith Component ws sms
-    alphs = map S.alphabet sms
+    comps = N.zipWith Component ws sms
+    alphs = N.map S.alphabet sms
+    allEqual [] = True
+    allEqual xs = all (== head xs) $ tail xs
 
 -- | Concatenate mixture models.
-concatenate :: S.Name -> NonEmpty MixtureModel -> MixtureModel
+concatenate :: S.Name -> N.NonEmpty MixtureModel -> MixtureModel
 concatenate n mms = fromSubstitutionModels n ws sms
   where
-    comps = sconcat $ map components mms
-    ws = map weight comps
-    sms = map substModel comps
+    comps = sconcat $ N.map components mms
+    ws = N.map weight comps
+    sms = N.map substModel comps
 
 scaleComponent :: Double -> Component -> Component
 scaleComponent s c = c {substModel = s'} where s' = S.scale s $ substModel c
@@ -112,16 +96,16 @@
 scale s m = m {components = cs'}
   where
     cs = components m
-    cs' = map (scaleComponent s) cs
+    cs' = N.map (scaleComponent s) cs
 
 -- | Globally normalize a mixture model so that on average one event happens per
 -- unit time.
 normalize :: MixtureModel -> MixtureModel
 normalize mm = scale (1 / c) mm
   where
-    c = sum $ zipWith (*) weights scales
+    c = sum $ N.zipWith (*) weights scales
     weights = getWeights mm
-    scales = map S.totalRate $ getSubstitutionModels mm
+    scales = N.map S.totalRate $ getSubstitutionModels mm
 
 appendNameComponent :: S.Name -> Component -> Component
 appendNameComponent n c = c {substModel = s'}
@@ -133,29 +117,4 @@
 appendNameComponents n m = m {components = cs'}
   where
     cs = components m
-    cs' = map (appendNameComponent n) cs
-
--- | Summarize a mixture model component; lines to be printed to screen or log.
-summarizeComponent :: Component -> [L.ByteString]
-summarizeComponent c =
-  L.pack "Weight: "
-    <> (L.toLazyByteString . L.doubleDec $ weight c) :
-  S.summarize (substModel c)
-
--- | Summarize a mixture model; lines to be printed to screen or log.
-summarize :: MixtureModel -> [L.ByteString]
-summarize m =
-  [ L.pack $ "Mixture model: " ++ name m ++ ".",
-    L.pack $ "Number of components: " ++ show n ++ "."
-  ]
-    ++ detail
-  where
-    n = length $ components m
-    detail =
-      if n <= 100
-        then
-          concat
-            [ L.pack ("Component " ++ show i ++ ":") : summarizeComponent c
-              | (i, c) <- zip [1 :: Int ..] (toList $ components m)
-            ]
-        else []
+    cs' = N.map (appendNameComponent n) cs
diff --git a/src/ELynx/Data/MarkovProcess/Nucleotide.hs b/src/ELynx/Data/MarkovProcess/Nucleotide.hs
--- a/src/ELynx/Data/MarkovProcess/Nucleotide.hs
+++ b/src/ELynx/Data/MarkovProcess/Nucleotide.hs
@@ -25,10 +25,10 @@
   )
 where
 
+import qualified Data.Vector.Storable as V
 import ELynx.Data.Alphabet.Alphabet
 import ELynx.Data.MarkovProcess.RateMatrix
 import ELynx.Data.MarkovProcess.SubstitutionModel
-import ELynx.Tools
 import Numeric.LinearAlgebra hiding (normalize)
 
 -- XXX: Another idea of structuring the code. This would probably be cleaner in
@@ -77,9 +77,12 @@
       0.0
     ]
 
+uniformVec :: Vector Double
+uniformVec = V.replicate n (1 / fromIntegral n)
+
 -- | JC substitution model.
 jc :: SubstitutionModel
-jc = substitutionModel DNA "JC" [] d jcExch where d = uniformVec n
+jc = substitutionModel DNA "JC" [] uniformVec jcExch
 
 -- | F81 substitution model.
 f81 :: StationaryDistribution -> SubstitutionModel
diff --git a/src/ELynx/Data/MarkovProcess/PhyloModel.hs b/src/ELynx/Data/MarkovProcess/PhyloModel.hs
--- a/src/ELynx/Data/MarkovProcess/PhyloModel.hs
+++ b/src/ELynx/Data/MarkovProcess/PhyloModel.hs
@@ -18,11 +18,9 @@
 module ELynx.Data.MarkovProcess.PhyloModel
   ( PhyloModel (..),
     getAlphabet,
-    summarize,
   )
 where
 
-import qualified Data.ByteString.Lazy.Char8 as L
 import ELynx.Data.Alphabet.Alphabet
 import qualified ELynx.Data.MarkovProcess.MixtureModel as M
 import qualified ELynx.Data.MarkovProcess.SubstitutionModel as S
@@ -34,10 +32,5 @@
 
 -- | Extract code from phylogenetic model.
 getAlphabet :: PhyloModel -> Alphabet
-getAlphabet (MixtureModel mm) = M.getAlphabet mm
+getAlphabet (MixtureModel mm) = M.alphabet mm
 getAlphabet (SubstitutionModel sm) = S.alphabet sm
-
--- | Summarize a phylogenetic model; lines to be printed to screen or log.
-summarize :: PhyloModel -> [L.ByteString]
-summarize (MixtureModel mm) = M.summarize mm
-summarize (SubstitutionModel sm) = S.summarize sm
diff --git a/src/ELynx/Data/MarkovProcess/RateMatrix.hs b/src/ELynx/Data/MarkovProcess/RateMatrix.hs
--- a/src/ELynx/Data/MarkovProcess/RateMatrix.hs
+++ b/src/ELynx/Data/MarkovProcess/RateMatrix.hs
@@ -34,7 +34,7 @@
   )
 where
 
-import ELynx.Tools
+import qualified Data.Vector.Storable as V
 import Numeric.LinearAlgebra hiding (normalize)
 import Numeric.SpecFunctions
 import Prelude hiding ((<>))
@@ -50,14 +50,21 @@
 -- | Stationary distribution of a rate matrix.
 type StationaryDistribution = Vector R
 
+epsRelaxed :: Double
+epsRelaxed = 1e-5
+
 -- | True if distribution sums to 1.0.
 isValid :: StationaryDistribution -> Bool
-isValid d = nearlyEqWith eps' (norm_1 d) 1.0
+isValid d = epsRelaxed > abs (norm_1 d - 1.0)
 
 -- | Normalize a stationary distribution so that the elements sum to 1.0.
 normalizeSD :: StationaryDistribution -> StationaryDistribution
 normalizeSD d = d / scalar (norm_1 d)
 
+matrixSetDiagToZero :: Matrix R -> Matrix R
+matrixSetDiagToZero m = m - diag (takeDiag m)
+{-# INLINE matrixSetDiagToZero #-}
+
 -- | Get average number of substitutions per unit time.
 totalRateWith :: StationaryDistribution -> RateMatrix -> Double
 totalRateWith d m = norm_1 $ d <# matrixSetDiagToZero m
@@ -95,6 +102,15 @@
   ExchangeabilityMatrix -> StationaryDistribution -> RateMatrix
 fromExchangeabilityMatrix em d = setDiagonal $ em <> diag d
 
+eps :: Double
+eps = 1e-12
+
+normalizeSumVec :: V.Vector Double -> V.Vector Double
+normalizeSumVec v = V.map (/ s) v
+  where
+    s = V.sum v
+{-# INLINE normalizeSumVec #-}
+
 -- | Get stationary distribution from 'RateMatrix'. Involves eigendecomposition.
 -- If the given matrix does not satisfy the required properties of transition
 -- rate matrices and no eigenvector with an eigenvalue nearly equal to 0 is
@@ -104,8 +120,8 @@
 -- function)?
 getStationaryDistribution :: RateMatrix -> StationaryDistribution
 getStationaryDistribution m =
-  if magnitude (eVals ! i) `nearlyEq` 0
-    then normalizeSumVec 1.0 distReal
+  if eps > abs (magnitude (eVals ! i))
+    then normalizeSumVec distReal
     else
       error
         "getStationaryDistribution: Could not retrieve stationary distribution."
diff --git a/src/ELynx/Data/MarkovProcess/SubstitutionModel.hs b/src/ELynx/Data/MarkovProcess/SubstitutionModel.hs
--- a/src/ELynx/Data/MarkovProcess/SubstitutionModel.hs
+++ b/src/ELynx/Data/MarkovProcess/SubstitutionModel.hs
@@ -17,9 +17,10 @@
     Params,
     SubstitutionModel,
 
-    -- * Lenses and other accessors
+    -- * Accessors
     alphabet,
     name,
+    params,
     stationaryDistribution,
     exchangeabilityMatrix,
     rateMatrix,
@@ -32,16 +33,12 @@
     scale,
     normalize,
     appendName,
-
-    -- * Output
-    summarize,
   )
 where
 
-import qualified Data.ByteString.Lazy.Char8 as L
+import qualified Data.Vector.Storable as V
 import ELynx.Data.Alphabet.Alphabet
 import qualified ELynx.Data.MarkovProcess.RateMatrix as R
-import ELynx.Tools
 import qualified Numeric.LinearAlgebra as LinAlg
 
 -- | Name of substitution model; abstracted and subject to change.
@@ -89,6 +86,12 @@
 totalRate :: SubstitutionModel -> Double
 totalRate sm = R.totalRate (rateMatrix sm)
 
+normalizeSumVec :: V.Vector Double -> V.Vector Double
+normalizeSumVec v = V.map (/ s) v
+  where
+    s = V.sum v
+{-# INLINE normalizeSumVec #-}
+
 -- | Create normalized 'SubstitutionModel'. See 'normalize'.
 substitutionModel ::
   Alphabet ->
@@ -99,13 +102,11 @@
   SubstitutionModel
 substitutionModel c n ps d e =
   if R.isValid d
-    then normalize $ SubstitutionModel c n ps d' e
+    then normalize $ SubstitutionModel c n ps (normalizeSumVec d) e
     else
       error $
         "substitionModel: Stationary distribution does not sum to 1.0: "
           ++ show d
-  where
-    d' = normalizeSumVec 1.0 d
 
 -- | Scale the rate of a substitution model by given factor.
 scale :: Double -> SubstitutionModel -> SubstitutionModel
@@ -121,29 +122,3 @@
 -- | Abbend to name.
 appendName :: Name -> SubstitutionModel -> SubstitutionModel
 appendName n sm = sm {name = n'} where n' = name sm <> n
-
--- | Summarize a substitution model; lines to be printed to screen or log.
-summarize :: SubstitutionModel -> [L.ByteString]
-summarize sm =
-  map L.pack $
-    (show (alphabet sm) ++ " substitution model: " ++ name sm ++ ".") :
-    ["Parameters: " ++ show (params sm) ++ "." | not (null (params sm))]
-      ++ case alphabet sm of
-        DNA ->
-          [ "Stationary distribution: "
-              ++ dispv precision (stationaryDistribution sm)
-              ++ ".",
-            "Exchangeability matrix:\n"
-              ++ dispmi 2 precision (exchangeabilityMatrix sm)
-              ++ ".",
-            "Scale: " ++ show (roundN precision $ totalRate sm) ++ "."
-          ]
-        Protein ->
-          [ "Stationary distribution: "
-              ++ dispv precision (stationaryDistribution sm)
-              ++ ".",
-            "Scale: " ++ show (roundN precision $ totalRate sm) ++ "."
-          ]
-        _ ->
-          error
-            "Extended character sets are not supported with substitution models."
diff --git a/src/ELynx/Import/MarkovProcess/EDMModelPhylobayes.hs b/src/ELynx/Import/MarkovProcess/EDMModelPhylobayes.hs
--- a/src/ELynx/Import/MarkovProcess/EDMModelPhylobayes.hs
+++ b/src/ELynx/Import/MarkovProcess/EDMModelPhylobayes.hs
@@ -1,3 +1,5 @@
+{-# LANGUAGE OverloadedStrings #-}
+
 -- |
 -- Module      :  ELynx.Import.MarkovProcess.EDMModelPhylobayes
 -- Description :  Import stationary distributions from Phylobayes format
@@ -10,63 +12,53 @@
 --
 -- Creation date: Tue Jan 29 12:12:55 2019.
 module ELynx.Import.MarkovProcess.EDMModelPhylobayes
-  ( Parser,
-    EDMComponent,
+  ( EDMComponent,
     phylobayes,
   )
 where
 
+import Control.Applicative
 import Control.Monad
-import qualified Data.ByteString.Lazy.Char8 as L
+import qualified Data.Attoparsec.ByteString as AS
+import qualified Data.Attoparsec.ByteString.Char8 as AC
+import qualified Data.ByteString.Lazy.Char8 as BL
 import qualified Data.Vector.Storable as V
-import Data.Void
 import ELynx.Data.MarkovProcess.MixtureModel
-  ( Weight,
-  )
-import ELynx.Tools
-import Text.Megaparsec
-import Text.Megaparsec.Byte
-import Text.Megaparsec.Byte.Lexer hiding (space)
 
--- | Shortcut.
-type Parser = Parsec Void L.ByteString
-
 -- | An empirical mixture model component has a weight and a stationary
 -- distribution.
 type EDMComponent = (Weight, V.Vector Double)
 
 -- | Parse stationary distributions from Phylobayes format.
-phylobayes :: Parser [EDMComponent]
-phylobayes = do
+phylobayes :: AS.Parser [EDMComponent]
+phylobayes = (AS.<?> "phylobayes") $ do
   n <- headerLine
   k <- kComponentsLine
-  cs <- count k $ dataLine n
-  _ <- many newline *> eof <?> "phylobayes"
+  cs <- AS.count k $ dataLine n
+  _ <- AC.skipWhile AC.isSpace
+  _ <- AS.endOfInput
   return cs
 
-horizontalSpace :: Parser ()
-horizontalSpace = skipMany $ char (c2w ' ') <|> tab
-
-headerLine :: Parser Int
+headerLine :: AS.Parser Int
 headerLine = do
-  n <- decimal
-  _ <- horizontalSpace
+  n <- AC.decimal
+  _ <- AS.skipWhile AC.isHorizontalSpace
   -- XXX: This should be more general, but then we also want to ensure that the
   -- order of states is correct.
   _ <-
-    chunk (L.pack "A C D E F G H I K L M N P Q R S T V W Y")
-      <|> chunk (L.pack "A C G T")
-  _ <- many newline <?> "headerLine"
+    AS.string (BL.toStrict "A C D E F G H I K L M N P Q R S T V W Y")
+      <|> AS.string (BL.toStrict "A C G T")
+  _ <- AC.skipWhile AC.isSpace AS.<?> "headerLine"
   return n
 
-kComponentsLine :: Parser Int
-kComponentsLine = decimal <* newline <?> "kComponentsLine"
+kComponentsLine :: AS.Parser Int
+kComponentsLine = AC.decimal <* AC.skipWhile AC.isSpace AS.<?> "kComponentsLine"
 
-dataLine :: Int -> Parser EDMComponent
-dataLine n = do
-  weight <- float
-  _ <- horizontalSpace
-  vals <- float `sepEndBy1` horizontalSpace
+dataLine :: Int -> AS.Parser EDMComponent
+dataLine n = (AS.<?> "dataLine") $ do
+  wght <- AC.double
+  _ <- AS.skipWhile AC.isHorizontalSpace
+  vals <- AC.double `AC.sepBy1` AS.skipWhile AC.isHorizontalSpace
+  _ <- AC.skipWhile AC.isSpace
   when (length vals /= n) (error "Did not find correct number of entries.")
-  _ <- space <?> "dataLine"
-  return (weight, V.fromList vals)
+  return (wght, V.fromList vals)
diff --git a/src/ELynx/Import/MarkovProcess/SiteprofilesPhylobayes.hs b/src/ELynx/Import/MarkovProcess/SiteprofilesPhylobayes.hs
--- a/src/ELynx/Import/MarkovProcess/SiteprofilesPhylobayes.hs
+++ b/src/ELynx/Import/MarkovProcess/SiteprofilesPhylobayes.hs
@@ -21,24 +21,21 @@
   )
 where
 
+import Control.Applicative
 import Control.Monad
+import qualified Data.Attoparsec.ByteString as AS
+import qualified Data.Attoparsec.ByteString.Char8 as AC
 import Data.Containers.ListUtils (nubInt)
 import qualified Data.Vector.Storable as V
 import ELynx.Import.MarkovProcess.EDMModelPhylobayes
-  ( EDMComponent,
-    Parser,
-  )
-import ELynx.Tools
-import Text.Megaparsec
-import Text.Megaparsec.Byte
-import Text.Megaparsec.Byte.Lexer
 
 -- | Parse stationary distributions from Phylobayes format.
-siteprofiles :: Parser [EDMComponent]
-siteprofiles = do
+siteprofiles :: AS.Parser [EDMComponent]
+siteprofiles = (AS.<?> "phylobayes siteprofiles") $ do
   _ <- headerLines
   cs <- many dataLine
-  _ <- many newline *> eof <?> "phylobayes siteprofiles"
+  _ <- AC.skipWhile AC.isSpace
+  _ <- AS.endOfInput
   let ls = map length cs
       nLs = length $ nubInt ls
   when
@@ -46,28 +43,20 @@
     (error "The site profiles have a different number of entries.")
   return cs
 
-horizontalSpace :: Parser ()
-horizontalSpace = skipMany $ char (c2w ' ') <|> tab
-
-line :: Parser ()
-line = do
-  _ <- many $ noneOf [c2w '\n']
-  pure ()
+line :: AS.Parser ()
+line = AS.skipWhile (not . AC.isEndOfLine)
 
 -- For now, just ignore the header.
-headerLines :: Parser ()
-headerLines = do
-  _ <- line
-  _ <- many newline <?> "headerLine"
-  pure ()
+headerLines :: AS.Parser ()
+headerLines = line *> AC.skipWhile AC.isSpace AS.<?> "headerLine"
 
-dataLine :: Parser EDMComponent
-dataLine = do
+dataLine :: AS.Parser EDMComponent
+dataLine = (AS.<?> "dataLine") $ do
   -- Ignore site number.
-  _ <- decimal :: Parser Integer
-  _ <- horizontalSpace
+  _ <- AC.decimal :: AS.Parser Int
+  _ <- AS.skipWhile AC.isHorizontalSpace
   -- Also ignore additional white space on line.
-  vals <- float `sepEndBy1` horizontalSpace
-  _ <- many newline <?> "dataLine"
+  vals <- AC.double `AS.sepBy1` AS.skipWhile AC.isHorizontalSpace
+  _ <- AC.skipWhile AC.isSpace
   -- Set the weight to 1.0 for all sites.
   return (1.0, V.fromList vals)
diff --git a/src/ELynx/Simulate/MarkovProcessAlongTree.hs b/src/ELynx/Simulate/MarkovProcessAlongTree.hs
--- a/src/ELynx/Simulate/MarkovProcessAlongTree.hs
+++ b/src/ELynx/Simulate/MarkovProcessAlongTree.hs
@@ -1,3 +1,5 @@
+{-# LANGUAGE ScopedTypeVariables #-}
+
 -- |
 --   Description :  Work with transition probability matrices on rooted trees
 --   Copyright   :  (c) Dominik Schrempf 2017
@@ -24,15 +26,18 @@
   )
 where
 
+import Control.Concurrent
+import Control.Concurrent.Async
 import Control.Monad
 import Control.Monad.Primitive
 import Control.Parallel.Strategies
 import Data.Tree
+import qualified Data.Vector.Storable as V
+import Data.Word (Word32)
 import ELynx.Data.MarkovProcess.RateMatrix
 import ELynx.Simulate.MarkovProcess
-import ELynx.Tools
 import Numeric.LinearAlgebra
-import System.Random.MWC (Gen, GenIO)
+import System.Random.MWC
 import System.Random.MWC.Distributions (categorical)
 
 toProbTree :: RateMatrix -> Tree Double -> Tree ProbMatrix
@@ -159,8 +164,28 @@
       concat
         <$> sequence [simulateAndFlattenMixtureModel' is' cs t g | t <- f]
 
--- | See 'simulateAndFlattenMixtureModel', parallel version;
--- needs to be run in IO monad.
+getChunks :: Int -> Int -> [Int]
+getChunks c n = ns
+  where
+    n' = n `div` c
+    r = n `mod` c
+    ns = replicate r (n' + 1) ++ replicate (c - r) n'
+
+splitGen :: PrimMonad m => Int -> Gen (PrimState m) -> m [Gen (PrimState m)]
+splitGen n gen
+  | n <= 0 = return []
+  | otherwise = do
+    seeds :: [V.Vector Word32] <- replicateM (n -1) $ uniformVector gen 256
+    fmap (gen :) (mapM initialize seeds)
+
+parComp :: Int -> (Int -> GenIO -> IO b) -> GenIO -> IO [b]
+parComp num fun gen = do
+  ncap <- getNumCapabilities
+  let chunks = getChunks ncap num
+  gs <- splitGen ncap gen
+  mapConcurrently (uncurry fun) (zip chunks gs)
+
+-- | See 'simulateAndFlattenMixtureModel', parallel version.
 simulateAndFlattenMixtureModelPar ::
   Int ->
   Vector R ->
diff --git a/test/ELynx/Data/MarkovProcess/AminoAcidSpec.hs b/test/ELynx/Data/MarkovProcess/AminoAcidSpec.hs
--- a/test/ELynx/Data/MarkovProcess/AminoAcidSpec.hs
+++ b/test/ELynx/Data/MarkovProcess/AminoAcidSpec.hs
@@ -13,6 +13,7 @@
   )
 where
 
+import qualified Data.Vector.Storable as V
 import ELynx.Data.MarkovProcess.AminoAcid
 import qualified ELynx.Data.MarkovProcess.RateMatrix as R
 import qualified ELynx.Data.MarkovProcess.SubstitutionModel as S
@@ -20,9 +21,15 @@
 import Numeric.LinearAlgebra
 import Test.Hspec
 
+normalizeSumVec :: Vector Double -> Vector Double
+normalizeSumVec v = V.map (/ s) v
+  where
+    s = V.sum v
+{-# INLINE normalizeSumVec #-}
+
 statDistLGPython :: R.StationaryDistribution
 statDistLGPython =
-  normalizeSumVec 1.0 $
+  normalizeSumVec $
     fromList
       [ 0.079066,
         0.012937,