diff --git a/ChangeLog.md b/ChangeLog.md
--- a/ChangeLog.md
+++ b/ChangeLog.md
@@ -5,6 +5,20 @@
 ## Unreleased changes
 
 
+## Version 0.4.1
+
+-   Improve `TimeSpec` (Point process).
+-   Parallel evaluation strategies.
+-   Change names of some functions involving partitions. For example, `mp` was
+    renamed to `pt`.
+-   Improve documentation for (bi)partitions.
+-   Bugfix `tlynx compare`; do not throw error when branch support values are not
+    set.
+-   Add `no-elynx-file` option.
+-   Also parse Nexus files with `tlynx` commands.
+-   Bugfix `subSample`; the sub sample was reversed.
+
+
 ## Version 0.4.0
 
 -   Major refactor of `elynx-tree`. All required function can now conveniently
diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -2,7 +2,7 @@
 
 # The ELynx Suite
 
-Version: 0.4.0.
+Version: 0.4.1.
 Reproducible evolution made easy.
 
 <p align="center"><img src="https://travis-ci.org/dschrempf/elynx.svg?branch=master"/></p>
@@ -90,12 +90,12 @@
 
     slynx --help | head -n -16
 
-    ELynx Suite version 0.4.0.
+    ELynx Suite version 0.4.1.
     Developed by Dominik Schrempf.
-    Compiled on September 4, 2020, at 13:37 pm, UTC.
+    Compiled on October 16, 2020, at 12:46 pm, UTC.
     
     Usage: slynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME] 
-                 [-f|--force] COMMAND
+                 [-f|--force] [--no-elynx-file] COMMAND
       Analyze, and simulate multi sequence alignments.
     
     Available options:
@@ -107,6 +107,7 @@
                                Specify base name of output file
       -f,--force               Ignore previous analysis and overwrite existing
                                output files.
+      --no-elynx-file          Do not write files for needed for reproducibility.
     
     Available commands:
       concatenate              Concatenate sequences found in input files.
@@ -136,9 +137,9 @@
 
     slynx concatenate --help
 
-    ELynx Suite version 0.4.0.
+    ELynx Suite version 0.4.1.
     Developed by Dominik Schrempf.
-    Compiled on September 4, 2020, at 13:37 pm, UTC.
+    Compiled on October 16, 2020, at 12:46 pm, UTC.
     
     Usage: slynx concatenate (-a|--alphabet NAME) INPUT-FILE
       Concatenate sequences found in input files.
@@ -157,9 +158,9 @@
 
     slynx examine --help
 
-    ELynx Suite version 0.4.0.
+    ELynx Suite version 0.4.1.
     Developed by Dominik Schrempf.
-    Compiled on September 4, 2020, at 13:37 pm, UTC.
+    Compiled on October 16, 2020, at 12:46 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 +180,9 @@
 
     slynx filter-rows --help
 
-    ELynx Suite version 0.4.0.
+    ELynx Suite version 0.4.1.
     Developed by Dominik Schrempf.
-    Compiled on September 4, 2020, at 13:37 pm, UTC.
+    Compiled on October 16, 2020, at 12:46 pm, UTC.
     
     Usage: slynx filter-rows (-a|--alphabet NAME) INPUT-FILE [--longer-than LENGTH] 
                              [--shorter-than LENGTH] [--standard-characters]
@@ -202,9 +203,9 @@
 
     slynx filter-columns --help
 
-    ELynx Suite version 0.4.0.
+    ELynx Suite version 0.4.1.
     Developed by Dominik Schrempf.
-    Compiled on September 4, 2020, at 13:37 pm, UTC.
+    Compiled on October 16, 2020, at 12:46 pm, UTC.
     
     Usage: slynx filter-columns (-a|--alphabet NAME) INPUT-FILE 
                                 [--standard-chars DOUBLE]
@@ -226,9 +227,9 @@
 
     slynx simulate --help
 
-    ELynx Suite version 0.4.0.
+    ELynx Suite version 0.4.1.
     Developed by Dominik Schrempf.
-    Compiled on September 4, 2020, at 13:37 pm, UTC.
+    Compiled on October 16, 2020, at 12:46 pm, UTC.
     
     Usage: slynx simulate (-t|--tree-file Name) [-s|--substitution-model MODEL] 
                           [-m|--mixture-model MODEL] [-e|--edm-file NAME] 
@@ -305,9 +306,9 @@
 
     slynx sub-sample --help
 
-    ELynx Suite version 0.4.0.
+    ELynx Suite version 0.4.1.
     Developed by Dominik Schrempf.
-    Compiled on September 4, 2020, at 13:37 pm, UTC.
+    Compiled on October 16, 2020, at 12:46 pm, UTC.
     
     Usage: slynx sub-sample (-a|--alphabet NAME) INPUT-FILE
                             (-n|--number-of-sites INT)
@@ -335,9 +336,9 @@
 
     slynx translate --help
 
-    ELynx Suite version 0.4.0.
+    ELynx Suite version 0.4.1.
     Developed by Dominik Schrempf.
-    Compiled on September 4, 2020, at 13:37 pm, UTC.
+    Compiled on October 16, 2020, at 12:46 pm, UTC.
     
     Usage: slynx translate (-a|--alphabet NAME) INPUT-FILE (-r|--reading-frame INT)
                            (-u|--universal-code CODE)
@@ -360,12 +361,12 @@
 
     tlynx --help | head -n -16
 
-    ELynx Suite version 0.4.0.
+    ELynx Suite version 0.4.1.
     Developed by Dominik Schrempf.
-    Compiled on September 4, 2020, at 13:37 pm, UTC.
+    Compiled on October 16, 2020, at 12:46 pm, UTC.
     
     Usage: tlynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME] 
-                 [-f|--force] COMMAND
+                 [-f|--force] [--no-elynx-file] COMMAND
       Compare, examine, and simulate phylogenetic trees.
     
     Available options:
@@ -377,6 +378,7 @@
                                Specify base name of output file
       -f,--force               Ignore previous analysis and overwrite existing
                                output files.
