elynx-tree 0.4.0 → 0.4.1

raw patch · 26 files changed

+336/−151 lines, 26 filesdep +parallelPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: parallel

API changes (from Hackage documentation)

Files

• ChangeLog.md +14/−0
• README.md +32/−29
• bench/Bench.hs +50/−4
• data/UltraMetric.log +0/−22
• elynx-tree.cabal +8/−3
• src/ELynx/Tree.hs +6/−0
• src/ELynx/Tree/Bipartition.hs +6/−4
• src/ELynx/Tree/Distance.hs +7/−9
• src/ELynx/Tree/Export/Newick.hs +8/−2
• src/ELynx/Tree/Import/Newick.hs +34/−8
• src/ELynx/Tree/Import/Nexus.hs +5/−5
• src/ELynx/Tree/Parallel.hs +62/−0
• src/ELynx/Tree/Partition.hs +25/−25
• src/ELynx/Tree/Rooted.hs +3/−1
• src/ELynx/Tree/Simulate/Coalescent.hs +1/−1
• src/ELynx/Tree/Simulate/PointProcess.hs +62/−25
• src/ELynx/Tree/Supported.hs +1/−1
• test/ELynx/Tree/BipartitionSpec.hs +1/−1
• test/ELynx/Tree/DistanceSpec.hs +1/−1
• test/ELynx/Tree/Export/NewickSpec.hs +1/−1
• test/ELynx/Tree/Export/NexusSpec.hs +1/−1
• test/ELynx/Tree/Import/NewickSpec.hs +1/−1
• test/ELynx/Tree/Import/NexusSpec.hs +1/−1
• test/ELynx/Tree/PartitionSpec.hs +4/−4
• test/ELynx/Tree/RootedSpec.hs +1/−1
• test/ELynx/Tree/SupportedSpec.hs +1/−1

@@ -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

@@ -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

@@ -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+        ]+    ]

@@ -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.-------------------------------------------------

@@ -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

@@ -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

@@ -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'.

@@ -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. --

@@ -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

@@ -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.

@@ -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)

@@ -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."

@@ -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

@@ -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

@@ -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 

@@ -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).

@@ -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'.

@@ -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

@@ -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 

@@ -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) 

@@ -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

@@ -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

@@ -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

@@ -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

@@ -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

@@ -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