+      --no-elynx-file          Do not write files for needed for reproducibility.
     
     Available commands:
       compare                  Compare two phylogenetic trees (compute distances and branch-wise differences).
@@ -390,6 +392,7 @@
     Available tree file formats:
       - Newick Standard: Branch support values are stored in square brackets after branch lengths.
       - Newick IqTree:   Branch support values are stored as node names after the closing bracket of forests.
+      - Newick RevBayes: Key-value pairs is provided in square brackets after node names as well as branch lengths. XXX: Key value pairs are ignored at the moment.
 
 
 ## Compare
@@ -398,9 +401,9 @@
 
     tlynx compare --help
 
-    ELynx Suite version 0.4.0.
+    ELynx Suite version 0.4.1.
     Developed by Dominik Schrempf.
-    Compiled on September 4, 2020, at 13:37 pm, UTC.
+    Compiled on October 16, 2020, at 12:46 pm, UTC.
     
     Usage: tlynx compare [-n|--normalize] [-b|--bipartitions] [-t|--intersect] 
                          [-f|--newick-format FORMAT] NAMES
@@ -427,9 +430,9 @@
 
     tlynx examine --help
 
-    ELynx Suite version 0.4.0.
+    ELynx Suite version 0.4.1.
     Developed by Dominik Schrempf.
-    Compiled on September 4, 2020, at 13:37 pm, UTC.
+    Compiled on October 16, 2020, at 12:46 pm, UTC.
     
     Usage: tlynx examine INPUT-FILE [-f|--newick-format FORMAT]
       Compute summary statistics of phylogenetic trees.
@@ -451,9 +454,9 @@
 
     tlynx simulate --help
 
-    ELynx Suite version 0.4.0.
+    ELynx Suite version 0.4.1.
     Developed by Dominik Schrempf.
-    Compiled on September 4, 2020, at 13:37 pm, UTC.
+    Compiled on October 16, 2020, at 12:46 pm, UTC.
     
     Usage: tlynx simulate (-t|--nTrees INT) (-n|--nLeaves INT) PROCESS 
                           [-u|--sub-sample DOUBLE] [-s|--summary-statistics] 
@@ -488,9 +491,9 @@
 
     elynx --help | head -n -16
 
-    ELynx Suite version 0.4.0.
+    ELynx Suite version 0.4.1.
     Developed by Dominik Schrempf.
-    Compiled on September 4, 2020, at 13:37 pm, UTC.
+    Compiled on October 16, 2020, at 12:46 pm, UTC.
     
     Usage: elynx COMMAND
       Validate and redo past ELynx analyses
diff --git a/bench/Bench.hs b/bench/Bench.hs
--- a/bench/Bench.hs
+++ b/bench/Bench.hs
@@ -1,3 +1,5 @@
+{-# LANGUAGE BangPatterns #-}
+
 -- |
 -- Module      :  Bench
 -- Description :  Various benchmarks
@@ -11,11 +13,15 @@
 -- Creation date: Mon Dec 16 13:33:27 2019.
 module Main where
 
+import Control.Parallel.Strategies
 import Criterion.Main
+import Data.Bifunctor
 import qualified Data.ByteString.Char8 as BS
-import ELynx.Tools
+import Data.Foldable
+import ELynx.Tools hiding (Random)
 import ELynx.Tree
-import ELynx.Tree.Import.Newick
+import ELynx.Tree.Simulate.PointProcess
+import System.Random.MWC
 
 treeFileMany :: FilePath
 treeFileMany = "data/Many.trees"
@@ -23,8 +29,48 @@
 getManyTrees :: IO (Forest Phylo BS.ByteString)
 getManyTrees = parseFileWith (someNewick Standard) treeFileMany
 
+hugeTree :: IO (Tree Length Int)
+hugeTree = create >>= simulateReconstructedTree 50000 Random 1.0 0.9
+
+sinN :: Int -> Double -> Double
+sinN n x = iterate sin x !! n
+
+func :: Double -> Double
+func = sinN 200
+
+hugeTreeCalcPar :: Int -> Tree Double Int -> Tree Double Int
+hugeTreeCalcPar n t = first func t `using` parTree n
+
 main :: IO ()
 main = do
-  ts <- getManyTrees
+  !ts <- getManyTrees
+  !ht <- first getLen <$> hugeTree
+  let mr1 = hugeTreeCalcPar 0 ht
+      mr2 = hugeTreeCalcPar 1 ht
+  if mr1 == mr2
+    then putStrLn "Map OK."
+    else do
+      print mr1
+      print mr2
+      error "Map wrong."
+  let fr1 = (foldl' (+) 0 . branches) ht
+      fr2 = parBranchFoldMap 1 id (+) ht
+  if 1e-8 > abs (fr1 - fr2)
+    then putStrLn "Fold OK."
+    else do
+      print fr1
+      print fr2
+      error "Fold wrong."
   defaultMain
-    [bgroup "bipartition" [bench "manyTrees" $ nf (map bipartitions) ts]]
+    [ bgroup "bipartition" [bench "manyTrees" $ nf (map bipartitions) ts],
+      bgroup
+        "map strategies"
+        [ bench "sequential" $ nf (hugeTreeCalcPar 0) ht,
+          bench "parallel 3" $ nf (hugeTreeCalcPar 3) ht
+        ],
+      bgroup
+        "fold strategies"
+        [ bench "sequential" $ nf (parBranchFoldMap 0 func (+)) ht,
+          bench "parallel 3" $ nf (parBranchFoldMap 3 func (+)) ht
+        ]
+    ]
diff --git a/data/UltraMetric.log b/data/UltraMetric.log
deleted file mode 100644
--- a/data/UltraMetric.log
+++ /dev/null
@@ -1,22 +0,0 @@
----------------------------------------------------------------
--- Simulate phylogenetic trees using birth and death processes.
-ELynx Suite version 0.0.1.
-Developed by Dominik Schrempf.
-Compiled on January 13, 2020, at 16:42 pm, UTC.
-Start time: January 14, 2020, at 09:35 am, UTC.
-Command line: tlynx -o UltraMetric simulate -t 1 -n 10 -H 1.0
--- Arguments
-Number of simulated trees: 1
-Number of leaves per tree: 10
-Height of trees: 1.0, conditioned on origin
-Birth rate: 1.0
-Death rate: 0.9
-Sampling probability: 1.0
-Perform sub-sampling: False
-Summary statistics only: False
-Seed: Random
--- Simulation
-Number of used cores: 4
-Write simulated trees to file 'UltraMetric.tree'.
--- End time: January 14, 2020, at 09:35 am, UTC.
-------------------------------------------------
diff --git a/elynx-tree.cabal b/elynx-tree.cabal
--- a/elynx-tree.cabal
+++ b/elynx-tree.cabal
@@ -1,6 +1,6 @@
 cabal-version:  2.2
 name:           elynx-tree
-version:        0.4.0
+version:        0.4.1
 synopsis:       Handle phylogenetic trees
 description:    Examine, compare, and simulate phylogenetic trees in a reproducible way. Please see the README on GitHub at <https://github.com/dschrempf/elynx>.
 category:       Bioinformatics
@@ -31,7 +31,6 @@
     data/RevBayes.trees
     data/SimpleTree.nex
     data/TreeDist.trees
-    data/UltraMetric.log
     data/UltraMetric.tree
 
 source-repository head
@@ -47,6 +46,7 @@
       ELynx.Tree.Distance
       ELynx.Tree.Measurable
       ELynx.Tree.Named
+      ELynx.Tree.Parallel
       ELynx.Tree.Partition
       ELynx.Tree.Phylogeny
       ELynx.Tree.Rooted
@@ -86,6 +86,7 @@
     , elynx-nexus
     , math-functions
     , mwc-random
+    , parallel
     , primitive
     , statistics
   default-language: Haskell2010
@@ -128,11 +129,15 @@
       Paths_elynx_tree
   hs-source-dirs:
       bench
-  ghc-options: -Wall
+  ghc-options: -Wall -threaded -rtsopts -with-rtsopts=-N
+  -- ghc-options: -Wall
   build-depends:
       base >=4.7 && <5
     , bytestring
     , criterion
+    , deepseq
     , elynx-tools
     , elynx-tree
+    , mwc-random
+    , parallel
   default-language: Haskell2010
diff --git a/src/ELynx/Tree.hs b/src/ELynx/Tree.hs
--- a/src/ELynx/Tree.hs
+++ b/src/ELynx/Tree.hs
@@ -15,8 +15,13 @@
 module ELynx.Tree
   ( -- * Rooted trees
     module ELynx.Tree.Rooted,
+
+    -- * Tree zipper
     module ELynx.Tree.Zipper,
 
+    -- * Parallel evaluation
+    module ELynx.Tree.Parallel,
+
     -- * Branch label classes
     module ELynx.Tree.Measurable,
     module ELynx.Tree.Splittable,
@@ -49,6 +54,7 @@
 import ELynx.Tree.Import.Nexus
 import ELynx.Tree.Measurable
 import ELynx.Tree.Named
+import ELynx.Tree.Parallel
 import ELynx.Tree.Partition
 import ELynx.Tree.Phylogeny
 import ELynx.Tree.Rooted
diff --git a/src/ELynx/Tree/Bipartition.hs b/src/ELynx/Tree/Bipartition.hs
--- a/src/ELynx/Tree/Bipartition.hs
+++ b/src/ELynx/Tree/Bipartition.hs
@@ -56,9 +56,13 @@
 -- I am proud of this awesome 'Comonad' usage here :).
 groups = extend leaves
 
--- | Each branch of a tree partitions the leaves of the tree into two subsets,
--- or a bipartition.
+-- | A bipartition of a tree is a grouping of the leaves of the tree into two
+-- non-overlapping, non-empty sub sets.
 --
+-- For example, each branch of a tree partitions the leaves of the tree into two
+-- subsets, or a bipartition. Also, a bifurcating root induces a bipartition;
+-- see 'bipartition'.
+--
 -- The order of the two subsets of a 'Bipartition' is meaningless. We ensure by
 -- construction that the smaller subset comes first, and hence, that equality
 -- checks are meaningful.
@@ -152,8 +156,6 @@
       [bipartitions' c s | (c, s) <- zip cs ts]
   where
     cs = getComplementaryLeaves p t
-
--- TODO: Unrooted? See module comment of Distance.hs.
 
 -- | Convert a tree into a 'Map' from each 'Bipartition' to the branch inducing
 -- the respective 'Bipartition'.
diff --git a/src/ELynx/Tree/Distance.hs b/src/ELynx/Tree/Distance.hs
--- a/src/ELynx/Tree/Distance.hs
+++ b/src/ELynx/Tree/Distance.hs
@@ -10,20 +10,18 @@
 --
 -- Creation date: Thu Jun 13 17:15:54 2019.
 --
--- Various distance functions for phylogenetic trees (and trees with branch
--- lengths in general).
+-- Various distance functions for trees.
 --
 -- The functions provided in this module return distances for __unrooted__
 -- trees. See comments of 'symmetric', 'branchScore', and 'bipartitionToBranch',
 -- as well as the documentation of
 -- [treedist](http://evolution.genetics.washington.edu/phylip/doc/treedist.html).
 --
--- It is a little unfortunate that 'Tree' data type, which represents rooted
--- trees, is also used in this module. However, rooted trees are much easier to
--- handle. In the future, a separate data type for unrooted trees may be
--- introduced. In theory, this is quite straight forward, for example, using
--- algebraic graphs. Difficulties may arise because the branches of an unrooted
--- tree are undirected.
+-- It is a little unfortunate that the 'Tree' data type represents rooted trees.
+-- However, rooted trees are much easier to handle computationally. In the
+-- future, a separate data type for unrooted trees may be introduced, for
+-- example, using algebraic graphs. Difficulties may arise because the branches
+-- of an unrooted tree are undirected.
 module ELynx.Tree.Distance
   ( symmetric,
     incompatibleSplits,
@@ -66,7 +64,7 @@
 
 countIncompatibilities :: (Show a, Ord a) => Set (Bipartition a) -> Set (Partition a) -> Int
 countIncompatibilities bs ms =
-  foldl' (\i b -> if any (compatible $ bpToMp b) ms then i else i + 1) 0 bs
+  foldl' (\i b -> if any (compatible $ bpToPt b) ms then i else i + 1) 0 bs
 
 -- | Number of incompatible splits.
 --
diff --git a/src/ELynx/Tree/Export/Newick.hs b/src/ELynx/Tree/Export/Newick.hs
--- a/src/ELynx/Tree/Export/Newick.hs
+++ b/src/ELynx/Tree/Export/Newick.hs
@@ -27,6 +27,12 @@
 import ELynx.Tree.Phylogeny
 import ELynx.Tree.Rooted
 
+buildBrSup :: Double -> BB.Builder
+buildBrSup bs = BB.char8 '[' <> BB.doubleDec bs <> BB.char8 ']'
+
+buildBrLen :: Double -> BB.Builder
+buildBrLen bl = BB.char8 ':' <> BB.doubleDec bl
+
 -- | See 'toNewick'.
 toNewickBuilder :: Named a => Tree Phylo a -> BB.Builder
 toNewickBuilder t = go t <> BB.char8 ';'
@@ -37,8 +43,8 @@
         <> mconcat (intersperse (BB.char8 ',') $ map go ts)
         <> BB.char8 ')'
         <> lbl b l
-    mBrSupBuilder x = maybe mempty (\bs -> BB.char8 '[' <> BB.doubleDec bs <> BB.char8 ']') (brSup x)
-    mBrLenBuilder x = maybe mempty (\bl -> BB.char8 ':' <> BB.doubleDec bl) (brLen x)
+    mBrSupBuilder x = maybe mempty buildBrSup (brSup x)
+    mBrLenBuilder x = maybe mempty buildBrLen (brLen x)
     lbl x y =
       BB.lazyByteString (getName y)
         <> mBrLenBuilder x
diff --git a/src/ELynx/Tree/Import/Newick.hs b/src/ELynx/Tree/Import/Newick.hs
--- a/src/ELynx/Tree/Import/Newick.hs
+++ b/src/ELynx/Tree/Import/Newick.hs
@@ -1,5 +1,6 @@
 {-# LANGUAGE DeriveGeneric #-}
 
+-- |
 -- Module      :  ELynx.Tree.Import.Newick
 -- Description :  Import Newick trees
 -- Copyright   :  (c) Dominik Schrempf 2020
@@ -19,15 +20,24 @@
 -- In particular, no conversion from _ to (space) is done right now.
 --
 -- For a description of rooted 'Tree's, please see the 'ELynx.Tree.Rooted'
-
--- |
--- module header.
+--
+-- Code snippet:
+--
+-- @
+-- import Data.Attoparsec.ByteString
+-- import ELynx.Tree
+--
+-- getOneNewick = either error id . parseOnly (oneNewick Standard)
+-- @
 module ELynx.Tree.Import.Newick
   ( NewickFormat (..),
     description,
     newick,
+    parseNewick,
     oneNewick,
+    parseOneNewick,
     someNewick,
+    parseSomeNewick,
   )
 where
 
@@ -64,24 +74,36 @@
 description RevBayes =
   "RevBayes: Key-value pairs is provided in square brackets after node names as well as branch lengths. XXX: Key value pairs are ignored at the moment."
 
--- | Parse a single Newick tree. Also succeeds when more trees follow.
+-- | Newick tree parser. Also succeeds when more trees follow.
 newick :: NewickFormat -> Parser (Tree Phylo BS.ByteString)
 newick Standard = newickStandard
 newick IqTree = newickIqTree
 newick RevBayes = newickRevBayes
 
--- | Parse a single Newick tree. Fails when end of file is not reached.
+-- | See 'newick'.
+parseNewick :: NewickFormat -> BS.ByteString -> Tree Phylo BS.ByteString
+parseNewick f = either error id . parseOnly (newick f)
+
+-- | One Newick tree parser. Fails when end of input is not reached.
 oneNewick :: NewickFormat -> Parser (Tree Phylo BS.ByteString)
 oneNewick Standard = oneNewickStandard
 oneNewick IqTree = oneNewickIqTree
 oneNewick RevBayes = oneNewickRevBayes
 
--- | Parse one or more Newick trees until end of file.
+-- | See 'oneNewick'.
+parseOneNewick :: NewickFormat -> BS.ByteString -> Tree Phylo BS.ByteString
+parseOneNewick f = either error id . parseOnly (oneNewick f)
+
+-- | One or more Newick trees parser.
 someNewick :: NewickFormat -> Parser (Forest Phylo BS.ByteString)
 someNewick Standard = someNewickStandard
 someNewick IqTree = someNewickIqTree
 someNewick RevBayes = someNewickRevBayes
 
+-- | See 'someNewick'.
+parseSomeNewick :: NewickFormat -> BS.ByteString -> [Tree Phylo BS.ByteString]
+parseSomeNewick f = either error id . parseOnly (someNewick f)
+
 -- Parse a single Newick tree. Also succeeds when more trees follow.
 newickStandard :: Parser (Tree Phylo BS.ByteString)
 newickStandard = skipWhile isSpace *> tree <* char ';' <* skipWhile isSpace <?> "newickStandard"
@@ -183,10 +205,14 @@
 -- about nodes and branches. Parse a single Newick tree. Also succeeds when more
 -- trees follow.
 --
--- XXX: Key value pairs are ignored at the moment.
+-- TODO: Key value pairs are ignored at the moment.
 newickRevBayes :: Parser (Tree Phylo BS.ByteString)
 newickRevBayes =
-  skipWhile isSpace *> optional brackets *> treeRevBayes <* char ';' <* skipWhile isSpace <?> "newickRevBayes"
+  skipWhile isSpace
+    *> optional brackets
+    *> treeRevBayes
+    <* char ';'
+    <* skipWhile isSpace <?> "newickRevBayes"
 
 -- See 'newickRevBayes'. Parse a single Newick tree. Fails when end of file is
 -- not reached.
diff --git a/src/ELynx/Tree/Import/Nexus.hs b/src/ELynx/Tree/Import/Nexus.hs
--- a/src/ELynx/Tree/Import/Nexus.hs
+++ b/src/ELynx/Tree/Import/Nexus.hs
@@ -27,7 +27,7 @@
 
 -- | Parse a Nexus files with a TREES block.
 nexusTrees :: NewickFormat -> Parser [(BS.ByteString, Tree Phylo BS.ByteString)]
-nexusTrees = nexus . trees
+nexusTrees = nexusBlock . trees
 
 trees :: NewickFormat -> Block [(BS.ByteString, Tree Phylo BS.ByteString)]
 trees f = Block "TREES" (some $ namedNewick f)
@@ -35,11 +35,11 @@
 namedNewick :: NewickFormat -> Parser (BS.ByteString, Tree Phylo BS.ByteString)
 namedNewick f = do
   _ <- skipWhile isSpace
-  _ <- string "TREE"
+  _ <- stringCI "TREE" <?> "namedNewickTreeStart"
   _ <- skipWhile isSpace
-  n <- takeWhile1 (\x -> isAlpha_ascii x || isDigit x)
+  n <- takeWhile1 (\x -> isAlpha_ascii x || isDigit x) <?> "namedNewickTreeName"
   _ <- skipWhile isSpace
-  _ <- char '='
+  _ <- char '=' <?> "namedNewickEqual"
   _ <- skipWhile isSpace
-  t <- newick f
+  t <- newick f <?> "namedNewickTree"
   return (n, t)
diff --git a/src/ELynx/Tree/Parallel.hs b/src/ELynx/Tree/Parallel.hs
new file mode 100644
--- /dev/null
+++ b/src/ELynx/Tree/Parallel.hs
@@ -0,0 +1,62 @@
+-- |
+-- Module      :  ELynx.Tree.Parallel
+-- Description :  Evaluation strategies for trees
+-- Copyright   :  (c) Dominik Schrempf, 2020
+-- License     :  GPL-3.0-or-later
+--
+-- Maintainer  :  dominik.schrempf@gmail.com
+-- Stability   :  unstable
+-- Portability :  portable
+--
+-- Creation date: Mon Sep  7 13:36:45 2020.
+module ELynx.Tree.Parallel
+  ( parTree,
+    parBranchFoldMap,
+    parNodeFoldMap,
+  )
+where
+
+import Control.Parallel.Strategies
+import Data.Foldable
+import ELynx.Tree.Rooted
+
+myParList :: Strategy a -> Strategy [a]
+myParList _ [] = return []
+myParList s xs = do
+  ys <- parList s $ tail xs
+  y <- s $ head xs
+  return $ y : ys
+
+-- | Parallel evaluation strategy for a tree into normal form.
+--
+-- Evaluate the sub trees up to given layer in parallel.
+parTree :: (NFData e, NFData a) => Int -> Strategy (Tree e a)
+parTree n t@(Node br lb ts)
+  | n == 0 = rdeepseq t
+  | n == 1 = do
+    ts' <- myParList rdeepseq ts
+    return $ Node br lb ts'
+  | n >= 2 = do
+    ts' <- myParList (parTree (n -1)) ts
+    return $ Node br lb ts'
+  | otherwise = error "parTree: n is negative."
+
+branchFoldMap :: (e -> f) -> (f -> f -> f) -> Tree e a -> f
+branchFoldMap f op (Node br _ ts) = foldl' op (f br) $ map (branchFoldMap f op) ts
+
+-- | Map and fold over branches. Evaluate the sub trees up to given layer in parallel.
+parBranchFoldMap :: NFData f => Int -> (e -> f) -> (f -> f -> f) -> Tree e a -> f
+parBranchFoldMap 0 f op t = branchFoldMap f op t
+parBranchFoldMap n f op (Node br _ ts)
+  | n >= 1 = foldl' op (f br) (map (parBranchFoldMap (n - 1) f op) ts `using` myParList rdeepseq)
+  | otherwise = error "parBranchFoldMap: n is negative."
+
+nodeFoldMap :: (a -> b) -> (b -> b -> b) -> Tree e a -> b
+nodeFoldMap f op (Node _ lb ts) = foldl' op (f lb) $ map (nodeFoldMap f op) ts
+
+-- | Map and fold over nodes. Evaluate the sub trees up to given layer in parallel.
+parNodeFoldMap :: NFData b => Int -> (a -> b) -> (b -> b -> b) -> Tree e a -> b
+parNodeFoldMap 0 f op t = nodeFoldMap f op t
+parNodeFoldMap n f op (Node _ lb ts)
+  | n >= 1 = foldl' op (f lb) (map (parNodeFoldMap (n - 1) f op) ts `using` myParList rdeepseq)
+  | otherwise = error "parNodeFoldMap: n is negative."
diff --git a/src/ELynx/Tree/Partition.hs b/src/ELynx/Tree/Partition.hs
--- a/src/ELynx/Tree/Partition.hs
+++ b/src/ELynx/Tree/Partition.hs
@@ -15,10 +15,10 @@
 module ELynx.Tree.Partition
   ( -- * Data type
     Partition (fromPartition),
-    mp,
-    mpUnsafe,
-    bpToMp,
-    mpHuman,
+    pt,
+    ptUnsafe,
+    bpToPt,
+    ptHuman,
 
     -- * Work with 'Partition's
     partition,
@@ -33,15 +33,17 @@
 import ELynx.Tree.Bipartition
 import ELynx.Tree.Rooted
 
--- | Each branch of a tree partitions the leaves of the tree into two subsets
--- (see 'ELynx.Tree.Bipartition'). In a similar way, each internal node
--- (excluding the root node) partitions the leaves into three (or more) subsets
--- which is called 'Partition'. If the tree is multifurcating, and a
--- specific node has more than two children, the number of subsets induced by
--- this node is larger than three. Partitions are interesting in that we
--- can use them for calculating incompatible splits, see
--- 'ELynx.Tree.Distance'.
+-- | A partition of a tree is a grouping of the leaves of the tree into
+-- non-overlapping, non-empty sub sets.
 --
+-- For example, each branch of a tree partitions the leaves of the tree into two
+-- subsets (see 'ELynx.Tree.Bipartition'). In a similar way, each bifurcating
+-- internal node (excluding the root node) partitions the leaves into three
+-- subsets called a 'Partition'. If the tree is multifurcating, and a specific
+-- node has more than two children, the number of subsets induced by this node
+-- is larger than three. Partitions are interesting in that we can use them for
+-- calculating incompatible splits, see 'ELynx.Tree.Distance'.
+--
 -- The order of the subsets of a 'Partition' is meaningless. We ensure by
 -- construction that the subsets are ordered, and hence, that equality checks
 -- are meaningful.
@@ -52,21 +54,19 @@
 
 -- TODO: Check that list is not empty after filtering.
 
--- TODO: Rename these functions; don't use 'multi'.
-
 -- | Create a partition.
-mp :: Ord a => [Set a] -> Either String (Partition a)
-mp xs = case filter (not . S.null) xs of
+pt :: Ord a => [Set a] -> Either String (Partition a)
+pt xs = case filter (not . S.null) xs of
   [] -> Left "mp: Empty list."
-  xs' -> Right $ mpUnsafe xs'
+  xs' -> Right $ ptUnsafe xs'
 
 -- | Create a partition.
-mpUnsafe :: Ord a => [Set a] -> Partition a
-mpUnsafe xs = Partition (S.fromList xs)
+ptUnsafe :: Ord a => [Set a] -> Partition a
+ptUnsafe xs = Partition (S.fromList xs)
 
 -- | Convert a bipartition to a partition.
-bpToMp :: Ord a => Bipartition a -> Partition a
-bpToMp = mpUnsafe . tupleToList . fromBipartition
+bpToPt :: Ord a => Bipartition a -> Partition a
+bpToPt = ptUnsafe . tupleToList . fromBipartition
   where
     -- Be careful with tuples, because 'toList' does something very weird. It only
     -- takes the second element of the tuple!
@@ -76,8 +76,8 @@
 
 -- | Show a partition in a human readable form. Use a provided function to
 -- extract the valuable information.
-mpHuman :: Show a => Partition a -> String
-mpHuman (Partition xs) =
+ptHuman :: Show a => Partition a -> String
+ptHuman (Partition xs) =
   "(" ++ intercalate "|" (map setShow (S.toList xs)) ++ ")"
 
 -- Show the elements of a set in a human readable format.
@@ -93,7 +93,7 @@
 partition (Node _ _ []) = Left "partition: Encountered a leaf."
 partition t@(Node _ _ ts)
   | duplicateLeaves t = Left "partition: Tree contains duplicate leaves."
-  | otherwise = mp $ map (S.fromList . leaves) ts
+  | otherwise = pt $ map (S.fromList . leaves) ts
 
 -- | Get all 'Partition's of a tree.
 --
@@ -108,7 +108,7 @@
 partitions' _ (Node _ _ []) = S.empty
 partitions' p t@(Node _ _ ts) =
   S.unions $
-    either (const S.empty) S.singleton (mp (p : map label ts)) :
+    either (const S.empty) S.singleton (pt (p : map label ts)) :
     zipWith partitions' cs ts
   where
     cs = getComplementaryLeaves p t
diff --git a/src/ELynx/Tree/Rooted.hs b/src/ELynx/Tree/Rooted.hs
--- a/src/ELynx/Tree/Rooted.hs
+++ b/src/ELynx/Tree/Rooted.hs
@@ -97,7 +97,8 @@
 -- | Rooted rose trees with branch labels.
 --
 -- Unary instances such as 'Functor' act on node labels, and not on branch
--- labels. Binary instances such as 'Bifunctor' act on both labels.
+-- labels. Binary instances such as 'Bifunctor' act on both labels (`first` acts
+-- on branches, `second` on node labels).
 --
 -- Lifted instances are not provided.
 data Tree e a = Node
@@ -115,6 +116,7 @@
   fmap f ~(Node br lb ts) = Node br (f lb) $ map (fmap f) ts
   x <$ ~(Node br _ ts) = Node br x (map (x <$) ts)
 
+-- | The function 'first' acts on branch labels, 'second' on node labels.
 instance Bifunctor Tree where
   bimap f g ~(Node br lb ts) = Node (f br) (g lb) $ map (bimap f g) ts
   first f ~(Node br lb ts) = Node (f br) lb $ map (first f) ts
diff --git a/src/ELynx/Tree/Simulate/Coalescent.hs b/src/ELynx/Tree/Simulate/Coalescent.hs
--- a/src/ELynx/Tree/Simulate/Coalescent.hs
+++ b/src/ELynx/Tree/Simulate/Coalescent.hs
@@ -15,10 +15,10 @@
 where
 
 import Control.Monad.Primitive
+import ELynx.Tree.Distribution.CoalescentContinuous
 import ELynx.Tree.Measurable
 import ELynx.Tree.Phylogeny
 import ELynx.Tree.Rooted
-import ELynx.Tree.Distribution.CoalescentContinuous
 import Statistics.Distribution
 import System.Random.MWC
 
diff --git a/src/ELynx/Tree/Simulate/PointProcess.hs b/src/ELynx/Tree/Simulate/PointProcess.hs
--- a/src/ELynx/Tree/Simulate/PointProcess.hs
+++ b/src/ELynx/Tree/Simulate/PointProcess.hs
@@ -19,7 +19,7 @@
 -- and death process.
 module ELynx.Tree.Simulate.PointProcess
   ( PointProcess (..),
-    TimeSpec,
+    TimeSpec (..),
     simulate,
     toReconstructedTree,
     simulateReconstructedTree,
@@ -89,33 +89,52 @@
   }
   deriving (Read, Show, Eq)
 
--- | If nothing, sample time of origin from respective distribution. If time is
--- given, we need to know if we condition on the time of origin, or the time of
--- the most recent common ancestor (MRCA).
-type TimeSpec = Maybe (Time, Bool)
+-- | Tree height specification.
+data TimeSpec
+  = -- | Sample time of origin from respective distribution.
+    Random
+  | -- | Condition on time of origin.
+    Origin Time
+  | -- | Condition on time of most recent common ancestor (MRCA).
+    Mrca Time
 
 -- | Sample a point process using the 'BirthDeathDistribution'. The names of the
 -- points will be integers.
 simulate ::
   (PrimMonad m) =>
-  -- | Number of points (samples)
+  -- | Number of points (samples).
   Int ->
-  -- | Time of origin or MRCA
+  -- | Time of origin or MRCA.
   TimeSpec ->
-  -- | Birth rate
+  -- | Birth rate.
   Rate ->
-  -- | Death rate
+  -- | Death rate.
   Rate ->
-  -- | Generator (see 'System.Random.MWC')
+  -- | Generator.
   Gen (PrimState m) ->
   m (PointProcess Int Double)
+simulate n ts l m g
+  | n < 1 = error "Number of samples needs to be one or larger."
+  | l < 0.0 = error "Birth rate needs to be positive."
+  | otherwise = case ts of
+    Random -> simulateRandom n l m g
+    Origin t -> simulateOrigin n t l m g
+    Mrca t -> simulateMrca n t l m g
+
 -- No time of origin given. We also don't need to take care of the conditioning
 -- (origin or MRCA).
-simulate n Nothing l m g
+simulateRandom ::
+  PrimMonad m =>
+  Int ->
+  Double ->
+  Double ->
+  Gen (PrimState m) ->
+  m (PointProcess Int Double)
+simulateRandom n l m g
   | -- XXX. There is no formula for the over-critical process.
     m > l =
     error
-      "Time of origin distribution formula not available when mu > lambda. Please specify height for the moment."
+      "simulateRandom: Please specify height if mu > lambda."
   | -- For the critical process, we have no idea about the time of origin, but can
     -- use a specially derived distribution.
     m =~= l =
@@ -128,17 +147,24 @@
     abs (m - l) <= epsNearCriticalTimeOfOrigin =
     do
       t <- D.genContVar (TONCD n l m) g
-      simulate n (Just (t, False)) l m g
+      simulateOrigin n t l m g
   | -- For a sub-critical branching process, we can use the formula from Tanja Stadler.
     otherwise =
     do
       t <- D.genContVar (TOD n l m) g
-      simulate n (Just (t, False)) l m g
+      simulateOrigin n t l m g
+
 -- Time of origin is given.
-simulate n (Just (t, c)) l m g
-  | n < 1 = error "Number of samples needs to be one or larger."
-  | t < 0.0 = error "Time of origin needs to be positive."
-  | l < 0.0 = error "Birth rate needs to be positive."
+simulateOrigin ::
+  PrimMonad m =>
+  Int ->
+  Time ->
+  Double ->
+  Double ->
+  Gen (PrimState m) ->
+  m (PointProcess Int Double)
+simulateOrigin n t l m g
+  | t < 0.0 = error "simulateOrigin: Time of origin needs to be positive."
   | -- See Stadler, T., & Steel, M. (2019). Swapping birth and death: symmetries
     -- and transformations in phylodynamic models. , (), .
     -- http://dx.doi.org/10.1101/494583. Should be possible now.
@@ -147,28 +173,39 @@
     -- 1. The critical branching process.
     -- 2. The near critical branching process.
     -- 3. Normal values :).
-    (m =~= l) && not c = do
+    m =~= l = do
     !vs <- replicateM (n - 1) (D.genContVar (BDCD t l) g)
     return $ PointProcess [0 .. (n - 1)] vs t
-  | (abs (m - l) <= epsNearCriticalPointProcess) && not c = do
+  | abs (m - l) <= epsNearCriticalPointProcess = do
     !vs <- replicateM (n - 1) (D.genContVar (BDNCD t l m) g)
     return $ PointProcess [0 .. (n - 1)] vs t
-  | not c = do
+  | otherwise = do
     !vs <- replicateM (n - 1) (D.genContVar (BDD t l m) g)
     return $ PointProcess [0 .. (n - 1)] vs t
-  | (m =~= l) && c = do
+
+-- Time of Mrca is given.
+simulateMrca ::
+  PrimMonad m =>
+  Int ->
+  Time ->
+  Double ->
+  Double ->
+  Gen (PrimState m) ->
+  m (PointProcess Int Double)
+simulateMrca n t l m g
+  | t < 0.0 = error "simulateMrca: Time of MRCA needs to be positive."
+  | m =~= l = do
     !vs <- replicateM (n - 2) (D.genContVar (BDCD t l) g)
     vs' <- randomInsertList t vs g
     return $ PointProcess [0 .. (n - 1)] vs' t
-  | (abs (m - l) <= epsNearCriticalPointProcess) && c = do
+  | abs (m - l) <= epsNearCriticalPointProcess = do
     !vs <- replicateM (n - 2) (D.genContVar (BDNCD t l m) g)
     vs' <- randomInsertList t vs g
     return $ PointProcess [0 .. (n - 1)] vs' t
-  | c = do
+  | otherwise = do
     !vs <- replicateM (n - 2) (D.genContVar (BDD t l m) g)
     vs' <- randomInsertList t vs g
     return $ PointProcess [0 .. (n - 1)] vs' t
-  | otherwise = error "simulate: Fell through guard, this should never happen."
 
 -- Sort the values of a point process and their indices to be (the indices
 -- that they will have while creating the tree).
diff --git a/src/ELynx/Tree/Supported.hs b/src/ELynx/Tree/Supported.hs
--- a/src/ELynx/Tree/Supported.hs
+++ b/src/ELynx/Tree/Supported.hs
@@ -55,7 +55,7 @@
 
 -- A leaf has full support.
 highP :: Supported e => Double -> Tree e a -> Bool
-highP _  (Node _ _ []) = True
+highP _ (Node _ _ []) = True
 highP th (Node br _ _) = getSup br >= th
 
 -- See 'collapse'.
diff --git a/test/ELynx/Tree/BipartitionSpec.hs b/test/ELynx/Tree/BipartitionSpec.hs
--- a/test/ELynx/Tree/BipartitionSpec.hs
+++ b/test/ELynx/Tree/BipartitionSpec.hs
@@ -20,8 +20,8 @@
 import qualified Data.Map as M
 import Data.Set (Set)
 import qualified Data.Set as S
-import ELynx.Tree
 import ELynx.Tools
+import ELynx.Tree
 import Test.Hspec
 
 sfrom :: [BS.ByteString] -> S.Set BS.ByteString
diff --git a/test/ELynx/Tree/DistanceSpec.hs b/test/ELynx/Tree/DistanceSpec.hs
--- a/test/ELynx/Tree/DistanceSpec.hs
+++ b/test/ELynx/Tree/DistanceSpec.hs
@@ -17,9 +17,9 @@
 where
 
 import qualified Data.ByteString.Char8 as BS
+import ELynx.Tools
 import ELynx.Tree
 import ELynx.Tree.Arbitrary ()
-import ELynx.Tools
 import Test.Hspec
 import Test.QuickCheck
 
diff --git a/test/ELynx/Tree/Export/NewickSpec.hs b/test/ELynx/Tree/Export/NewickSpec.hs
--- a/test/ELynx/Tree/Export/NewickSpec.hs
+++ b/test/ELynx/Tree/Export/NewickSpec.hs
@@ -18,9 +18,9 @@
 import Data.Bifunctor
 import qualified Data.ByteString.Char8 as BS
 import qualified Data.ByteString.Lazy.Char8 as BL
+import ELynx.Tools
 import ELynx.Tree
 import ELynx.Tree.Arbitrary ()
-import ELynx.Tools
 import Test.Hspec
 import Test.QuickCheck hiding (label)
 
diff --git a/test/ELynx/Tree/Export/NexusSpec.hs b/test/ELynx/Tree/Export/NexusSpec.hs
--- a/test/ELynx/Tree/Export/NexusSpec.hs
+++ b/test/ELynx/Tree/Export/NexusSpec.hs
@@ -17,8 +17,8 @@
 where
 
 import qualified Data.ByteString.Char8 as BS
-import ELynx.Tree
 import ELynx.Tools
+import ELynx.Tree
 import Test.Hspec
 
 tree :: Tree Phylo BS.ByteString
diff --git a/test/ELynx/Tree/Import/NewickSpec.hs b/test/ELynx/Tree/Import/NewickSpec.hs
--- a/test/ELynx/Tree/Import/NewickSpec.hs
+++ b/test/ELynx/Tree/Import/NewickSpec.hs
@@ -18,8 +18,8 @@
 import Data.Attoparsec.ByteString.Char8
 import qualified Data.ByteString.Char8 as BS
 import Data.Either
-import ELynx.Tree
 import ELynx.Tools
+import ELynx.Tree
 import Test.Hspec
 
 -- sampleLabelByteString :: ByteString
diff --git a/test/ELynx/Tree/Import/NexusSpec.hs b/test/ELynx/Tree/Import/NexusSpec.hs
--- a/test/ELynx/Tree/Import/NexusSpec.hs
+++ b/test/ELynx/Tree/Import/NexusSpec.hs
@@ -17,8 +17,8 @@
 where
 
 import qualified Data.ByteString.Char8 as BS
-import ELynx.Tree
 import ELynx.Tools
+import ELynx.Tree
 import Test.Hspec
 
 file :: FilePath
diff --git a/test/ELynx/Tree/PartitionSpec.hs b/test/ELynx/Tree/PartitionSpec.hs
--- a/test/ELynx/Tree/PartitionSpec.hs
+++ b/test/ELynx/Tree/PartitionSpec.hs
@@ -26,8 +26,8 @@
 sol1 :: Set (Partition Int)
 sol1 =
   fromList
-    [ mpUnsafe [fromList [1], fromList [3], fromList [4, 5, 6]],
-      mpUnsafe [fromList [1, 3], fromList [4], fromList [5], fromList [6]]
+    [ ptUnsafe [fromList [1], fromList [3], fromList [4, 5, 6]],
+      ptUnsafe [fromList [1, 3], fromList [4], fromList [5], fromList [6]]
     ]
 
 ex2 :: Tree () Int
@@ -36,8 +36,8 @@
 sol2 :: Set (Partition Int)
 sol2 =
   fromList
-    [ mpUnsafe [fromList [1], fromList [2], fromList [3, 4], fromList [5]],
-      mpUnsafe [fromList [1, 2, 5], fromList [3], fromList [4]]
+    [ ptUnsafe [fromList [1], fromList [2], fromList [3, 4], fromList [5]],
+      ptUnsafe [fromList [1, 2, 5], fromList [3], fromList [4]]
     ]
 
 spec :: Spec
diff --git a/test/ELynx/Tree/RootedSpec.hs b/test/ELynx/Tree/RootedSpec.hs
--- a/test/ELynx/Tree/RootedSpec.hs
+++ b/test/ELynx/Tree/RootedSpec.hs
@@ -19,8 +19,8 @@
 import qualified Data.ByteString.Char8 as BS
 import qualified Data.ByteString.Lazy.Char8 as BL
 import Data.Maybe
-import ELynx.Tree
 import ELynx.Tools
+import ELynx.Tree
 import Test.Hspec
 
 node :: Int -> Tree () Int
diff --git a/test/ELynx/Tree/SupportedSpec.hs b/test/ELynx/Tree/SupportedSpec.hs
--- a/test/ELynx/Tree/SupportedSpec.hs
+++ b/test/ELynx/Tree/SupportedSpec.hs
@@ -17,8 +17,8 @@
 where
 
 import qualified Data.ByteString.Char8 as BS
-import ELynx.Tree
 import ELynx.Tools
+import ELynx.Tree
 import Test.Hspec
 
 collapseTree :: Tree Phylo BS.ByteString