elynx-tree 0.3.4 → 0.4.0

raw patch · 82 files changed

+4556/−4554 lines, 82 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

• ChangeLog.md +6/−0
• README.md +27/−27
• bench/Bench.hs +2/−2
• elynx-tree.cabal +40/−40
• src/ELynx/Data/Topology/Phylogeny.hs +0/−115
• src/ELynx/Data/Topology/Rooted.hs +0/−207
• src/ELynx/Data/Tree.hs +0/−72
• src/ELynx/Data/Tree/Bipartition.hs +0/−198
• src/ELynx/Data/Tree/Distance.hs +0/−139
• src/ELynx/Data/Tree/Measurable.hs +0/−124
• src/ELynx/Data/Tree/Named.hs +0/−43
• src/ELynx/Data/Tree/Partition.hs +0/−137
• src/ELynx/Data/Tree/Phylogeny.hs +0/−487
• src/ELynx/Data/Tree/Rooted.hs +0/−355
• src/ELynx/Data/Tree/Splittable.hs +0/−29
• src/ELynx/Data/Tree/Supported.hs +0/−66
• src/ELynx/Data/Tree/Zipper.hs +0/−160
• src/ELynx/Distribution/BirthDeath.hs +0/−95
• src/ELynx/Distribution/BirthDeathCritical.hs +0/−84
• src/ELynx/Distribution/BirthDeathCriticalNoTime.hs +0/−76
• src/ELynx/Distribution/BirthDeathNearlyCritical.hs +0/−109
• src/ELynx/Distribution/CoalescentContinuous.hs +0/−26
• src/ELynx/Distribution/TimeOfOrigin.hs +0/−95
• src/ELynx/Distribution/TimeOfOriginNearCritical.hs +0/−96
• src/ELynx/Distribution/Types.hs +0/−22
• src/ELynx/Export/Tree/Newick.hs +0/−53
• src/ELynx/Export/Tree/Nexus.hs +0/−31
• src/ELynx/Import/Tree/Newick.hs +0/−235
• src/ELynx/Import/Tree/Nexus.hs +0/−45
• src/ELynx/Simulate/Coalescent.hs +0/−60
• src/ELynx/Simulate/PointProcess.hs +0/−289
• src/ELynx/Topology/Phylogeny.hs +115/−0
• src/ELynx/Topology/Rooted.hs +207/−0
• src/ELynx/Tree.hs +82/−0
• src/ELynx/Tree/Bipartition.hs +198/−0
• src/ELynx/Tree/Distance.hs +139/−0
• src/ELynx/Tree/Distribution/BirthDeath.hs +95/−0
• src/ELynx/Tree/Distribution/BirthDeathCritical.hs +84/−0
• src/ELynx/Tree/Distribution/BirthDeathCriticalNoTime.hs +76/−0
• src/ELynx/Tree/Distribution/BirthDeathNearlyCritical.hs +109/−0
• src/ELynx/Tree/Distribution/CoalescentContinuous.hs +26/−0
• src/ELynx/Tree/Distribution/TimeOfOrigin.hs +95/−0
• src/ELynx/Tree/Distribution/TimeOfOriginNearCritical.hs +96/−0
• src/ELynx/Tree/Distribution/Types.hs +22/−0
• src/ELynx/Tree/Export/Newick.hs +53/−0
• src/ELynx/Tree/Export/Nexus.hs +31/−0
• src/ELynx/Tree/Import/Newick.hs +235/−0
• src/ELynx/Tree/Import/Nexus.hs +45/−0
• src/ELynx/Tree/Measurable.hs +124/−0
• src/ELynx/Tree/Named.hs +43/−0
• src/ELynx/Tree/Partition.hs +137/−0
• src/ELynx/Tree/Phylogeny.hs +487/−0
• src/ELynx/Tree/Rooted.hs +355/−0
• src/ELynx/Tree/Simulate/Coalescent.hs +60/−0
• src/ELynx/Tree/Simulate/PointProcess.hs +288/−0
• src/ELynx/Tree/Splittable.hs +29/−0
• src/ELynx/Tree/Supported.hs +66/−0
• src/ELynx/Tree/Zipper.hs +160/−0
• test/ELynx/Data/Topology/RootedSpec.hs +0/−23
• test/ELynx/Data/Tree/Arbitrary.hs +0/−66
• test/ELynx/Data/Tree/BipartitionSpec.hs +0/−88
• test/ELynx/Data/Tree/DistanceSpec.hs +0/−259
• test/ELynx/Data/Tree/PartitionSpec.hs +0/−53
• test/ELynx/Data/Tree/PhylogenySpec.hs +0/−155
• test/ELynx/Data/Tree/RootedSpec.hs +0/−68
• test/ELynx/Data/Tree/SupportedSpec.hs +0/−42
• test/ELynx/Export/Tree/NewickSpec.hs +0/−52
• test/ELynx/Export/Tree/NexusSpec.hs +0/−46
• test/ELynx/Import/Tree/NewickSpec.hs +0/−134
• test/ELynx/Import/Tree/NexusSpec.hs +0/−51
• test/ELynx/Topology/RootedSpec.hs +23/−0
• test/ELynx/Tree/Arbitrary.hs +66/−0
• test/ELynx/Tree/BipartitionSpec.hs +87/−0
• test/ELynx/Tree/DistanceSpec.hs +258/−0
• test/ELynx/Tree/Export/NewickSpec.hs +50/−0
• test/ELynx/Tree/Export/NexusSpec.hs +43/−0
• test/ELynx/Tree/Import/NewickSpec.hs +133/−0
• test/ELynx/Tree/Import/NexusSpec.hs +49/−0
• test/ELynx/Tree/PartitionSpec.hs +53/−0
• test/ELynx/Tree/PhylogenySpec.hs +155/−0
• test/ELynx/Tree/RootedSpec.hs +67/−0
• test/ELynx/Tree/SupportedSpec.hs +40/−0

@@ -5,6 +5,12 @@ ## Unreleased changes  +## Version 0.4.0++-   Major refactor of `elynx-tree`. All required function can now conveniently+    reexported by `ELynx.Tree`.++ ## Version 0.3.4  -   Improve `slynx examine`; show hamming distance; show constant sites.

@@ -2,7 +2,7 @@  # The ELynx Suite -Version: 0.3.4.+Version: 0.4.0. 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.4.+    ELynx Suite version 0.4.0.     Developed by Dominik Schrempf.-    Compiled on August 21, 2020, at 09:40 am, UTC.+    Compiled on September 4, 2020, at 13:37 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.4.+    ELynx Suite version 0.4.0.     Developed by Dominik Schrempf.-    Compiled on August 21, 2020, at 09:40 am, UTC.+    Compiled on September 4, 2020, at 13:37 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.4.+    ELynx Suite version 0.4.0.     Developed by Dominik Schrempf.-    Compiled on August 21, 2020, at 09:40 am, UTC.+    Compiled on September 4, 2020, at 13:37 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.4.+    ELynx Suite version 0.4.0.     Developed by Dominik Schrempf.-    Compiled on August 21, 2020, at 09:40 am, UTC.+    Compiled on September 4, 2020, at 13:37 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.4.+    ELynx Suite version 0.4.0.     Developed by Dominik Schrempf.-    Compiled on August 21, 2020, at 09:40 am, UTC.+    Compiled on September 4, 2020, at 13:37 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.4.+    ELynx Suite version 0.4.0.     Developed by Dominik Schrempf.-    Compiled on August 21, 2020, at 09:40 am, UTC.+    Compiled on September 4, 2020, at 13:37 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.4.+    ELynx Suite version 0.4.0.     Developed by Dominik Schrempf.-    Compiled on August 21, 2020, at 09:40 am, UTC.+    Compiled on September 4, 2020, at 13:37 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.4.+    ELynx Suite version 0.4.0.     Developed by Dominik Schrempf.-    Compiled on August 21, 2020, at 09:40 am, UTC.+    Compiled on September 4, 2020, at 13:37 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.4.+    ELynx Suite version 0.4.0.     Developed by Dominik Schrempf.-    Compiled on August 21, 2020, at 09:40 am, UTC.+    Compiled on September 4, 2020, at 13:37 pm, UTC.          Usage: tlynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME]                   [-f|--force] COMMAND@@ -398,9 +398,9 @@      tlynx compare --help -    ELynx Suite version 0.3.4.+    ELynx Suite version 0.4.0.     Developed by Dominik Schrempf.-    Compiled on August 21, 2020, at 09:40 am, UTC.+    Compiled on September 4, 2020, at 13:37 pm, UTC.          Usage: tlynx compare [-n|--normalize] [-b|--bipartitions] [-t|--intersect]                           [-f|--newick-format FORMAT] NAMES@@ -427,9 +427,9 @@      tlynx examine --help -    ELynx Suite version 0.3.4.+    ELynx Suite version 0.4.0.     Developed by Dominik Schrempf.-    Compiled on August 21, 2020, at 09:40 am, UTC.+    Compiled on September 4, 2020, at 13:37 pm, UTC.          Usage: tlynx examine INPUT-FILE [-f|--newick-format FORMAT]       Compute summary statistics of phylogenetic trees.@@ -451,9 +451,9 @@      tlynx simulate --help -    ELynx Suite version 0.3.4.+    ELynx Suite version 0.4.0.     Developed by Dominik Schrempf.-    Compiled on August 21, 2020, at 09:40 am, UTC.+    Compiled on September 4, 2020, at 13:37 pm, UTC.          Usage: tlynx simulate (-t|--nTrees INT) (-n|--nLeaves INT) PROCESS                            [-u|--sub-sample DOUBLE] [-s|--summary-statistics] @@ -488,9 +488,9 @@      elynx --help | head -n -16 -    ELynx Suite version 0.3.4.+    ELynx Suite version 0.4.0.     Developed by Dominik Schrempf.-    Compiled on August 21, 2020, at 09:40 am, UTC.+    Compiled on September 4, 2020, at 13:37 pm, UTC.          Usage: elynx COMMAND       Validate and redo past ELynx analyses

@@ -13,9 +13,9 @@  import Criterion.Main import qualified Data.ByteString.Char8 as BS-import ELynx.Data.Tree-import ELynx.Import.Tree.Newick import ELynx.Tools+import ELynx.Tree+import ELynx.Tree.Import.Newick  treeFileMany :: FilePath treeFileMany = "data/Many.trees"

@@ -1,6 +1,6 @@ cabal-version:  2.2 name:           elynx-tree-version:        0.3.4+version:        0.4.0 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@@ -40,33 +40,33 @@  library   exposed-modules:-      ELynx.Data.Topology.Phylogeny-      ELynx.Data.Topology.Rooted-      ELynx.Data.Tree-      ELynx.Data.Tree.Bipartition-      ELynx.Data.Tree.Distance-      ELynx.Data.Tree.Measurable-      ELynx.Data.Tree.Named-      ELynx.Data.Tree.Partition-      ELynx.Data.Tree.Phylogeny-      ELynx.Data.Tree.Rooted-      ELynx.Data.Tree.Splittable-      ELynx.Data.Tree.Supported-      ELynx.Data.Tree.Zipper-      ELynx.Distribution.BirthDeath-      ELynx.Distribution.BirthDeathCritical-      ELynx.Distribution.BirthDeathCriticalNoTime-      ELynx.Distribution.BirthDeathNearlyCritical-      ELynx.Distribution.CoalescentContinuous-      ELynx.Distribution.TimeOfOrigin-      ELynx.Distribution.TimeOfOriginNearCritical-      ELynx.Distribution.Types-      ELynx.Export.Tree.Newick-      ELynx.Export.Tree.Nexus-      ELynx.Import.Tree.Newick-      ELynx.Import.Tree.Nexus-      ELynx.Simulate.Coalescent-      ELynx.Simulate.PointProcess+      ELynx.Topology.Phylogeny+      ELynx.Topology.Rooted+      ELynx.Tree+      ELynx.Tree.Bipartition+      ELynx.Tree.Distance+      ELynx.Tree.Measurable+      ELynx.Tree.Named+      ELynx.Tree.Partition+      ELynx.Tree.Phylogeny+      ELynx.Tree.Rooted+      ELynx.Tree.Splittable+      ELynx.Tree.Supported+      ELynx.Tree.Zipper+      ELynx.Tree.Distribution.BirthDeath+      ELynx.Tree.Distribution.BirthDeathCritical+      ELynx.Tree.Distribution.BirthDeathCriticalNoTime+      ELynx.Tree.Distribution.BirthDeathNearlyCritical+      ELynx.Tree.Distribution.CoalescentContinuous+      ELynx.Tree.Distribution.TimeOfOrigin+      ELynx.Tree.Distribution.TimeOfOriginNearCritical+      ELynx.Tree.Distribution.Types+      ELynx.Tree.Export.Newick+      ELynx.Tree.Export.Nexus+      ELynx.Tree.Import.Newick+      ELynx.Tree.Import.Nexus+      ELynx.Tree.Simulate.Coalescent+      ELynx.Tree.Simulate.PointProcess   other-modules:       Paths_elynx_tree   autogen-modules:@@ -94,18 +94,18 @@   type: exitcode-stdio-1.0   main-is: Spec.hs   other-modules:-      ELynx.Data.Topology.RootedSpec-      ELynx.Data.Tree.Arbitrary-      ELynx.Data.Tree.BipartitionSpec-      ELynx.Data.Tree.DistanceSpec-      ELynx.Data.Tree.PartitionSpec-      ELynx.Data.Tree.PhylogenySpec-      ELynx.Data.Tree.RootedSpec-      ELynx.Data.Tree.SupportedSpec-      ELynx.Export.Tree.NewickSpec-      ELynx.Export.Tree.NexusSpec-      ELynx.Import.Tree.NewickSpec-      ELynx.Import.Tree.NexusSpec+      ELynx.Topology.RootedSpec+      ELynx.Tree.Arbitrary+      ELynx.Tree.BipartitionSpec+      ELynx.Tree.DistanceSpec+      ELynx.Tree.PartitionSpec+      ELynx.Tree.PhylogenySpec+      ELynx.Tree.RootedSpec+      ELynx.Tree.SupportedSpec+      ELynx.Tree.Export.NewickSpec+      ELynx.Tree.Export.NexusSpec+      ELynx.Tree.Import.NewickSpec+      ELynx.Tree.Import.NexusSpec       Paths_elynx_tree   hs-source-dirs:       test

@@ -1,115 +0,0 @@--- |--- Module      :  ELynx.Data.Topology.Phylogeny--- Description :  Phylogenetic topologies--- Copyright   :  (c) Dominik Schrempf, 2020--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Sat Jul 18 13:15:49 2020.------ A topology, as it is used in phylogenetics is a 'Topology' with unique leaf--- labels, and the order of the topologies in the sub-forest is considered to be--- meaningless.------ Internally, however, the underlying 'Topology' data structure stores the--- sub-forest as a (non-empty) list, which has a specific order. Hence, we have--- to do some tricks when comparing topologies, and topology comparison is slow.------ Also, the uniqueness of the leaves is not ensured by the data type, but has--- to be checked at runtime. Functions relying on the tree to have unique leaves--- do perform this check, and return 'Left' with an error message, if the tree--- has duplicate leaves.------ Note: Topologies are rooted.------ Note: Topologies encoded in Newick format correspond to rooted topologies. By--- convention only, a topology parsed from Newick format is usually thought to--- be unrooted, when the root node is multifurcating and has three children.--- This convention is not enforced here. Newick topologies are just parsed as--- they are, and a rooted topology is returned.------ The bifurcating root of a topology can be changed with 'roots', or 'rootAt'.------ Topologies with multifurcating root nodes can be properly rooted using--- 'outgroup'.-module ELynx.Data.Topology.Phylogeny-  ( outgroup,-    roots,-    rootAt,-  )-where--import Data.Set (Set)-import ELynx.Data.Topology.Rooted-import ELynx.Data.Tree.Bipartition---- TODO.---- -- | Remove multifurcations.--- ----- -- A caterpillar like bifurcating tree is used to resolve all multifurcations on--- -- a tree. The multifurcating nodes are copied.--- ----- -- Branch labels are not handled.--- resolve :: Tree () a -> Tree () a--- resolve t@(Node _ _ []) = t--- resolve (Node _ l [x]) = Node () l [resolve x]--- resolve (Node _ l [x, y]) = Node () l $ map resolve [x, y]--- resolve (Node _ l (x : xs)) = Node () l $ map resolve [x, Node () l xs]---- | Resolve a multifurcation at the root using an outgroup.-outgroup :: Ord a => Set a -> Topology a -> Either String (Topology a)-outgroup = undefined---- -- | For a rooted tree with a bifurcating root node, get all possible rooted--- -- trees.--- ----- -- The root node is moved.--- ----- -- For a tree with @l=2@ leaves, there is one rooted tree. For a bifurcating--- -- tree with @l>2@ leaves, there are @(2l-3)@ rooted trees. For a general tree--- -- with a bifurcating root node, and a total number of @n>2@ nodes, there are--- -- (n-2) rooted trees.--- ----- -- Moving a multifurcating root node to another branch would change the--- -- topology, and so, a bifurcating root is required. To resolve a multifurcating--- -- root, please see and use TODO.--- ----- -- Branch labels are not handled, but see 'rootsBranch'.--- ----- -- 'rootAt' roots the tree at a specific position.--- ----- -- Return 'Left' if the root node is not 'bifurcating'.--- roots :: Tree () a -> Either String (Forest () a)--- roots (Node _ _ []) = Left "roots: Root node is a leaf."--- roots (Node _ _ [_]) = Left "roots: Root node has degree two."--- roots t@(Node _ c [tL, tR]) = Right $ t : descend id () c tR tL ++ descend id () c tL tR--- roots _ = Left "roots: Root node is multifurcating."---- | For a rooted topology with a bifurcating root node, get all possible rooted--- topologies.-roots :: Topology a -> Either String (Forest a)-roots = undefined---- -- | Root a tree at a specific position.--- ----- -- Root the tree at the branch defined by the given bipartition. The original--- -- root node is moved to the new position.--- ----- -- The root node must be bifurcating (see 'roots').--- ----- -- Branch labels are not handled, but see 'rootAtBranch'.--- ----- -- Return 'Left', if:--- -- - the root node is not bifurcating;--- -- - the tree has duplicate leaves;--- -- - the bipartition does not match the leaves of the tree.--- rootAt :: Ord a => Bipartition a -> Tree () a -> Either String (Tree () a)--- rootAt = rootAtBranch id---- | Root a tree at a specific position.-rootAt :: Ord a => Bipartition a -> Topology a -> Either String (Forest a)-rootAt = undefined

@@ -1,207 +0,0 @@-{-# LANGUAGE DeriveDataTypeable #-}-{-# LANGUAGE DeriveGeneric #-}---- |--- Module      :  ELynx.Data.Topology.Rooted--- Description :  Topologies--- Copyright   :  (c) Dominik Schrempf, 2020--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Sat Jul 11 10:28:28 2020.------ A 'Topology' differs from a classical rose 'Data.Tree.Tree' in that it does--- not have internal node labels. The leaves have labels.------ For rooted trees, please see 'ELynx.Data.Tree.Rooted'.------ In phylogenetics, the order of children of a topology node is arbitrary.--- Internally, however, the underlying 'Topology' data structure stores the--- sub-forest as a (non-empty) list, which has a specific order. Hence, we have--- to do some tricks when comparing topologies, and topology comparison is slow.-module ELynx.Data.Topology.Rooted-  ( -- * Data type-    Topology (..),-    Forest,-    fromTree,-    fromLabeledTree,--    -- * Functions-    degree,-    leaves,-    flatten,-    identify,-    prune,-    dropLeavesWith,-    zipTreesWith,-    zipTrees,-    duplicateLeaves,-  )-where--import Control.Applicative-import Control.DeepSeq-import Control.Monad-import Data.Aeson-import Data.Data-import Data.Foldable-import Data.List.NonEmpty (NonEmpty)-import qualified Data.List.NonEmpty as N-import Data.Maybe-import qualified Data.Set as S-import Data.Traversable-import qualified Data.Tree as T-import qualified ELynx.Data.Tree.Rooted as R-import GHC.Generics--singleton :: NonEmpty a -> Bool-singleton xs = 1 == length (N.take 2 xs)---- | Rooted topologies with leaf labels.-data Topology a-  = Node {forest :: Forest a}-  | Leaf {label :: a}-  deriving (Eq, Read, Show, Data, Generic)---- | A shortcut.-type Forest a = NonEmpty (Topology a)--instance Functor Topology where-  fmap f (Node ts) = Node $ fmap (fmap f) ts-  fmap f (Leaf lb) = Leaf $ f lb--instance Foldable Topology where-  foldMap f (Node ts) = foldMap (foldMap f) ts-  foldMap f (Leaf lb) = f lb--  null _ = False-  {-# INLINE null #-}--  toList = flatten-  {-# INLINE toList #-}--instance Traversable Topology where-  traverse g (Node ts) = Node <$> traverse (traverse g) ts-  traverse g (Leaf lb) = Leaf <$> g lb---- TODO: This type checks, but I doubt the implementation is bug-free.-instance Applicative Topology where-  pure = Leaf--  (Node tsF) <*> tx = Node $ fmap (<*> tx) tsF-  (Leaf lbF) <*> tx = lbF <$> tx--  liftA2 f (Node tsX) ty = Node $ fmap (\tx -> liftA2 f tx ty) tsX-  liftA2 f (Leaf lbX) (Node tsY) = Node $ fmap (f lbX <$>) tsY-  liftA2 f (Leaf lbX) (Leaf lbY) = Leaf $ f lbX lbY--  (Node tsX) *> ty@(Node tsY) = Node $ tsY <> fmap (*> ty) tsX-  (Leaf _) *> (Node tsY) = Node tsY-  _ *> (Leaf lbY) = Leaf lbY--  (Node tsX) <* ty = Node $ fmap (<* ty) tsX-  (Leaf lbX) <* _ = Leaf lbX---- TODO: This type checks, but I doubt the implementation is bug-free.-instance Monad Topology where-  (Node ts) >>= f = Node $ fmap (>>= f) ts-  (Leaf lb) >>= f = case f lb of-    Node ts' -> Node ts'-    Leaf lb' -> Leaf lb'--instance NFData a => NFData (Topology a) where-  rnf (Node ts) = rnf ts-  rnf (Leaf lb) = rnf lb--instance ToJSON a => ToJSON (Topology a)--instance FromJSON a => FromJSON (Topology a)---- | The degree of the root node.-degree :: Topology a -> Int-degree (Node ts) = (+ 1) $ length ts-degree (Leaf _) = 1---- | Set of leaves.-leaves :: Ord a => Topology a -> [a]-leaves (Leaf lb) = [lb]-leaves (Node ts) = concatMap leaves ts---- | Return leaf labels in pre-order.-flatten :: Topology a -> [a]-flatten t = squish t []-  where-    squish (Node ts) xs = foldr squish xs ts-    squish (Leaf lb) xs = lb : xs---- TODO: Provide and fix tests, provide arbitrary instances.---- | Convert a rooted rose tree to a rooted topology. Internal node labels are lost.-fromTree :: T.Tree a -> Topology a-fromTree (T.Node lb []) = Leaf lb-fromTree (T.Node _ xs) = Node $ fromTree <$> N.fromList xs---- | Convert a rooted, labeled rose tree to a rooted topology. Branch labels and--- internal node labels are lost.-fromLabeledTree :: R.Tree e a -> Topology a-fromLabeledTree (R.Node _ lb []) = Leaf lb-fromLabeledTree (R.Node _ _ xs) = Node $ fromLabeledTree <$> N.fromList xs---- | Label the leaves with unique integers starting at 0.-identify :: Traversable t => t a -> t Int-identify = snd . mapAccumL (\i _ -> (i + 1, i)) (0 :: Int)---- | Prune degree two nodes.-prune :: Topology a -> Topology a-prune (Node ts)-  | singleton ts = Node $ fmap prune $ forest $ N.head ts-  | otherwise = Node $ fmap prune ts-prune (Leaf lb) = Leaf lb---- | Drop leaves satisfying predicate.------ Degree two nodes may arise.------ Return 'Nothing' if all leaves satisfy the predicate.-dropLeavesWith :: (a -> Bool) -> Topology a -> Maybe (Topology a)-dropLeavesWith p (Leaf lb)-  | p lb = Nothing-  | otherwise = Just $ Leaf lb-dropLeavesWith p (Node ts) =-  if null ts'-    then Nothing-    else -- XXX: May be slow, unnecessary conversion to and from list.-      Just $ Node $ N.fromList ts'-  where-    ts' = catMaybes $ N.toList $ fmap (dropLeavesWith p) ts---- | Zip leaves of two equal topologies.------ Return 'Nothing' if the topologies are different.-zipTreesWith :: (a1 -> a2 -> a) -> Topology a1 -> Topology a2 -> Maybe (Topology a)-zipTreesWith f (Node tsL) (Node tsR) =-  if N.length tsL == N.length tsR-    then -- XXX: May be slow, unnecessary conversion to and from list.-      zipWithM (zipTreesWith f) (N.toList tsL) (N.toList tsR) >>= Just . Node . N.fromList-    else Nothing-zipTreesWith f (Leaf lbL) (Leaf lbR) = Just $ Leaf $ f lbL lbR-zipTreesWith _ _ _ = Nothing---- | Zip leaves of two equal topologies.------ Return 'Nothing' if the topologies are different.-zipTrees :: Topology a1 -> Topology a2 -> Maybe (Topology (a1, a2))-zipTrees = zipTreesWith (,)--duplicates :: Ord a => [a] -> Bool-duplicates = go S.empty-  where-    go _ [] = False-    go seen (x : xs) = x `S.member` seen || go (S.insert x seen) xs---- | Check if a topology has duplicate leaves.-duplicateLeaves :: Ord a => Topology a -> Bool-duplicateLeaves = duplicates . leaves

@@ -1,72 +0,0 @@--- TODO: Topology data type.--- data Topology a = Node (NonEmptySet (Topology a)) | Leaf a---- |--- Module      :  ELynx.Data.Tree--- Description :  Phylogenetic trees--- Copyright   :  (c) Dominik Schrempf 2020--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Sat Mar 21 16:27:20 2020.-module ELynx.Data.Tree-  ( -- * Rooted trees-    module ELynx.Data.Tree.Rooted,-    module ELynx.Data.Tree.Zipper,--    -- * Branch label classes-    module ELynx.Data.Tree.Measurable,-    module ELynx.Data.Tree.Splittable,-    module ELynx.Data.Tree.Supported,--    -- * Node label classes-    module ELynx.Data.Tree.Named,--    -- * Phylogenies-    module ELynx.Data.Tree.Phylogeny,--    -- * Partitions and distances-    module ELynx.Data.Tree.Bipartition,-    module ELynx.Data.Tree.Partition,-    module ELynx.Data.Tree.Distance,-  )-where--import ELynx.Data.Tree.Bipartition-import ELynx.Data.Tree.Distance-import ELynx.Data.Tree.Measurable-import ELynx.Data.Tree.Named-import ELynx.Data.Tree.Partition-import ELynx.Data.Tree.Phylogeny-import ELynx.Data.Tree.Rooted-import ELynx.Data.Tree.Splittable-import ELynx.Data.Tree.Supported-import ELynx.Data.Tree.Zipper---- -- | An evolutionary label has some information about where the corresponding--- -- node is on the tree, and if the node is 'extant', 'extinct', 'internal', or--- -- 'external'. The latter two could also be determined from the tree. This could--- -- be species, genes or individuals; probably more.--- class EvoLabel n where---   extant          :: n -> Bool---   extinct         :: n -> Bool----   internal        :: n -> Bool---   internal n = not $ extant n || extinct n----   external        :: n -> Bool---   external   = not . internal---- -- -- | Glue branches together, so that one new tree emerges. It's root node is--- -- -- new, the sub-forest has to be given (a list of trees).--- -- glue :: (NodeType c)--- --      => PhyloLabel a b c       -- ^ New root node.--- --      -> [PhyloTree a b c]      -- ^ Sub-forest.--- --      -> PhyloTree a b c--- -- glue s@(PhyloLabel _ _ n) ts--- --   | extant n  = error "Root node cannot be of type 'Exant'."--- --   | extinct n = error "Root node cannot be of type 'Extinct'."--- --   | otherwise = Node s ts

@@ -1,198 +0,0 @@-{-# LANGUAGE GeneralizedNewtypeDeriving #-}---- |--- Module      :  ELynx.Data.Tree.Bipartition--- Description :  Bipartitions on trees--- Copyright   :  (c) Dominik Schrempf 2020--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Fri Aug 30 15:28:17 2019.------ 'Bipartition's are weird in that--- > Bipartition x y == Bipartition y x--- is True.------ Also,--- > Bipartition x y > Bipartition y x--- is False, even when @x > y@.------ That's why we have to make sure that for--- > Bipartition x y--- we always have @x >= y@.-module ELynx.Data.Tree.Bipartition-  ( groups,--    -- * Data type-    Bipartition (fromBipartition),-    bp,-    bpUnsafe,-    toSet,-    bpHuman,--    -- * Work with 'Bipartition's-    bipartition,-    bipartitions,-    getComplementaryLeaves,-    bipartitionToBranch,-  )-where--import Control.Comonad-import Control.DeepSeq-import Data.List hiding (partition)-import Data.Map (Map)-import qualified Data.Map as M-import Data.Set (Set)-import qualified Data.Set as S-import ELynx.Data.Tree.Rooted---- | Each node of a tree is root of an induced subtree. Set the node labels to--- the leaves of the induced subtrees.-groups :: Tree e a -> Tree e [a]--- 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.------ 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.-newtype Bipartition a = Bipartition-  { fromBipartition :: (Set a, Set a)-  }-  deriving (Eq, Ord, Show, Read, NFData)---- | Create a bipartition from two sets.------ Ensure that the smaller set comes first.------ Return 'Left' if one set is empty.-bp :: Ord a => Set a -> Set a -> Either String (Bipartition a)-bp xs ys-  | S.null xs = Left "bp: Left set empty."-  | S.null ys = Left "bp: Right set empty."-  | otherwise = Right $ bpUnsafe xs ys---- | Create a bipartition from two sets.------ Ensure that the smaller set comes first.-bpUnsafe :: Ord a => Set a -> Set a -> Bipartition a-bpUnsafe xs ys = if xs >= ys then Bipartition (xs, ys) else Bipartition (ys, xs)---- | Conversion to a set containing both partitions.-toSet :: Ord a => Bipartition a -> Set a-toSet (Bipartition (x, y)) = S.union x y---- I decided not to provide a human readable show instance because I need the--- following identity to hold:------ > read . show = id------ This identity is met by the derived instance anyways. A more human readable--- instance would most likely violate the identity.---- | Show a bipartition in a human readable format. Use a provided function to--- extract information of interest.-bpHuman :: Show a => Bipartition a -> String-bpHuman (Bipartition (x, y)) = "(" ++ setShow x ++ "|" ++ setShow y ++ ")"---- Show the elements of a set in a human readable format.-setShow :: Show a => Set a -> String-setShow = intercalate "," . map show . S.toList---- -- | Map a function over all elements in the 'Bipartition'.--- bpMap :: Ord b => (a -> b) -> Bipartition a -> Bipartition b--- bpMap f (Bipartition (x, y)) = bp (S.map f x) (S.map f y)---- | For a bifurcating root, get the bipartition induced by the root node.------ Return 'Left' if--- - the root node is not bifurcating;--- - a leave set is empty.-bipartition :: Ord a => Tree e a -> Either String (Bipartition a)-bipartition (Node _ _ [x, y]) = bp (S.fromList $ leaves x) (S.fromList $ leaves y)-bipartition _ = Left "bipartition: Root node is not bifurcating."---- | Get all bipartitions of the tree.------ Return 'Left' if the tree contains duplicate leaves.-bipartitions :: Ord a => Tree e a -> Either String (Set (Bipartition a))-bipartitions t-  | duplicateLeaves t = Left "bipartitions: Tree contains duplicate leaves."-  | otherwise = Right $ bipartitions' S.empty $ S.fromList <$> groups t---- | Report the complementary leaves for each child.-getComplementaryLeaves ::-  (Ord a) =>-  -- Complementary leaves.-  Set a ->-  -- Tree with node labels storing leaves.-  Tree e (Set a) ->-  [Set a]-getComplementaryLeaves p (Node _ _ ts) =-  [ S.unions $ p : take i lvsChildren ++ drop (i + 1) lvsChildren-    | i <- [0 .. (n -1)]-  ]-  where-    n = length ts-    lvsChildren = map label ts---- See 'bipartitions', but do not check if leaves are unique, nor if--- bipartitions are valid.-bipartitions' :: Ord a => Set a -> Tree e (Set a) -> Set (Bipartition a)-bipartitions' p (Node _ p' []) = either (const S.empty) S.singleton $ bp p p'-bipartitions' p t@(Node _ p' ts) =-  S.unions $-    either (const S.empty) S.singleton (bp p p') :-      [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'.------ Since the induced bipartitions of the daughter branches of a bifurcating root--- node are equal, the branches leading to the root have to be combined in this--- case. See http://evolution.genetics.washington.edu/phylip/doc/treedist.html--- and how unrooted trees should be handled.------ Further, branches connected to degree two nodes also induce the same--- bipartitions and have to be combined.------ For combining branches, a binary function is required. This requirement is--- encoded in the 'Semigroup' type class constraint (see 'prune').------ Return 'Left' if the tree contains duplicate leaves.-bipartitionToBranch ::-  (Semigroup e, Ord a) =>-  Tree e a ->-  Either String (Map (Bipartition a) e)-bipartitionToBranch t-  | duplicateLeaves t = Left "bipartitionToBranch: Tree contains duplicate leaves."-  | otherwise = Right $ bipartitionToBranch' S.empty pTree-  where-    pTree = S.fromList <$> groups t---- When calculating the map, branches separated by various degree two nodes have--- to be combined. Hence, not only the complementary leaves, but also the branch--- label itself have to be passed along.-bipartitionToBranch' ::-  (Semigroup e, Ord a) =>-  -- Complementary leaves.-  Set a ->-  -- Partition tree.-  Tree e (Set a) ->-  Map (Bipartition a) e-bipartitionToBranch' p t@(Node b p' ts) =-  M.unionsWith (<>) $-    either (const M.empty) (`M.singleton` b) (bp p p') :-      [bipartitionToBranch' c s | (c, s) <- zip cs ts]-  where-    cs = getComplementaryLeaves p t

@@ -1,139 +0,0 @@--- |--- Module      :  ELynx.Data.Tree.Distance--- Description :  Compute distances between trees--- Copyright   :  (c) Dominik Schrempf 2020--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Thu Jun 13 17:15:54 2019.------ Various distance functions for phylogenetic trees (and trees with branch--- lengths in general).------ 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.-module ELynx.Data.Tree.Distance-  ( symmetric,-    incompatibleSplits,-    branchScore,-  )-where---- adjacent,--import Data.Bifunctor-import Data.List-import qualified Data.Map as M-import Data.Monoid-import Data.Set (Set)-import qualified Data.Set as S-import ELynx.Data.Tree.Bipartition-import ELynx.Data.Tree.Measurable-import ELynx.Data.Tree.Partition-import ELynx.Data.Tree.Rooted---- Symmetric difference between two 'Set's.-symmetricDifference :: Ord a => Set a -> Set a -> Set a-symmetricDifference xs ys = S.difference xs ys `S.union` S.difference ys xs---- | Symmetric (Robinson-Foulds) distance between two trees.------ Although a rooted tree data type is used, the distance between the unrooted--- trees is returned.------ Return 'Nothing' if the trees contain different leaves.------ XXX: Comparing a list of trees may recompute bipartitions.-symmetric :: Ord a => Tree e1 a -> Tree e2 a -> Either String Int-symmetric t1 t2-  | S.fromList (leaves t1) /= S.fromList (leaves t2) = Left "symmetric: Trees contain different leaves."-  | otherwise = do-    bps1 <- bipartitions t1-    bps2 <- bipartitions t2-    return $ length $ symmetricDifference bps1 bps2--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---- | Number of incompatible splits.------ Similar to 'symmetric' but all bipartitions induced by multifurcations are--- considered. For a detailed description of how the distance is calculated, see--- 'ELynx.Data.Tree.Bipartition.bipartitionCompatible'.------ A multifurcation on a tree may (but not necessarily does) represent missing--- information about the order of bifurcations. In this case, it is interesting--- to get a set of compatible bifurcations of the tree. For example, the star tree------ > (A,B,C,D);------ induces the following bipartitions:------ > A|BCD--- > B|ACD--- > C|ABD--- > D|ABC------ However, the tree is additionally compatible with the following hidden--- bipartitions:------ > AB|CD--- > AC|BD--- > AD|BC------ For an explanation of how compatibility of partitions is checked, see--- 'compatible'. Before using 'compatible', bipartitions are simply converted to--- partitions with two subsets.------ A bipartition is incompatible with a tree if it is incompatible with all--- induced multifurcations of the tree.------ XXX: Comparing a list of trees with this function recomputes bipartitions.-incompatibleSplits :: (Show a, Ord a) => Tree e1 a -> Tree e2 a -> Either String Int-incompatibleSplits t1 t2-  | S.fromList (leaves t1) /= S.fromList (leaves t2) =-    Left "incompatibleSplits: Trees do not have equal leaf sets."-  | otherwise = do-    -- Bipartitions.-    bs1 <- bipartitions t1-    bs2 <- bipartitions t2-    -- traceShowM $ "bs1" ++ show (S.map bpHuman bs1)-    -- traceShowM $ "bs2" ++ show (S.map bpHuman bs2)-    let -- Putative incompatible bipartitions of trees one and two, respectively.-        putIncBs1 = bs1 S.\\ bs2-        putIncBs2 = bs2 S.\\ bs1-    -- Partitions.-    ms1 <- partitions t1-    ms2 <- partitions t2-    -- traceShowM $ "putIncBs1 " ++ show (S.map bpHuman putIncBs1)-    -- traceShowM $ "putIncBs2 " ++ show (S.map bpHuman putIncBs2)-    return $ countIncompatibilities putIncBs1 ms2 + countIncompatibilities putIncBs2 ms1---- | Compute branch score distance between two trees.------ Although a rooted tree data type is used, the distance between the unrooted--- trees is returned.------ XXX: Comparing a list of trees with this function recomputes bipartitions.-branchScore :: (Measurable e1, Measurable e2, Ord a) => Tree e1 a -> Tree e2 a -> Either String Double-branchScore t1 t2-  | S.fromList (leaves t1) /= S.fromList (leaves t2) = Left "branchScoreWith: Trees do not have equal leaf sets."-  | otherwise = do-    bpToBr1 <- bipartitionToBranch $ first (Sum . getLen) t1-    bpToBr2 <- bipartitionToBranch $ first (Sum . getLen) t2-    let dBs = M.unionWith (-) bpToBr1 bpToBr2-        dsSquared = foldl' (\acc e -> acc + e * e) 0 dBs-    return $ sqrt $ getSum dsSquared

@@ -1,124 +0,0 @@--- |--- Module      :  ELynx.Data.Tree.Measurable--- Description :  Measurable branch labels--- Copyright   :  (c) Dominik Schrempf 2020--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Thu Jan 17 14:16:34 2019.------ Non-negativity of branch lengths is not (yet) ensured. To ensure--- non-negativity, a newtype wrapper could be used, but this would be a major--- refactor.-module ELynx.Data.Tree.Measurable-  ( BranchLength,-    Measurable (..),-    applyStem,-    getStem,-    setStem,-    height,-    rootHeight,-    distancesOriginLeaves,-    totalBranchLength,-    normalizeBranchLengths,-    normalizeHeight,-    ultrametric,-    makeUltrametric,-  )-where--import Data.Bifoldable-import Data.Bifunctor-import ELynx.Data.Tree.Rooted---- | Branch length.-type BranchLength = Double---- | A branch label with measurable and modifiable branch length.-class Measurable e where-  -- | Length of attached branch.-  getLen :: e -> BranchLength--  -- | Set attached branch length.-  setLen :: BranchLength -> e -> e--instance Measurable Double where-  getLen = id-  setLen = const---- Apply a function to a branch support label.-apply :: Measurable e => (BranchLength -> BranchLength) -> e -> e-apply f l = setLen (f s) l where s = getLen l---- | Lengthen the stem of a tree.-applyStem :: Measurable e => (BranchLength -> BranchLength) -> Tree e a -> Tree e a-applyStem f t = t {branch = apply f b}-  where-    b = branch t---- | Get the length of the stem of a tree.-getStem :: Measurable e => Tree e a -> BranchLength-getStem (Node br _ _) = getLen br---- | Set the length of the stem of a tree.-setStem :: Measurable e => BranchLength -> Tree e a -> Tree e a-setStem x = applyStem (const x)---- | The maximum distance between origin and leaves.------ The height includes the length of the stem.-height :: Measurable e => Tree e a -> BranchLength-height = maximum . distancesOriginLeaves---- | The maximum distance between root node and leaves.-rootHeight :: Measurable e => Tree e a -> BranchLength-rootHeight (Node _ _ []) = 0-rootHeight t = maximum $ concatMap distancesOriginLeaves (forest t)---- | Distances from the origin of a tree to the leaves.------ The distances include the length of the stem.-distancesOriginLeaves :: Measurable e => Tree e a -> [BranchLength]-distancesOriginLeaves (Node br _ []) = [getLen br]-distancesOriginLeaves (Node br _ ts) = map (getLen br +) (concatMap distancesOriginLeaves ts)---- | Total branch length of a tree.-totalBranchLength :: Measurable e => Tree e a -> BranchLength-totalBranchLength = bifoldl' (+) const 0 . first getLen---- | Normalize branch lengths so that the sum is 1.0.-normalizeBranchLengths :: Measurable e => Tree e a -> Tree e a-normalizeBranchLengths t = first (apply (/ s)) t-  where-    s = totalBranchLength t---- | Normalize height of tree to 1.0.-normalizeHeight :: Measurable e => Tree e a -> Tree e a-normalizeHeight t = first (apply (/ h)) t-  where-    h = height t--eps :: Double-eps = 1e-12--allNearlyEqual :: [Double] -> Bool-allNearlyEqual [] = True-allNearlyEqual xs = all (\y -> eps > abs (x - y)) (tail xs)-  where-    x = head xs---- | Check if a tree is ultrametric.-ultrametric :: Measurable e => Tree e a -> Bool-ultrametric = allNearlyEqual . distancesOriginLeaves---- | Elongate terminal branches such that the tree becomes ultrametric.-makeUltrametric :: Measurable e => Tree e a -> Tree e a-makeUltrametric t = go 0 t-  where-    h = height t-    go :: Measurable e => BranchLength -> Tree e a -> Tree e a-    go h' (Node br lb []) = let dh = h - h' - getLen br in Node (apply (+ dh) br) lb []-    go h' (Node br lb ts) = let h'' = h' + getLen br in Node br lb $ map (go h'') ts

@@ -1,43 +0,0 @@--- |--- Module      :  ELynx.Data.Tree.Named--- Description :  Trees with named nodes--- Copyright   :  (c) Dominik Schrempf 2020--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Thu Jan 24 20:09:20 2019.-module ELynx.Data.Tree.Named-  ( Named (..),-  )-where--import qualified Data.ByteString.Builder as BB-import qualified Data.ByteString.Char8 as BS-import qualified Data.ByteString.Lazy.Char8 as BL-import Data.Double.Conversion.ByteString as BC---- | Data types with names.-class Named a where-  -- Use lazy byte strings because Newick strings are built using chunks.-  getName :: a -> BL.ByteString--instance Named () where-  getName = const BL.empty--instance Named Int where-  getName = BB.toLazyByteString . BB.intDec--instance Named Double where-  getName = BL.fromStrict . toShortest--instance Named Char where-  getName = BB.toLazyByteString . BB.char8--instance Named BL.ByteString where-  getName = id--instance Named BS.ByteString where-  getName = BL.fromStrict

@@ -1,137 +0,0 @@--- |--- Module      :  ELynx.Data.Tree.Partition--- Description :  Partitions on rose trees--- Copyright   :  (c) Dominik Schrempf 2020--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Thu Dec 12 12:58:49 2019.------ A multifurcation induces a 'Partition', similar to branches inducing--- 'ELynx.Data.Tree.Bipartition's.-module ELynx.Data.Tree.Partition-  ( -- * Data type-    Partition (fromPartition),-    mp,-    mpUnsafe,-    bpToMp,-    mpHuman,--    -- * Work with 'Partition's-    partition,-    partitions,-    compatible,-  )-where--import Data.List hiding (partition)-import Data.Set (Set)-import qualified Data.Set as S-import ELynx.Data.Tree.Bipartition-import ELynx.Data.Tree.Rooted---- | Each branch of a tree partitions the leaves of the tree into two subsets--- (see 'ELynx.Data.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.Data.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.-newtype Partition a = Partition-  { fromPartition :: Set (Set a)-  }-  deriving (Eq, Ord, Show, Read)---- 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-  [] -> Left "mp: Empty list."-  xs' -> Right $ mpUnsafe xs'---- | Create a partition.-mpUnsafe :: Ord a => [Set a] -> Partition a-mpUnsafe xs = Partition (S.fromList xs)---- | Convert a bipartition to a partition.-bpToMp :: Ord a => Bipartition a -> Partition a-bpToMp = mpUnsafe . tupleToList . fromBipartition-  where-    -- Be careful with tuples, because 'toList' does something very weird. It only-    -- takes the second element of the tuple!-    ---    -- toList :: Foldable t => t a -> [a]-    tupleToList (x, y) = [x, y]---- | 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) =-  "(" ++ intercalate "|" (map setShow (S.toList xs)) ++ ")"---- Show the elements of a set in a human readable format.-setShow :: Show a => Set a -> String-setShow = intercalate "," . map show . S.toList---- | Get partition defined by the root of the tree.------ Return 'Left' if:--- - the tree is a leaf;--- - the tree contains duplicate leaves.-partition :: Ord a => Tree e a -> Either String (Partition a)-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---- | Get all 'Partition's of a tree.------ Return 'Left' if tree contains duplicate leaves.-partitions :: Ord a => Tree e a -> Either String (Set (Partition a))-partitions t-  | duplicateLeaves t = Left "partitions: Tree contains duplicate leaves."-  | otherwise = Right $ partitions' S.empty $ S.fromList <$> groups t---- See 'partitions', but do not check if leaves are unique.-partitions' :: Ord a => Set a -> Tree e (Set a) -> Set (Partition a)-partitions' _ (Node _ _ []) = S.empty-partitions' p t@(Node _ _ ts) =-  S.unions $-    either (const S.empty) S.singleton (mp (p : map label ts)) :-    zipWith partitions' cs ts-  where-    cs = getComplementaryLeaves p t---- | 'Partition's are compatible if they do not contain conflicting--- information. This function checks if two partitions are compatible with--- each other. Thereby, a variation of the following algorithm is used:------ @--- mp1 `compatible` mp2--- for set1 in mp1:---   for set2 in mp2:---     if set1 `S.isSubSetOf` set2:---       remove set1 from mp1---     if set2 `S.isSubSetOf` set1:---       remove set2 from mp2--- if either mp2 or mp2 is empty, they are compatible--- @-compatible :: (Show a, Ord a) => Partition a -> Partition a -> Bool-compatible l r = S.null (S.filter (`remove` rs) ls) || S.null (S.filter (`remove` ls) rs)-  where-    ls = fromPartition l-    rs = fromPartition r--remove :: Ord a => Set a -> Set (Set a) -> Bool-remove s = not . any (s `S.isSubsetOf`)

@@ -1,487 +0,0 @@-{-# LANGUAGE DeriveAnyClass #-}-{-# LANGUAGE DeriveGeneric #-}-{-# LANGUAGE DerivingVia #-}---- |--- Module      :  ELynx.Data.Tree.Phylogeny--- Description :  Phylogenetic trees--- Copyright   :  (c) Dominik Schrempf 2020--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Thu Jan 17 16:08:54 2019.------ A phylogeny is a 'Tree' with unique leaf labels, and the order of the trees--- in the sub-forest is considered to be meaningless.------ Internally, however, the underlying 'Tree' data structure stores the--- sub-forest as a list, which has a specific order. Hence, we have to do some--- tricks when comparing trees, and tree comparison is slow.------ Also, the uniqueness of the leaves is not ensured by the data type, but has--- to be checked at runtime. Functions relying on the tree to have unique leaves--- do perform this check, and return 'Left' with an error message, if the tree--- has duplicate leaves.------ Note: 'Tree's are rooted.------ Note: 'Tree's encoded in Newick format correspond to rooted trees. By--- convention only, a tree parsed from Newick format is usually thought to be--- unrooted, when the root node is multifurcating and has three children. This--- convention is not used here. Newick trees are just parsed as they are, and a--- rooted tree is returned.------ The bifurcating root of a tree can be changed with 'rootAt' or 'midpoint'; a--- list of all rooted trees is returned by 'roots'.------ Trees with multifurcating root nodes can be rooted using 'outgroup'.-module ELynx.Data.Tree.Phylogeny-  ( -- * Functions-    equal,-    intersect,-    bifurcating,-    outgroup,-    midpoint,-    roots,-    rootAt,--    -- * Branch labels-    Phylo (..),-    toPhyloTree,-    measurableToPhyloTree,-    supportedToPhyloTree,-    Length (..),-    phyloToLengthTree,-    Support (..),-    phyloToSupportTree,-    phyloToSupportTreeUnsafe,-    PhyloExplicit (..),-    toExplicitTree,-  )-where--import Control.DeepSeq-import Data.Aeson-import Data.Bifoldable-import Data.Bifunctor-import Data.Bitraversable-import Data.List hiding (intersect)-import Data.Maybe-import Data.Monoid-import Data.Semigroup-import Data.Set (Set)-import qualified Data.Set as S-import ELynx.Data.Tree.Bipartition-import ELynx.Data.Tree.Measurable-import ELynx.Data.Tree.Rooted-import ELynx.Data.Tree.Splittable-import ELynx.Data.Tree.Supported-import GHC.Generics---- | The equality check is slow because the order of children is considered to--- be arbitrary.-equal :: (Eq e, Eq a) => Tree e a -> Tree e a -> Bool-equal ~(Node brL lbL tsL) ~(Node brR lbR tsR) =-  (brL == brR)-    && (lbL == lbR)-    && (length tsL == length tsR)-    && all (`elem` tsR) tsL---- | Compute the intersection of trees.------ The intersections are the largest subtrees sharing the same leaf set.------ Degree two nodes are pruned with 'prune'.------ Return 'Left' if:--- - the intersection of leaves is empty.-intersect ::-  (Semigroup e, Eq e, Ord a) => Forest e a -> Either String (Forest e a)-intersect ts-  | S.null lvsCommon = Left "intersect: Intersection of leaves is empty."-  | otherwise = case sequence [dropLeavesWith (predicate ls) t | (ls, t) <- zip leavesToDrop ts] of-    Nothing -> Left "intersect: A tree is empty."-    Just ts' -> Right ts'-  where-    -- Leaf sets.-    lvss = map (S.fromList . leaves) ts-    -- Common leaf set.-    lvsCommon = foldl1' S.intersection lvss-    -- Leaves to drop for each tree in the forest.-    leavesToDrop = map (S.\\ lvsCommon) lvss-    -- Predicate.-    predicate lvsToDr l = l `S.member` lvsToDr---- | Check if a tree is bifurcating.------ A Bifurcating tree only contains degree one (leaves) and degree three nodes--- (internal bifurcating nodes).-bifurcating :: Tree e a -> Bool-bifurcating (Node _ _ []) = True-bifurcating (Node _ _ [x, y]) = bifurcating x && bifurcating y-bifurcating _ = False---- I believe that manual treatment with 'outgroup' is preferable.---- -- | Remove multifurcations.--- ----- -- A caterpillar like bifurcating structure is used to resolve all--- -- multifurcations on a tree.--- ----- -- Multifurcating nodes are copied and branches are 'split'.--- resolve :: Splittable e => Tree e a -> Tree e a--- resolve t@(Node _ _ []) = t--- resolve (Node br lb [x]) = Node br lb [resolve x]--- resolve (Node br lb [x, y]) = Node br lb $ map resolve [x, y]--- resolve (Node br lb (Node brL lbL xsL : xs)) = Node br lb [Node brL' lbL (map resolve xsL), Node brL' lb (map resolve xs)]---   where brL' = split brL---- | Resolve a multifurcating root using an outgroup.------ A bifurcating root node with the provided label is introduced. The affected--- branch is 'split'.------ Note, the degree of the former root node is decreased by one.------ If the root node is bifurcating, use 'rootAt'.------ Return 'Left' if--- - the tree has duplicate leaves;--- - the root node is not multifurcating;--- - the provided outgroup is not found on the tree or is polyphyletic.-outgroup :: (Semigroup e, Splittable e, Ord a) => Set a -> a -> Tree e a -> Either String (Tree e a)-outgroup _ _ (Node _ _ []) = Left "outgroup: Root node is a leaf."-outgroup _ _ (Node _ _ [_]) = Left "outgroup: Root node has degree two."-outgroup _ _ (Node _ _ [_, _]) = Left "outgroup: Root node is bifurcating."-outgroup o r t@(Node b l ts)-  | duplicateLeaves t = Left "outgroup: Tree has duplicate leaves."-  | otherwise = do-    bip <- bp o (S.fromList lvs S.\\ o)-    rootAt bip t'-  where-    lvs = leaves t-    (Node brO lbO tsO) = head ts-    -- Introduce a bifurcating root node.-    t' = Node b r [Node (split brO) lbO tsO, Node (split brO) l (tail ts)]---- The 'midpoint' algorithm is pretty stupid because it calculates all rooted--- trees and then finds the one minimizing the difference between the heights of--- the left and right sub tree. Actually, one just needs to move left or right,--- with the aim to minimize the height difference between the left and right sub--- tree.---- | Root tree at the midpoint.------ Return 'Left' if--- - the root node is not bifurcating.-midpoint :: (Semigroup e, Splittable e, Measurable e) => Tree e a -> Either String (Tree e a)-midpoint (Node _ _ []) = Left "midpoint: Root node is a leaf."-midpoint (Node _ _ [_]) = Left "midpoint: Root node has degree two."-midpoint t@(Node _ _ [_, _]) = getMidpoint <$> roots t-midpoint _ = Left "midpoint: Root node is multifurcating."--findMinIndex :: Ord a => [a] -> Int-findMinIndex (x : xs) = go (0, x) 1 xs-  where-    go (i, _) _ [] = i-    go (i, z) j (y : ys) = if z < y then go (i, z) (j + 1) ys else go (j, y) (j + 1) ys-findMinIndex [] = error "findMinIndex: Empty list."--getMidpoint :: Measurable e => [Tree e a] -> Tree e a-getMidpoint ts = case t of-  (Node br lb [l, r]) ->-    let hl = height l-        hr = height r-        dh = (hl - hr) / 2-     in Node br lb [applyStem (subtract dh) l, applyStem (+ dh) r]-  -- Explicitly use 'error' here, because roots is supposed to return trees with-  -- bifurcating root nodes.-  _ -> error "getMidpoint: Root node is not bifurcating."-  where-    dhs = map getDeltaHeight ts-    i = findMinIndex dhs-    t = ts !! i---- find index of minimum; take this tree and move root to the midpoint of the branch---- Get delta height of left and right sub tree.-getDeltaHeight :: Measurable e => Tree e a -> Double-getDeltaHeight (Node _ _ [l, r]) = abs $ height l - height r--- Explicitly use 'error' here, because roots is supposed to return trees with--- bifurcating root nodes.-getDeltaHeight _ = error "getDeltaHeight: Root node is not bifurcating."---- | For a rooted tree with a bifurcating root node, get all possible rooted--- trees.------ The root node is moved.------ For a tree with @l=2@ leaves, there is one rooted tree. For a bifurcating--- tree with @l>2@ leaves, there are @(2l-3)@ rooted trees. For a general tree--- with a bifurcating root node, and a total number of @n>2@ nodes, there are--- (n-2) rooted trees.------ Moving a multifurcating root node to another branch would change the degree--- of the root node. Hence, a bifurcating root is required. To resolve a--- multifurcating root, please use 'outgroup'.------ Connect branches according to the provided 'Semigroup' instance.------ Upon insertion of the root, split the affected branch into one out of two--- equal entities according to a given function.------ Return 'Left' if the root node is not 'bifurcating'.-roots :: (Semigroup e, Splittable e) => Tree e a -> Either String (Forest e a)-roots (Node _ _ []) = Left "roots: Root node is a leaf."-roots (Node _ _ [_]) = Left "roots: Root node has degree two."-roots t@(Node b c [tL, tR]) = Right $ t : descend b c tR tL ++ descend b c tL tR-roots _ = Left "roots: Root node is multifurcating."--complementaryForests :: Tree e a -> Forest e a -> [Forest e a]-complementaryForests t ts = [t : take i ts ++ drop (i + 1) ts | i <- [0 .. (n -1)]]-  where-    n = length ts---- From the bifurcating root, descend into one of the two pits.------ descend splitFunction rootBranch rootLabel complementaryTree downwardsTree-descend :: (Semigroup e, Splittable e) => e -> a -> Tree e a -> Tree e a -> Forest e a-descend _ _ _ (Node _ _ []) = []-descend brR lbR tC (Node brD lbD tsD) =-  [ Node brR lbR [Node (split brDd) lbD f, Node (split brDd) lbDd tsDd]-    | (Node brDd lbDd tsDd, f) <- zip tsD cfs-  ]-    ++ concat-      [ descend brR lbR (Node (split brDd) lbD f) (Node (split brDd) lbDd tsDd)-        | (Node brDd lbDd tsDd, f) <- zip tsD cfs-      ]-  where-    brC' = branch tC <> brD-    tC' = tC {branch = brC'}-    cfs = complementaryForests tC' tsD---- | Root a tree at a specific position.------ Root the tree at the branch defined by the given bipartition. The original--- root node is moved to the new position.------ The root node must be bifurcating (see 'roots' and 'outgroup').------ Connect branches according to the provided 'Semigroup' instance.------ Upon insertion of the root, split the affected branch according to the--- provided 'Splittable' instance.------ Return 'Left', if:--- - the root node is not bifurcating;--- - the tree has duplicate leaves;--- - the bipartition does not match the leaves of the tree.-rootAt ::-  (Semigroup e, Splittable e, Eq a, Ord a) =>-  Bipartition a ->-  Tree e a ->-  Either String (Tree e a)-rootAt b t-  -- Tree is checked for being bifurcating in 'roots'.-  ---  -- Do not use 'duplicateLeaves' here, because we also need to compare the leaf-  -- set with the bipartition.-  | length lvLst /= S.size lvSet = Left "rootAt: Tree has duplicate leaves."-  | toSet b /= lvSet = Left "rootAt: Bipartition does not match leaves of tree."-  | otherwise = rootAt' b t-  where-    lvLst = leaves t-    lvSet = S.fromList $ leaves t---- Assume the leaves of the tree are unique.-rootAt' ::-  (Semigroup e, Splittable e, Ord a) =>-  Bipartition a ->-  Tree e a ->-  Either String (Tree e a)-rootAt' b t = do-  ts <- roots t-  case find (\x -> Right b == bipartition x) ts of-    Nothing -> Left "rootAt': Bipartition not found on tree."-    Just t' -> Right t'---- | Branch label for phylogenetic trees.------ Branches may have a length and a support value.-data Phylo = Phylo-  { brLen :: Maybe BranchLength,-    brSup :: Maybe BranchSupport-  }-  deriving (Read, Show, Eq, Ord, Generic, NFData)--instance Semigroup Phylo where-  Phylo mBL mSL <> Phylo mBR mSR =-    Phylo-      (getSum <$> (Sum <$> mBL) <> (Sum <$> mBR))-      (getMin <$> (Min <$> mSL) <> (Min <$> mSR))--instance ToJSON Phylo--instance FromJSON Phylo---- | Set all branch length and support values to 'Just' the value.------ Useful to export a tree with branch lengths in Newick format.-toPhyloTree :: (Measurable e, Supported e) => Tree e a -> Tree Phylo a-toPhyloTree = first toPhyloLabel--toPhyloLabel :: (Measurable e, Supported e) => e -> Phylo-toPhyloLabel x = Phylo (Just $ getLen x) (Just $ getSup x)---- | Set all branch support values to 'Nothing'.------ Useful to export a tree with branch lengths to Newick format.-measurableToPhyloTree :: Measurable e => Tree e a -> Tree Phylo a-measurableToPhyloTree = first measurableToPhyloLabel--measurableToPhyloLabel :: Measurable e => e -> Phylo-measurableToPhyloLabel x = Phylo (Just $ getLen x) Nothing---- | Set all branch lengths to 'Nothing'.------ Useful to export a tree with branch support to Newick format.-supportedToPhyloTree :: Supported e => Tree e a -> Tree Phylo a-supportedToPhyloTree = first supportedToPhyloLabel--supportedToPhyloLabel :: Supported e => e -> Phylo-supportedToPhyloLabel x = Phylo Nothing (Just $ getSup x)---- | Branch length label.------ For conversion, see 'phyloToLengthTree' and 'lengthToPhyloTree'.-newtype Length = Length {fromLength :: BranchLength}-  deriving (Read, Show, Eq, Ord, Generic, NFData)-  deriving (Num, Fractional, Floating) via Double-  deriving (Semigroup, Monoid) via Sum Double--instance Measurable Length where-  getLen = fromLength-  setLen b _ = Length b--instance Splittable Length where-  split = Length . (/ 2.0) . fromLength--instance ToJSON Length--instance FromJSON Length---- | If root branch length is not available, set it to 0.------ Return 'Left' if any other branch length is unavailable.-phyloToLengthTree :: Tree Phylo a -> Either String (Tree Length a)-phyloToLengthTree =-  maybe (Left "phyloToLengthTree: Length unavailable for some branches.") Right-    . bitraverse toLength pure-    . cleanRootLength--cleanRootLength :: Tree Phylo a -> Tree Phylo a-cleanRootLength (Node (Phylo Nothing s) l f) = Node (Phylo (Just 0) s) l f-cleanRootLength t = t--toLength :: Phylo -> Maybe Length-toLength p = Length <$> brLen p---- | Branch support label.------ For conversion, see 'phyloToSupportTree'.-newtype Support = Support {fromSupport :: BranchSupport}-  deriving (Read, Show, Eq, Ord, Generic, NFData)-  deriving (Num, Fractional, Floating) via Double-  deriving (Semigroup) via Min Double--instance Supported Support where-  getSup = fromSupport-  setSup s _ = Support s--instance Splittable Support where-  split = id--instance ToJSON Support--instance FromJSON Support---- | Set branch support values of branches leading to the leaves and of the root--- branch to maximum support.------ Return 'Left' if any other branch has no available support value.-phyloToSupportTree :: Tree Phylo a -> Either String (Tree Support a)-phyloToSupportTree t =-  maybe-    (Left "phyloToSupportTree: Support unavailable for some branches.")-    Right-    $ bitraverse toSupport pure $-      cleanLeafSupport m $-        cleanRootSupport m t-  where-    m = getMaxSupport t---- | Set all unavailable branch support values to maximum support.-phyloToSupportTreeUnsafe :: Tree Phylo a -> Tree Support a-phyloToSupportTreeUnsafe t = cleanSupport m t-  where-    m = getMaxSupport t---- If all branch support values are below 1.0, set the max support to 1.0.-getMaxSupport :: Tree Phylo a -> BranchSupport-getMaxSupport = fromJust . max (Just 1.0) . bimaximum . bimap brSup (const Nothing)--cleanRootSupport :: BranchSupport -> Tree Phylo a -> Tree Phylo a-cleanRootSupport maxSup (Node (Phylo b Nothing) l xs) = Node (Phylo b (Just maxSup)) l xs-cleanRootSupport _ t = t--cleanLeafSupport :: BranchSupport -> Tree Phylo a -> Tree Phylo a-cleanLeafSupport s (Node (Phylo b Nothing) l []) = Node (Phylo b (Just s)) l []-cleanLeafSupport s (Node b l xs) = Node b l $ map (cleanLeafSupport s) xs--toSupport :: Phylo -> Maybe Support-toSupport (Phylo _ Nothing) = Nothing-toSupport (Phylo _ (Just s)) = Just $ Support s--cleanSupport :: BranchSupport -> Tree Phylo a -> Tree Support a-cleanSupport maxSup (Node (Phylo _ s) l xs) = Node (Support $ fromMaybe maxSup s) l $ map (cleanSupport maxSup) xs---- | Explicit branch label for phylogenetic trees.-data PhyloExplicit = PhyloExplicit-  { sBrLen :: BranchLength,-    sBrSup :: BranchSupport-  }-  deriving (Read, Show, Eq, Ord, Generic)--instance Semigroup PhyloExplicit where-  PhyloExplicit bL sL <> PhyloExplicit bR sR = PhyloExplicit (bL + bR) (min sL sR)--instance Measurable PhyloExplicit where-  getLen = sBrLen-  setLen b l = l {sBrLen = b}--instance Splittable PhyloExplicit where-  split l = l {sBrLen = b'}-    where-      b' = sBrLen l / 2.0--instance Supported PhyloExplicit where-  getSup = sBrSup-  setSup s l = l {sBrSup = s}--instance ToJSON PhyloExplicit--instance FromJSON PhyloExplicit---- | Conversion to a 'PhyloExplicit' tree.------ See 'phyloToLengthTree' and 'phyloToSupportTree'.-toExplicitTree :: Tree Phylo a -> Either String (Tree PhyloExplicit a)-toExplicitTree t = do-  lt <- first fromLength <$> phyloToLengthTree t-  st <- first fromSupport <$> phyloToSupportTree t-  case zipTreesWith PhyloExplicit const lt st of-    Nothing -> error "toExplicitTree: This is a bug. Can not zip two trees with the same topology."-    Just zt -> return zt

@@ -1,355 +0,0 @@-{-# LANGUAGE DeriveDataTypeable #-}-{-# LANGUAGE DeriveGeneric #-}---- |--- Module      :  ELynx.Data.Tree.Rooted--- Description :  Rooted trees with labeled branches--- Copyright   :  (c) Dominik Schrempf 2020--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Thu Jan 17 09:57:29 2019.------ Rooted 'Tree's differes from a classical rose 'Data.Tree.Tree' in that it has--- labeled branches.------ For rooted topologies, please see 'ELynx.Data.Topology.Rooted'.------ A 'Tree' is defined as:------ @--- data Tree e a = Node---   { branch :: e,---     label :: a,---     forest :: Forest e a---   }--- @------ where------ @--- type Forest e a = [Tree e a]--- @------ This means, that the word 'Node' is reserved for the constructor of a tree,--- and that a 'Node' has an attached 'branch', a 'label', and a sub-'forest'.--- The value constructor /Node/ and the record function /label/ are not to be--- confused. The elements of the sub-forest are often called /children/.------ With respect to phylogenetic analyses, using the 'Tree' data type has some--- disadvantages:------ 1. All trees are rooted. Unrooted trees can be treated with a rooted data--- structure, as it is used here. However, some functions may be meaningless.------ 2. Changing branch labels, node labels, or the topology of the tree are slow--- operations, especially, when the changes are close to the leaves of the tree.------ In mathematical terms: A 'Tree' is a directed acyclic graph without loops,--- with vertex labels, with edge labels. Let me know if this definition is--- incomplete.-module ELynx.Data.Tree.Rooted-  ( -- * Data type-    Tree (..),-    Forest,-    toTreeBranchLabels,-    toTreeNodeLabels,--    -- * Access leaves, branches and labels-    leaves,-    duplicateLeaves,-    branches,-    setBranches,-    labels,-    setLabels,-    identify,--    -- * Change structure-    degree,-    prune,-    dropNodesWith,-    dropLeavesWith,-    zipTreesWith,-    zipTrees,-  )-where--import Control.Applicative-import Control.Comonad-import Control.DeepSeq-import Control.Monad-import Control.Monad.Fix-import Data.Aeson-import Data.Bifoldable-import Data.Bifunctor-import Data.Bitraversable-import Data.Data-import Data.Foldable-import Data.List-import Data.Maybe-import qualified Data.Set as S-import qualified Data.Tree as T-import GHC.Generics---- | 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.------ Lifted instances are not provided.-data Tree e a = Node-  { branch :: e,-    label :: a,-    forest :: Forest e a-  }-  deriving (Eq, Read, Show, Data, Generic)---- | A shorthand.-type Forest e a = [Tree e a]---- | Map over node labels.-instance Functor (Tree e) where-  fmap f ~(Node br lb ts) = Node br (f lb) $ map (fmap f) ts-  x <$ ~(Node br _ ts) = Node br x (map (x <$) ts)--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-  second g ~(Node br lb ts) = Node br (g lb) $ map (second g) ts---- | Combine node labels in pre-order.-instance Foldable (Tree e) where-  foldMap f ~(Node _ lb ts) = f lb <> foldMap (foldMap f) ts-  null _ = False-  {-# INLINE null #-}-  toList = labels-  {-# INLINE toList #-}--instance Bifoldable Tree where-  bifoldMap f g ~(Node br lb ts) = f br <> g lb <> foldMap (bifoldMap f g) ts--instance Traversable (Tree e) where-  traverse g ~(Node br lb ts) = Node br <$> g lb <*> traverse (traverse g) ts--instance Bitraversable Tree where-  bitraverse f g ~(Node br lb ts) = Node <$> f br <*> g lb <*> traverse (bitraverse f g) ts---- The following code provides a zip-like applicative instance. However,--- the zip-like instance makes the Monad instance meaningless. So, either we--- provide only 'Applicative' in zip-like form, or we use the classic instance--- for 'Applicative' and 'Monad'.---- -- | Note: The 'Applicative' instance of 'Tree' is similar to the one of--- -- 'Control.Applicative.ZipList', and differs from the instance of--- -- 'Data.Tree.Tree'!--- ----- -- >>> let t = Node "" 0 [Node "" 1 [], Node "" 2 []] :: Tree String Int--- -- >>> let f = Node "+3" (+3) [Node "*5" (*5) [], Node "+10" (+10) []] :: Tree String (Int -> Int)--- -- >>> f <*> t--- -- Node {branch = "+3", label = 3, forest = [Node {branch = "*5", label = 5, forest = []},Node {branch = "+10", label = 12, forest = []}]}--- ----- -- Note: The 'Monoid' instance of the branch labels determines how the branches--- -- are combined. For example, distances can be summed using the--- -- 'Data.Monoid.Sum' monoid.--- instance Monoid e => Applicative (Tree e) where---   pure lb = Node mempty lb []---   ~(Node brF lbF tsF) <*> ~(Node brX lbX tsX) =---     Node (brF <> brX) (lbF lbX) (zipWith (<*>) tsF tsX)---   liftA2 f ~(Node brX lbX tsX) ~(Node brY lbY tsY) =---     Node (brX <> brY) (f lbX lbY) (zipWith (liftA2 f) tsX tsY)---   ~(Node brX _ tsX) *> ~(Node brY lbY tsY) =---     Node (brX <> brY) lbY (zipWith (*>) tsX tsY)---   ~(Node brX lbX tsX) <* ~(Node brY _ tsY) =---     Node (brX <> brY) lbX (zipWith (<*) tsX tsY)---- | The 'Semigroup' instance of the branch labels determines how the--- branches are combined. For example, distances can be summed using--- 'Data.Semigroup.Sum'.------ The 'Monoid' instance of the branch labels determines the default branch--- label when using 'pure'.-instance Monoid e => Applicative (Tree e) where-  pure lb = Node mempty lb []-  ~(Node brF lbF tsF) <*> ~tx@(Node brX lbX tsX) =-    Node (brF <> brX) (lbF lbX) (map (lbF <$>) tsX ++ map (<*> tx) tsF)-  liftA2 f ~(Node brX lbX tsX) ~ty@(Node brY lbY tsY) =-    Node (brX <> brY) (f lbX lbY) (map (f lbX <$>) tsY ++ map (\tx -> liftA2 f tx ty) tsX)-  ~(Node brX _ tsX) *> ~ty@(Node brY lbY tsY) =-    Node (brX <> brY) lbY (tsY ++ map (*> ty) tsX)-  ~(Node brX lbX tsX) <* ~ty@(Node brY _ tsY) =-    Node (brX <> brY) lbX (map (lbX <$) tsY ++ map (<* ty) tsX)---- | The 'Semigroup' instance of the branch labels determines how the branches--- are combined. For example, distances can be summed using--- 'Data.Semigroup.Sum'.-instance Monoid e => Monad (Tree e) where-  ~(Node br lb ts) >>= f = case f lb of-    Node br' lb' ts' -> Node (br <> br') lb' (ts' ++ map (>>= f) ts)---- -- Cannot provide MonadZip instance because branch labels cannot be--- -- recovered from combined label.--- instance Monoid e => MonadZip (Tree e) where---   mzipWith f (Node brL lbL tsL) (Node brR lbR tsR) =---     Node (brL <> brR) (f lbL lbR) (mzipWith (mzipWith f) tsL tsR)------   munzip (Node br (lbL, lbR) ts) = (Node ? lbL tsL, Node ? lbR tsR)---     where---       (tsL, tsR) = munzip (map munzip ts)--instance Monoid e => MonadFix (Tree e) where-  mfix = mfixTree--mfixTree :: (a -> Tree e a) -> Tree e a-mfixTree f-  | Node br lb ts <- fix (f . label) =-    Node-      br-      lb-      ( zipWith-          (\i _ -> mfixTree ((!! i) . forest . f))-          [0 ..]-          ts-      )--instance Comonad (Tree e) where-  duplicate t@(Node br _ ts) = Node br t (map duplicate ts)-  extract (Node _ lb _) = lb-  {-# INLINE extract #-}--instance (NFData e, NFData a) => NFData (Tree e a) where-  rnf (Node br lb ts) = rnf br `seq` rnf lb `seq` rnf ts--instance (ToJSON e, ToJSON a) => ToJSON (Tree e a)--instance (FromJSON e, FromJSON a) => FromJSON (Tree e a)---- | Conversion to 'T.Tree' using branch labels.-toTreeBranchLabels :: Tree e a -> T.Tree e-toTreeBranchLabels (Node br _ ts) = T.Node br (map toTreeBranchLabels ts)---- | Conversion to 'T.Tree' using node labels.-toTreeNodeLabels :: Tree e a -> T.Tree a-toTreeNodeLabels (Node _ lb ts) = T.Node lb (map toTreeNodeLabels ts)---- | Get leaves.-leaves :: Tree e a -> [a]-leaves (Node _ lb []) = [lb]-leaves (Node _ _ ts) = concatMap leaves ts---- | Check if a tree has duplicate leaves.-duplicateLeaves :: Ord a => Tree e a -> Bool-duplicateLeaves = duplicates . leaves---- | Get branch labels in pre-order.-branches :: Tree e a -> [e]-branches t = squish t []-  where-    squish (Node br _ ts) xs = br : foldr squish xs ts---- | Set branch labels in pre-order.------ Return 'Nothing' if the provided list of branch labels is too short.-setBranches :: Bitraversable t => [f] -> t e a -> Maybe (t f a)-setBranches xs = bisequenceA . snd . bimapAccumL setBranch noChange xs-  where-    setBranch [] _ = ([], Nothing)-    setBranch (y : ys) _ = (ys, Just y)-    noChange ys z = (ys, Just z)---- | Return node labels in pre-order.-labels :: Tree e a -> [a]-labels t = squish t []-  where-    squish (Node _ lb ts) xs = lb : foldr squish xs ts---- | Set node labels in pre-order.------ Return 'Nothing' if the provided list of node labels is too short.-setLabels :: Traversable t => [b] -> t a -> Maybe (t b)-setLabels xs = sequenceA . snd . mapAccumL setLabel xs-  where-    setLabel [] _ = ([], Nothing)-    setLabel (y : ys) _ = (ys, Just y)---- | Label the nodes with unique integers starting at the root with 0.-identify :: Traversable t => t a -> t Int-identify = snd . mapAccumL (\i _ -> (i + 1, i)) (0 :: Int)---- | The degree of the root node.-degree :: Tree e a -> Int-degree = (+ 1) . length . forest---- | Prune degree two nodes.------ The information stored in a pruned node is lost. The branches are combined--- according to their 'Semigroup' instance of the form @\daughterBranch--- parentBranch -> combinedBranch@.-prune :: Semigroup e => Tree e a -> Tree e a-prune t@(Node _ _ []) = t-prune (Node paBr _ [Node daBr daLb daTs]) = Node (daBr <> paBr) daLb daTs-prune (Node paBr paLb paTs) = Node paBr paLb $ map prune paTs---- | Drop nodes satisfying predicate.------ Degree two nodes may arise.------ Also drop parent nodes of which all daughter nodes are dropped.------ Return 'Nothing' if the root node satisfies the predicate.-dropNodesWith :: (a -> Bool) -> Tree e a -> Maybe (Tree e a)-dropNodesWith p (Node br lb ts)-  | p lb = Nothing-  | otherwise =-    if null ts'-      then Nothing-      else Just $ Node br lb ts'-  where-    ts' = mapMaybe (dropNodesWith p) ts---- | Drop leaves satisfying predicate.------ Degree two nodes may arise.------ Also drop parent nodes of which all leaves are dropped.------ Return 'Nothing' if all leaves satisfy the predicate.-dropLeavesWith :: (a -> Bool) -> Tree e a -> Maybe (Tree e a)-dropLeavesWith p (Node br lb [])-  | p lb = Nothing-  | otherwise = Just $ Node br lb []-dropLeavesWith p (Node br lb ts) =-  if null ts'-    then Nothing-    else Just $ Node br lb ts'-  where-    ts' = mapMaybe (dropLeavesWith p) ts---- | Zip two trees with the same topology.------ Return 'Nothing' if the topologies are different.-zipTreesWith ::-  (e1 -> e2 -> e) ->-  (a1 -> a2 -> a) ->-  Tree e1 a1 ->-  Tree e2 a2 ->-  Maybe (Tree e a)-zipTreesWith f g (Node brL lbL tsL) (Node brR lbR tsR) =-  if length tsL == length tsR-    then -- I am proud of that :)).-      zipWithM (zipTreesWith f g) tsL tsR >>= Just . Node (f brL brR) (g lbL lbR)-    else Nothing---- | Zip two trees with the same topology.------ Return 'Nothing' if the topologies are different.-zipTrees :: Tree e1 a1 -> Tree e2 a2 -> Maybe (Tree (e1, e2) (a1, a2))-zipTrees = zipTreesWith (,) (,)--duplicates :: Ord a => [a] -> Bool-duplicates = go S.empty-  where-    go _ [] = False-    go seen (x : xs) = x `S.member` seen || go (S.insert x seen) xs

@@ -1,29 +0,0 @@--- |--- Module      :  ELynx.Data.Tree.Splittable--- Description :  Splittable branch labels--- Copyright   :  (c) Dominik Schrempf, 2020--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Sat Jul 18 13:52:22 2020.-module ELynx.Data.Tree.Splittable-  ( Splittable (..),-  )-where---- | A data type that can be combined using '<>' and split into one out of two--- equal entities.------ The following equality should hold:------ @--- split x <> split x = x--- @-class Splittable e where-  split :: e -> e--instance Splittable Double where-  split = (/ 2)

@@ -1,66 +0,0 @@--- |--- Module      :  ELynx.Data.Tree.Supported--- Description :  Branch label with support value--- Copyright   :  (c) Dominik Schrempf 2020--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Thu Jun 13 14:06:45 2019.------ Non-negativity of branch support values is not (yet) ensured. To ensure--- non-negativity, a newtype wrapper could be used, but this would be a major--- refactor.-module ELynx.Data.Tree.Supported-  ( BranchSupport,-    Supported (..),-    normalizeBranchSupport,-    collapse,-  )-where--import Data.Bifoldable-import Data.Bifunctor-import Data.List-import ELynx.Data.Tree.Rooted---- | Branch support.-type BranchSupport = Double---- | A branch label that supports extraction and setting of branch support values.-class Supported e where-  getSup :: e -> BranchSupport-  setSup :: BranchSupport -> e -> e---- Apply a function to a branch support label.-apply :: Supported e => (BranchSupport -> BranchSupport) -> e -> e-apply f l = setSup (f s) l where s = getSup l---- | Normalize branch support values. The maximum branch support value will be--- set to 1.0.-normalizeBranchSupport :: Supported e => Tree e a -> Tree e a-normalizeBranchSupport t = first (apply (/ m)) t-  where-    m = bimaximum $ bimap getSup (const 0) t---- | Collapse branches with support lower than given value.------ The branch and node labels of the collapsed branches are discarded.-collapse :: (Eq e, Eq a, Supported e) => BranchSupport -> Tree e a -> Tree e a-collapse th tr =-  let tr' = collapse' th tr-   in if tr == tr' then tr else collapse th tr'---- A leaf has full support.-highP :: Supported e => Double -> Tree e a -> Bool-highP _  (Node _ _ []) = True-highP th (Node br _ _) = getSup br >= th---- See 'collapse'.-collapse' :: Supported e => BranchSupport -> Tree e a -> Tree e a-collapse' th (Node br lb ts) = Node br lb $ map (collapse' th) (highSupport ++ lowSupportForest)-  where-    (highSupport, lowSupport) = partition (highP th) ts-    lowSupportForest = concatMap forest lowSupport

@@ -1,160 +0,0 @@--- |--- Module      :  ELynx.Data.Tree.Zipper--- Description :  Zippers on rooted rose trees with branch labels--- Copyright   :  (c) Dominik Schrempf, 2020--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Thu Jul 23 08:42:37 2020.-module ELynx.Data.Tree.Zipper-  ( -- * Data type-    TreePos (..),--    -- * Conversion-    fromTree,-    toTree,--    -- * Movement-    goUp,-    goRoot,-    goLeft,-    goRight,-    goChild,-    goPath,-    unsafeGoPath,--    -- * Modification-    insertTree,-    insertBranch,-    insertLabel,-  )-where--import Data.Foldable-import ELynx.Data.Tree.Rooted---- | Tree zipper. For reference, please see http://hackage.haskell.org/package/rosezipper.-data TreePos e a = Pos-  { -- | The currently selected tree.-    current :: Tree e a,-    -- | Forest to the left in reversed order.-    before :: Forest e a,-    -- | Forest to the right-    after :: Forest e a,-    -- | Finger to the selected tree-    parents :: [([Tree e a], e, a, [Tree e a])]-  }-  deriving (Show, Eq)---- | Get a zipper pointing to the root.-fromTree :: Tree e a -> TreePos e a-fromTree t = Pos {current = t, before = [], after = [], parents = []}---- | Get the complete tree of the zipper.-toTree :: TreePos e a -> Tree e a-toTree = current . goRoot--getForest :: TreePos e a -> Forest e a-getForest pos = foldl (flip (:)) (current pos : after pos) (before pos)---- | Go to parent.-goUp :: TreePos e a -> Maybe (TreePos e a)-goUp pos = case parents pos of-  (ls, br, lb, rs) : ps ->-    Just-      Pos-        { current = Node br lb $ getForest pos,-          before = ls,-          after = rs,-          parents = ps-        }-  [] -> Nothing---- | Go to root.-goRoot :: TreePos e a -> TreePos e a-goRoot pos = maybe pos goRoot (goUp pos)---- | Go to left sibling in current forest.-goLeft :: TreePos e a -> Maybe (TreePos e a)-goLeft pos =-  case before pos of-    t : ts ->-      Just-        pos-          { current = t,-            before = ts,-            after = current pos : after pos-          }-    [] -> Nothing---- | Go to right sibling in current forest.-goRight :: TreePos e a -> Maybe (TreePos e a)-goRight pos =-  case after pos of-    t : ts ->-      Just-        pos-          { current = t,-            before = current pos : before pos,-            after = ts-          }-    [] -> Nothing---- | Go to child with given index in forest.-goChild :: Int -> TreePos e a -> Maybe (TreePos e a)-goChild n pos = case current pos of-  (Node br lb ts)-    | null ts -> Nothing-    | length ts <= n -> Nothing-    | otherwise ->-      Just $-        Pos-          { current = head rs',-            before = reverse ls',-            after = tail rs',-            parents = (before pos, br, lb, after pos) : parents pos-          }-    where-      (ls', rs') = splitAt n ts---- | Go to node with given path.-goPath :: [Int] -> TreePos e a -> Maybe (TreePos e a)-goPath pos pth = foldlM (flip goChild) pth pos---- | Go to child with given index in forest. Call 'error' if child does not--- exist.-unsafeGoChild :: Int -> TreePos e a -> TreePos e a-unsafeGoChild n pos = case current pos of-  (Node br lb ts)-    | null ts -> error "unsafeGoChild: Forest is empty."-    | length ts <= n -> error "unsafeGoChild: Forest is too short."-    | otherwise ->-      Pos-        { current = head rs',-          before = reverse ls',-          after = tail rs',-          parents = (before pos, br, lb, after pos) : parents pos-        }-    where-      (ls', rs') = splitAt n ts---- | Got to node with given path. Call 'error' if path is invalid.-unsafeGoPath :: [Int] -> TreePos e a -> TreePos e a-unsafeGoPath pos pth = foldl (flip unsafeGoChild) pth pos---- | Insert a new tree into the current focus of the zipper.-insertTree :: Tree e a -> TreePos e a -> TreePos e a-insertTree t pos = pos {current = t}---- | Insert a new branch label into the current focus of the zipper.-insertBranch :: e -> TreePos e a -> TreePos e a-insertBranch br pos = case current pos of-  Node _ lb ts -> pos {current = Node br lb ts}---- | Insert a new node label into the current focus of the zipper.-insertLabel :: a -> TreePos e a -> TreePos e a-insertLabel lb pos = case current pos of-  Node br _ ts -> pos {current = Node br lb ts}

@@ -1,95 +0,0 @@-{-# LANGUAGE DeriveDataTypeable #-}-{-# LANGUAGE DeriveGeneric #-}---- |---   Module      :  ELynx.Distribution.BirthDeath---   Description :  Birth and death distribution---   Copyright   :  (c) Dominik Schrempf 2018---   License     :  GPL-3.0-or-later------   Maintainer  :  dominik.schrempf@gmail.com---   Stability   :  unstable---   Portability :  portable------ Creation date: Tue Feb 13 13:16:18 2018.------ See Gernhard, T. (2008). The conditioned reconstructed process. Journal of--- Theoretical Biology, 253(4), 769–778. http://doi.org/10.1016/j.jtbi.2008.04.005.------ Distribution of the values of the point process such that it corresponds to--- reconstructed trees under the birth and death process.-module ELynx.Distribution.BirthDeath-  ( BirthDeathDistribution (..),-    cumulative,-    density,-    quantile,-  )-where--import Data.Data-  ( Data,-    Typeable,-  )-import ELynx.Distribution.Types-import GHC.Generics (Generic)-import qualified Statistics.Distribution as D---- | Distribution of the values of the point process such that it corresponds to--- a reconstructed tree of the birth and death process.-data BirthDeathDistribution = BDD-  { -- | Time to origin of the tree.-    bddTOr :: Time,-    -- | Birth rate.-    bddLa :: Rate,-    -- | Death rate.-    bddMu :: Rate-  }-  deriving (Eq, Typeable, Data, Generic)--instance D.Distribution BirthDeathDistribution where-  cumulative = cumulative---- | Cumulative distribution function Eq. (3).-cumulative :: BirthDeathDistribution -> Time -> Double-cumulative (BDD t l m) x-  | x <= 0 = 0-  | x > t = 1-  | otherwise = t1 * t2-  where-    d = l - m-    t1 = (1.0 - exp (- d * x)) / (l - m * exp (- d * x))-    t2 = (l - m * exp (- d * t)) / (1.0 - exp (- d * t))--instance D.ContDistr BirthDeathDistribution where-  density = density-  quantile = quantile---- | Density function Eq. (2).-density :: BirthDeathDistribution -> Time -> Double-density (BDD t l m) x-  | x < 0 = 0-  | x > t = 0-  | otherwise = d ** 2 * t1 * t2-  where-    d = l - m-    t1 = exp (- d * x) / ((l - m * exp (- d * x)) ** 2)-    t2 = (l - m * exp (- d * t)) / (1.0 - exp (- d * t))---- | Inverted cumulative probability distribution 'cumulative'. See also--- 'D.ContDistr'.-quantile :: BirthDeathDistribution -> Double -> Time-quantile (BDD t l m) p-  | p >= 0 && p <= 1 =-    res-  | otherwise =-    error $-      "PointProcess.quantile: p must be in range [0,1] but got "-        ++ show p-        ++ "."-  where-    d = l - m-    t2 = (l - m * exp (- d * t)) / (1.0 - exp (- d * t))-    res = (-1.0 / d) * log ((1.0 - p * l / t2) / (1.0 - p * m / t2))--instance D.ContGen BirthDeathDistribution where-  genContVar = D.genContinuous

@@ -1,84 +0,0 @@-{-# LANGUAGE DeriveDataTypeable #-}-{-# LANGUAGE DeriveGeneric #-}---- |---   Module      :  ELynx.Distribution.BirthDeathCritical---   Description :  Birth and death distribution---   Copyright   :  (c) Dominik Schrempf 2018---   License     :  GPL-3.0-or-later------   Maintainer  :  dominik.schrempf@gmail.com---   Stability   :  unstable---   Portability :  portable------ Creation date: Tue Feb 13 13:16:18 2018.------ See Gernhard, T. (2008). The conditioned reconstructed process. Journal of--- Theoretical Biology, 253(4), 769–778. http://doi.org/10.1016/j.jtbi.2008.04.005.------ Distribution of the values of the point process such that it corresponds to--- reconstructed trees under the birth and death process; critical birth and death--- process with lambda=mu.-module ELynx.Distribution.BirthDeathCritical-  ( BirthDeathCriticalDistribution (..),-    cumulative,-    density,-    quantile,-  )-where--import Data.Data-  ( Data,-    Typeable,-  )-import ELynx.Distribution.Types-import GHC.Generics (Generic)-import qualified Statistics.Distribution as D---- | Distribution of the values of the point process such that it corresponds to--- a reconstructed tree of the birth and death process.-data BirthDeathCriticalDistribution = BDCD-  { -- | Time to origin of the tree.-    bdcdTOr :: Time,-    -- | Birth and death rate.-    bdcdLa :: Rate-  }-  deriving (Eq, Typeable, Data, Generic)--instance D.Distribution BirthDeathCriticalDistribution where-  cumulative = cumulative---- | Cumulative distribution function section 2.1.2, second formula.-cumulative :: BirthDeathCriticalDistribution -> Time -> Double-cumulative (BDCD t l) x-  | x <= 0 = 0-  | x > t = 1-  | otherwise = x / (1.0 + l * x) * (1.0 + l * t) / t--instance D.ContDistr BirthDeathCriticalDistribution where-  density = density-  quantile = quantile---- | Density function section 2.1.2, first formula.-density :: BirthDeathCriticalDistribution -> Time -> Double-density (BDCD t l) x-  | x < 0 = 0-  | x > t = 0-  | otherwise = (1.0 + l * t) / (t * (1.0 + l * x) ** 2)---- | Inverted cumulative probability distribution 'cumulative'. See also--- 'D.ContDistr'.-quantile :: BirthDeathCriticalDistribution -> Double -> Time-quantile (BDCD t l) p-  | p >= 0 && p <= 1 =-    res-  | otherwise =-    error $-      "PointProcess.quantile: p must be in [0,1] range. Got: "-        ++ show p-        ++ "."-  where-    res = p * t / (1 + l * t - l * p * t)--instance D.ContGen BirthDeathCriticalDistribution where-  genContVar = D.genContinuous

@@ -1,76 +0,0 @@-{-# LANGUAGE DeriveDataTypeable #-}-{-# LANGUAGE DeriveGeneric #-}---- |---   Module      :  ELynx.Distribution.BirthDeathCriticalNoTime---   Description :  Birth and death distribution---   Copyright   :  (c) Dominik Schrempf 2018---   License     :  GPL-3.0-or-later------   Maintainer  :  dominik.schrempf@gmail.com---   Stability   :  unstable---   Portability :  portable------ Creation date: Tue Feb 13 13:16:18 2018.------ See Gernhard, T. (2008). The conditioned reconstructed process. Journal of--- Theoretical Biology, 253(4), 769–778. http://doi.org/10.1016/j.jtbi.2008.04.005.------ Distribution of the values of the point process such that it corresponds to--- reconstructed trees under the birth and death process; critical birth and death--- process with lambda=mu.-module ELynx.Distribution.BirthDeathCriticalNoTime-  ( BirthDeathCriticalNoTimeDistribution (..),-    cumulative,-    density,-    quantile,-  )-where--import Data.Data-  ( Data,-    Typeable,-  )-import ELynx.Distribution.Types-import GHC.Generics (Generic)-import qualified Statistics.Distribution as D---- | Distribution of the values of the point process such that it corresponds to--- a reconstructed tree of the birth and death process.-newtype BirthDeathCriticalNoTimeDistribution = BDCNTD-  { -- | Birth and death rate.-    bdcntdLa :: Rate-  }-  deriving (Eq, Typeable, Data, Generic)--instance D.Distribution BirthDeathCriticalNoTimeDistribution where-  cumulative = cumulative---- | Cumulative distribution function section 2.1.2, second formula.-cumulative :: BirthDeathCriticalNoTimeDistribution -> Time -> Double-cumulative (BDCNTD l) x-  | x <= 0 = 0-  | otherwise = x * l / (1.0 + x * l)--instance D.ContDistr BirthDeathCriticalNoTimeDistribution where-  density = density-  quantile = quantile---- | Density function section 2.1.2, first formula; t cancels out because it is--- expected to be much larger than 1.0; because t \in [0, \infty].-density :: BirthDeathCriticalNoTimeDistribution -> Time -> Double-density (BDCNTD l) x-  | x < 0 = 0-  | otherwise = l / ((1.0 + x * l) ** 2)---- | Inverted cumulative probability distribution 'cumulative'. See also--- 'D.ContDistr'.-quantile :: BirthDeathCriticalNoTimeDistribution -> Double -> Time-quantile (BDCNTD l) p-  | p >= 0 && p <= 1 =-    p / (l - l * p)-  | otherwise =-    error $ "PointProcess.quantile: p must be in [0,1]. Got: " ++ show p ++ "."--instance D.ContGen BirthDeathCriticalNoTimeDistribution where-  genContVar = D.genContinuous

@@ -1,109 +0,0 @@-{-# LANGUAGE DeriveDataTypeable #-}-{-# LANGUAGE DeriveGeneric #-}---- |---   Module      :  ELynx.Distribution.BirthDeathNearlyCritical---   Description :  Birth and death distribution---   Copyright   :  (c) Dominik Schrempf 2018---   License     :  GPL-3.0-or-later------   Maintainer  :  dominik.schrempf@gmail.com---   Stability   :  unstable---   Portability :  portable------ Creation date: Tue Feb 13 13:16:18 2018.------ See Gernhard, T. (2008). The conditioned reconstructed process. Journal of--- Theoretical Biology, 253(4), 769–778. http://doi.org/10.1016/j.jtbi.2008.04.005.------ Distribution of the values of the point process such that it corresponds to--- reconstructed trees under the birth and death process; nearly critical birth and--- death process with lambda~mu.------ Basically, this is a Taylor expansion of Eq. (2) and Eq. (3).-module ELynx.Distribution.BirthDeathNearlyCritical-  ( BirthDeathNearlyCriticalDistribution (..),-    cumulative,-    density,-    quantile,-  )-where--import Data.Data-  ( Data,-    Typeable,-  )-import ELynx.Distribution.Types-import GHC.Generics (Generic)-import qualified Statistics.Distribution as D---- | Distribution of the values of the point process such that it corresponds to--- a reconstructed tree of the birth and death process.-data BirthDeathNearlyCriticalDistribution = BDNCD-  { -- | Time to origin of the tree.-    bdncdTOr :: Time,-    -- | Birth and death rate.-    bdncdLa :: Rate,-    -- | Birth and death rate.-    bdncdMu :: Rate-  }-  deriving (Eq, Typeable, Data, Generic)--instance D.Distribution BirthDeathNearlyCriticalDistribution where-  cumulative = cumulative---- | Cumulative distribution function section 2.1.2, second formula.-cumulative :: BirthDeathNearlyCriticalDistribution -> Time -> Double-cumulative (BDNCD t l m) s-  | s <= 0 = 0-  | s > t = 1-  | otherwise = o0 + o1-  where-    o0 = s * (1.0 + t * l) / t / (1.0 + s * l)-    o1 = (- s * s + s * t) * (m - l) / (2.0 * t * (1.0 + s * l) ** 2)--instance D.ContDistr BirthDeathNearlyCriticalDistribution where-  density = density-  quantile = quantile---- | Density function section 2.1.2, first formula.-density :: BirthDeathNearlyCriticalDistribution -> Time -> Double-density (BDNCD t l m) s-  | s < 0 = 0-  | s > t = 0-  | otherwise = o0 + o1-  where-    o0 = (1.0 + t * l) / (t * (1.0 + s * l) ** 2)-    o1 = (-2.0 * s + t - s * t * l) * (m - l) / (2.0 * t * (1.0 + s * l) ** 3)---- | Inverted cumulative probability distribution 'cumulative'. See also--- 'D.ContDistr'.-quantile :: BirthDeathNearlyCriticalDistribution -> Double -> Time-quantile (BDNCD t l m) p-  | p >= 0 && p <= 1 =-    res-  | otherwise =-    error $-      "PointProcess.quantile: p must be in [0,1] range. Got: "-        ++ show p-        ++ "."-  where-    den = l * (-3.0 + 2.0 * t * (-1.0 + p) * l) + m-    t1 = (2.0 + t * (l - 4.0 * p * l + m)) / den-    t2Nom =-      4.0-        + t-        * ( l-              * (4.0 + t * l + 8.0 * p * (1.0 + t * l))-              + 2.0-              * (2.0 + t * l - 4.0 * p * (1.0 + t * l))-              * m-              + t-              * m-              * m-          )-    t2 = t2Nom / (den ** 2)-    res = 0.5 * (t1 + sqrt t2)--instance D.ContGen BirthDeathNearlyCriticalDistribution where-  genContVar = D.genContinuous

@@ -1,26 +0,0 @@--- |--- Module      :  ELynx.Distribution.CoalescentContinuous--- Description :  Distribution of coalescent times--- Copyright   :  (c) Dominik Schrempf 2018--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Wed May 16 12:40:45 2018.-module ELynx.Distribution.CoalescentContinuous-  ( coalescentDistributionCont,-  )-where--import Numeric.SpecFunctions (choose)-import Statistics.Distribution.Exponential---- | Distribution of the next coalescent event for a number of samples @n@. The--- time is measured in units of effective number of population size.-coalescentDistributionCont ::-  -- | Sample size.-  Int ->-  ExponentialDistribution-coalescentDistributionCont n = exponential (choose n 2)

@@ -1,95 +0,0 @@-{-# LANGUAGE DeriveDataTypeable #-}-{-# LANGUAGE DeriveGeneric #-}---- |---   Module      :  ELynx.Distribution.TimeOfOrigin---   Description :  Distribution of time of origin for birth and death trees---   Copyright   :  (c) Dominik Schrempf 2018---   License     :  GPL-3.0-or-later------   Maintainer  :  dominik.schrempf@gmail.com---   Stability   :  unstable---   Portability :  portable------ Creation date: Tue Feb 13 13:16:18 2018.------ See Gernhard, T. (2008). The conditioned reconstructed process. Journal of--- Theoretical Biology, 253(4), 769–778. http://doi.org/10.1016/j.jtbi.2008.04.005.------ Distribution of the time of origin for birth and death trees. See corollary 3.3--- in the paper cited above.-module ELynx.Distribution.TimeOfOrigin-  ( TimeOfOriginDistribution (..),-    cumulative,-    density,-    quantile,-  )-where--import Data.Data-  ( Data,-    Typeable,-  )-import ELynx.Distribution.Types-import GHC.Generics (Generic)-import qualified Statistics.Distribution as D---- | Distribution of the time of origin for a phylogenetic tree evolving under--- the birth and death process and conditioned on observing n leaves today.-data TimeOfOriginDistribution = TOD-  { -- | Number of leaves of the tree.-    todTN :: Int,-    -- | Birth rate.-    todLa :: Rate,-    -- | Death rate.-    todMu :: Rate-  }-  deriving (Eq, Typeable, Data, Generic)--instance D.Distribution TimeOfOriginDistribution where-  cumulative = cumulative---- | Cumulative distribution function Corollary 3.3.-cumulative :: TimeOfOriginDistribution -> Time -> Double-cumulative (TOD n l m) x-  | x <= 0 = 0-  | otherwise = te ** fromIntegral n-  where-    d = l - m-    te = l * (1.0 - exp (- d * x)) / (l - m * exp (- d * x))--instance D.ContDistr TimeOfOriginDistribution where-  density = density-  quantile = quantile---- | The density function Eq. (5).-density :: TimeOfOriginDistribution -> Time -> Double-density (TOD nn l m) x-  | x < 0 = 0-  | otherwise = n * l ** n * d ** 2 * t1 ** (n - 1.0) * ex / (t2 ** (n + 1.0))-  where-    d = l - m-    n = fromIntegral nn-    ex = exp (- d * x)-    t1 = 1.0 - ex-    t2 = l - m * ex---- | The inverted cumulative probability distribution 'cumulative'. See also--- 'D.ContDistr'.-quantile :: TimeOfOriginDistribution -> Double -> Time-quantile (TOD n' l m) p-  | p >= 0 && p <= 1 =-    -1.0 / d * log (t1 / t2)-  | otherwise =-    error $-      "PointProcess.quantile: p must be in [0,1] range. Got: "-        ++ show p-        ++ "."-  where-    d = l - m-    n = fromIntegral n'-    t1 = l * (1.0 - p ** (1.0 / n))-    t2 = l - p ** (1.0 / n) * m--instance D.ContGen TimeOfOriginDistribution where-  genContVar = D.genContinuous

@@ -1,96 +0,0 @@-{-# LANGUAGE DeriveDataTypeable #-}-{-# LANGUAGE DeriveGeneric #-}---- |---   Module      :  ELynx.Distribution.TimeOfOriginNearCritical---   Description :  Distribution of time of origin for birth and death trees---   Copyright   :  (c) Dominik Schrempf 2018---   License     :  GPL-3.0-or-later------   Maintainer  :  dominik.schrempf@gmail.com---   Stability   :  unstable---   Portability :  portable------ Creation date: Tue Feb 13 13:16:18 2018.------ See Gernhard, T. (2008). The conditioned reconstructed process. Journal of--- Theoretical Biology, 253(4), 769–778. http://doi.org/10.1016/j.jtbi.2008.04.005.------ Distribution of the time of origin for birth and death trees. See corollary 3.3--- in the paper cited above.-module ELynx.Distribution.TimeOfOriginNearCritical-  ( TimeOfOriginNearCriticalDistribution (..),-    cumulative,-    density,-    quantile,-  )-where--import Data.Data-  ( Data,-    Typeable,-  )-import ELynx.Distribution.Types-import GHC.Generics (Generic)-import qualified Statistics.Distribution as D---- | Distribution of the time of origin for a phylogenetic tree evolving under--- the birth and death process and conditioned on observing n leaves today.-data TimeOfOriginNearCriticalDistribution = TONCD-  { -- | Number of leaves of the tree.-    todTN :: Int,-    -- | Birth rate.-    todLa :: Rate,-    -- | Death rate.-    todMu :: Rate-  }-  deriving (Eq, Typeable, Data, Generic)--instance D.Distribution TimeOfOriginNearCriticalDistribution where-  cumulative = cumulative---- | Cumulative distribution function; see Mathematica notebook.-cumulative :: TimeOfOriginNearCriticalDistribution -> Time -> Double-cumulative (TONCD n' l m) t-  | t <= 0 = 0-  | otherwise = t1 + t2-  where-    d = l - m-    n = fromIntegral n'-    t1 = (t * l / (1.0 + t * l)) ** n-    t2 = (n * t * t1) * d / (2.0 * (1.0 + t * l))--instance D.ContDistr TimeOfOriginNearCriticalDistribution where-  density = density-  quantile = quantile---- | The density function Eq. (5).-density :: TimeOfOriginNearCriticalDistribution -> Time -> Double-density (TONCD n' l m) t-  | t < 0 = 0-  | otherwise = nom / den-  where-    n = fromIntegral n'-    nom =-      n * (t * l / (1 + t * l)) ** n * (2 + (3 + n) * t * l - (1 + n) * t * m)-    den = 2 * t * (1 + t * l) ** 2---- | The inverted cumulative probability distribution 'cumulative'. See also--- 'D.ContDistr'.-quantile :: TimeOfOriginNearCriticalDistribution -> Double -> Time-quantile (TONCD n' l m) p-  | p >= 0 && p <= 1 =-    t1 + t2nom / t2den-  | otherwise =-    error $-      "PointProcess.quantile: p must be in [0,1] range. Got: "-        ++ show p-        ++ "."-  where-    n = fromIntegral n'-    t1 = - p ** (1 / n) / ((-1 + p ** (1 / n)) * l)-    t2nom = p ** (2 / n) * (m - l)-    t2den = 2 * (-1 + p ** (1 / n)) ** 2 * l ** 2--instance D.ContGen TimeOfOriginNearCriticalDistribution where-  genContVar = D.genContinuous

@@ -1,22 +0,0 @@--- |--- Module      :  ELynx.Distribution.Types--- Description :  Data types for distributions on trees--- Copyright   :  (c) Dominik Schrempf 2018--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Wed May 16 12:21:57 2018.-module ELynx.Distribution.Types-  ( Time,-    Rate,-  )-where---- | Branch lengths are measured in time.-type Time = Double---- | Birth or death rates.-type Rate = Double

@@ -1,53 +0,0 @@--- |--- Module      :  ELynx.Export.Tree.Newick--- Description :  Export tree objects to Newick format--- Copyright   :  (c) Dominik Schrempf 2020--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Thu Jan 17 13:51:47 2019.------ Some functions are inspired by--- [Biobase.Newick.Import](https://hackage.haskell.org/package/BiobaseNewick).------ See nomenclature in 'ELynx.Data.Tree.Tree'.-module ELynx.Export.Tree.Newick-  ( toNewick,-    toNewickBuilder,-  )-where--import qualified Data.ByteString.Builder as BB-import qualified Data.ByteString.Lazy.Char8 as BL-import Data.List (intersperse)-import ELynx.Data.Tree.Named-import ELynx.Data.Tree.Phylogeny-import ELynx.Data.Tree.Rooted---- | See 'toNewick'.-toNewickBuilder :: Named a => Tree Phylo a -> BB.Builder-toNewickBuilder t = go t <> BB.char8 ';'-  where-    go (Node b l []) = lbl b l-    go (Node b l ts) =-      BB.char8 '('-        <> 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)-    lbl x y =-      BB.lazyByteString (getName y)-        <> mBrLenBuilder x-        -- After reading several discussion, I go for the "more semantical-        -- form" with branch support values in square brackets.-        <> mBrSupBuilder x---- | General conversion of a tree into a Newick 'BL.Bytestring'. Use provided--- functions to extract node labels and branch lengths builder objects. See also--- Biobase.Newick.Export.-toNewick :: Named a => Tree Phylo a -> BL.ByteString-toNewick = BB.toLazyByteString . toNewickBuilder

@@ -1,31 +0,0 @@-{-# LANGUAGE OverloadedStrings #-}---- |--- Module      :  ELynx.Export.Tree.Nexus--- Description :  Export trees to Nexus files--- Copyright   :  (c) Dominik Schrempf 2020--- License     :  GPL-3------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Tue Apr 28 20:24:19 2020.-module ELynx.Export.Tree.Nexus-  ( toNexusTrees,-  )-where--import qualified Data.ByteString.Lazy.Char8 as BL-import ELynx.Data.Tree.Named-import ELynx.Data.Tree.Phylogeny-import ELynx.Data.Tree.Rooted-import ELynx.Export.Nexus-import ELynx.Export.Tree.Newick---- | Export a list of (NAME, TREE) to a Nexus file.-toNexusTrees :: Named a => [(BL.ByteString, Tree Phylo a)] -> BL.ByteString-toNexusTrees ts = toNexus "TREES" (map tree ts)--tree :: Named a => (BL.ByteString, Tree Phylo a) -> BL.ByteString-tree (n, t) = "  TREE " <> n <> " = " <> toNewick t

@@ -1,235 +0,0 @@-{-# LANGUAGE DeriveGeneric #-}---- Module      :  ELynx.Import.Tree.Newick--- Description :  Import Newick trees--- Copyright   :  (c) Dominik Schrempf 2020--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Thu Jan 17 14:56:27 2019.------ Some functions are inspired by--- [Biobase.Newick.Import](https://hackage.haskell.org/package/BiobaseNewick).------ [Specifications](http://evolution.genetics.washington.edu/phylip/newicktree.html)------ In particular, no conversion from _ to (space) is done right now.------ For a description of rooted 'Tree's, please see the 'ELynx.Data.Tree.Rooted'---- |--- module header.-module ELynx.Import.Tree.Newick-  ( NewickFormat (..),-    description,-    newick,-    oneNewick,-    someNewick,-  )-where--import Control.Applicative-import Data.Aeson (FromJSON, ToJSON)-import Data.Attoparsec.ByteString.Char8-import qualified Data.ByteString.Char8 as BS-import ELynx.Data.Tree.Measurable-import ELynx.Data.Tree.Phylogeny-import ELynx.Data.Tree.Rooted hiding (forest, label)-import ELynx.Data.Tree.Supported-import GHC.Generics-import Prelude hiding (takeWhile)---- | Newick tree format.------ >>> unlines $ map (("- " <>) . description) (allValues :: [NewickFormat])--- - Standard: Branch support values are stored in square brackets after branch lengths.--- - IqTree:   Branch support values are stored as node names after the closing bracket of forests.--- - 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.-data NewickFormat = Standard | IqTree | RevBayes-  deriving (Eq, Show, Read, Bounded, Enum, Generic)--instance FromJSON NewickFormat--instance ToJSON NewickFormat---- | Short description of the supported Newick formats.-description :: NewickFormat -> String-description Standard =-  "Standard: Branch support values are stored in square brackets after branch lengths."-description IqTree =-  "IqTree:   Branch support values are stored as node names after the closing bracket of forests."-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 :: 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.-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.-someNewick :: NewickFormat -> Parser (Forest Phylo BS.ByteString)-someNewick Standard = someNewickStandard-someNewick IqTree = someNewickIqTree-someNewick RevBayes = someNewickRevBayes---- 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"---- Parse a single Newick tree. Fails when end of file is not reached.-oneNewickStandard :: Parser (Tree Phylo BS.ByteString)-oneNewickStandard = newickStandard <* endOfInput <?> "oneNewickStandard"---- Parse one ore more Newick trees until end of file.-someNewickStandard :: Parser (Forest Phylo BS.ByteString)-someNewickStandard = some newickStandard <* endOfInput <?> "someNewickStandard"--tree :: Parser (Tree Phylo BS.ByteString)-tree = branched <|> leaf <?> "tree"--branched :: Parser (Tree Phylo BS.ByteString)-branched = (<?> "branched") $ do-  f <- forest-  n <- name-  p <- phylo-  return $ Node p n f---- A 'forest' is a set of trees separated by @,@ and enclosed by parentheses.-forest :: Parser (Forest Phylo BS.ByteString)-forest = char '(' *> (tree `sepBy1` char ',') <* char ')' <?> "forest"---- A 'leaf' has a 'name' and a 'phylo' branch.-leaf :: Parser (Tree Phylo BS.ByteString)-leaf = (<?> "leaf") $ do-  n <- name-  p <- phylo-  return $ Node p n []--nameChar :: Char -> Bool-nameChar c = c `notElem` " :;()[],"---- A name can be any string of printable characters except blanks, colons,--- semicolons, parentheses, and square brackets (and commas).-name :: Parser BS.ByteString-name = takeWhile nameChar <?> "name"--phylo :: Parser Phylo-phylo = Phylo <$> optional branchLength <*> optional branchSupport <?> "phylo"---- Branch length.-branchLength :: Parser BranchLength-branchLength = char ':' *> double <?> "branchLength"--branchSupport :: Parser BranchSupport-branchSupport = (<?> "branchSupport") $-  do-    _ <- char '['-    s <- double-    _ <- char ']'-    return s------------------------------------------------------------------------------------- IQ-TREE.---- IQ-TREE stores the branch support as node names after the closing bracket of--- a forest. Parse a single Newick tree. Also succeeds when more trees follow.-newickIqTree :: Parser (Tree Phylo BS.ByteString)-newickIqTree = skipWhile isSpace *> treeIqTree <* char ';' <* skipWhile isSpace <?> "newickIqTree"---- See 'newickIqTree'. Parse a single Newick tree. Fails when end of file is not--- reached.-oneNewickIqTree :: Parser (Tree Phylo BS.ByteString)-oneNewickIqTree = newickIqTree <* endOfInput <?> "oneNewickIqTree"---- See 'newickIqTree'. Parse one ore more Newick trees until end of file.-someNewickIqTree :: Parser (Forest Phylo BS.ByteString)-someNewickIqTree = some newickIqTree <* endOfInput <?> "someNewickIqTree"---- IQ-TREE stores the branch support as node names after the closing bracket of a forest.-treeIqTree :: Parser (Tree Phylo BS.ByteString)-treeIqTree = branchedIqTree <|> leaf <?> "treeIqTree"---- IQ-TREE stores the branch support as node names after the closing bracket of a forest.-branchedIqTree :: Parser (Tree Phylo BS.ByteString)-branchedIqTree = (<?> "branchedIqTree") $ do-  f <- forestIqTree-  s <- optional double-  n <- name-  b <- optional branchLength-  return $ Node (Phylo b s) n f---- IQ-TREE stores the branch support as node names after the closing bracket of a forest.-forestIqTree :: Parser (Forest Phylo BS.ByteString)-forestIqTree = (<?> "forestIqTree") $ do-  _ <- char '('-  f <- treeIqTree `sepBy1` char ','-  _ <- char ')'-  return f------------------------------------------------------------------------------------- RevBayes.---- RevBayes uses square brackets and key-value pairs to define information--- about nodes and branches. Parse a single Newick tree. Also succeeds when more--- trees follow.------ XXX: Key value pairs are ignored at the moment.-newickRevBayes :: Parser (Tree Phylo BS.ByteString)-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.-oneNewickRevBayes :: Parser (Tree Phylo BS.ByteString)-oneNewickRevBayes = newickRevBayes <* endOfInput <?> "oneNewickRevBayes"---- See 'newickRevBayes'. Parse one ore more Newick trees until end of file.-someNewickRevBayes :: Parser (Forest Phylo BS.ByteString)-someNewickRevBayes = some newickRevBayes <* endOfInput <?> "someNewickRevBayes"--treeRevBayes :: Parser (Tree Phylo BS.ByteString)-treeRevBayes = branchedRevBayes <|> leafRevBayes <?> "treeRevBayes"--branchedRevBayes :: Parser (Tree Phylo BS.ByteString)-branchedRevBayes = (<?> "branchedRevgBayes") $ do-  f <- forestRevBayes-  n <- nameRevBayes-  b <- optional branchLengthRevBayes-  return $ Node (Phylo b Nothing) n f--forestRevBayes :: Parser (Forest Phylo BS.ByteString)-forestRevBayes = (<?> "forestRevBayes") $ do-  _ <- char '('-  f <- treeRevBayes `sepBy1` char ','-  _ <- char ')'-  return f--nameRevBayes :: Parser BS.ByteString-nameRevBayes = name <* optional brackets <?> "nameRevBayes"--branchLengthRevBayes :: Parser BranchLength-branchLengthRevBayes = branchLength <* optional brackets <?> "branchLengthRevBayes"--leafRevBayes :: Parser (Tree Phylo BS.ByteString)-leafRevBayes = (<?> "leafRevBayes") $ do-  n <- nameRevBayes-  b <- optional branchLengthRevBayes-  return $ Node (Phylo b Nothing) n []---- Drop anything between brackets.-brackets :: Parser ()-brackets = (<?> "brackets") $ do-  _ <- char '['-  _ <- takeWhile (/= ']')-  _ <- char ']'-  return ()

@@ -1,45 +0,0 @@-{-# LANGUAGE OverloadedStrings #-}---- |--- Module      :  ELynx.Import.Tree.Nexus--- Description :  Import trees from Nexus files--- Copyright   :  (c) Dominik Schrempf 2020--- License     :  GPL-3------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Tue Apr 28 17:44:13 2020.-module ELynx.Import.Tree.Nexus-  ( nexusTrees,-  )-where--import Control.Applicative-import Data.Attoparsec.ByteString.Char8-import qualified Data.ByteString.Char8 as BS-import ELynx.Data.Tree.Phylogeny-import ELynx.Data.Tree.Rooted-import ELynx.Import.Nexus-import ELynx.Import.Tree.Newick-import Prelude hiding (takeWhile)---- | Parse a Nexus files with a TREES block.-nexusTrees :: NewickFormat -> Parser [(BS.ByteString, Tree Phylo BS.ByteString)]-nexusTrees = nexus . trees--trees :: NewickFormat -> Block [(BS.ByteString, Tree Phylo BS.ByteString)]-trees f = Block "TREES" (some $ namedNewick f)--namedNewick :: NewickFormat -> Parser (BS.ByteString, Tree Phylo BS.ByteString)-namedNewick f = do-  _ <- skipWhile isSpace-  _ <- string "TREE"-  _ <- skipWhile isSpace-  n <- takeWhile1 (\x -> isAlpha_ascii x || isDigit x)-  _ <- skipWhile isSpace-  _ <- char '='-  _ <- skipWhile isSpace-  t <- newick f-  return (n, t)

@@ -1,60 +0,0 @@--- |--- Module      :  ELynx.Simulate.Coalescent--- Description :  Generate coalescent trees--- Copyright   :  (c) Dominik Schrempf 2018--- License     :  GPL-3.0-or-later------ Maintainer  :  dominik.schrempf@gmail.com--- Stability   :  unstable--- Portability :  portable------ Creation date: Wed May 16 13:13:11 2018.-module ELynx.Simulate.Coalescent-  ( simulate,-  )-where--import Control.Monad.Primitive-import ELynx.Data.Tree.Measurable-import ELynx.Data.Tree.Phylogeny-import ELynx.Data.Tree.Rooted-import ELynx.Distribution.CoalescentContinuous-import Statistics.Distribution-import System.Random.MWC---- | Simulate a coalescent tree with @n@ leaves. The branch lengths are in units--- of effective population size.-simulate ::-  (PrimMonad m) =>-  -- | Number of leaves.-  Int ->-  Gen (PrimState m) ->-  m (Tree Length Int)-simulate n = simulate' n 0 trs-  where-    trs = [Node (Length 0) i [] | i <- [0 .. n - 1]]--simulate' ::-  (PrimMonad m) =>-  Int ->-  Int ->-  Forest Length Int ->-  Gen (PrimState m) ->-  m (Tree Length Int)-simulate' n a trs g-  | n <= 0 = error "Cannot construct trees without leaves."-  | n == 1 && length trs /= 1 = error "Too many trees provided."-  | n == 1 && length trs == 1 = return $ head trs-  | otherwise = do-    -- Indices of the leaves to join will be i-1 and i.-    i <- uniformR (1, n - 1) g-    -- The time of the coalescent event.-    t <- genContVar (coalescentDistributionCont n) g-    let trs' = map (applyStem (+ t)) trs -- Move time 't' up on the tree.-        tl = trs' !! (i - 1)-        tr = trs' !! i-        -- Join the two chosen trees.-        tm = Node (Length 0) a [tl, tr]-        -- Take the trees on the left, the merged tree, and the trees on the right.-        trs'' = take (i - 1) trs' ++ [tm] ++ drop (i + 1) trs'-    simulate' (n - 1) a trs'' g

@@ -1,289 +0,0 @@-{-# LANGUAGE BangPatterns #-}---- |---   Module      :  ELynx.Simulate.PointProcess---   Description :  Point process and functions---   Copyright   :  (c) Dominik Schrempf 2018---   License     :  GPL-3.0-or-later------   Maintainer  :  dominik.schrempf@gmail.com---   Stability   :  unstable---   Portability :  portable------ Creation date: Tue Feb 13 13:16:18 2018.------ See Gernhard, T. (2008). The conditioned reconstructed process. Journal of--- Theoretical Biology, 253(4), 769–778. http://doi.org/10.1016/j.jtbi.2008.04.005.------ The point process can be used to simulate reconstructed trees under the birth--- and death process.-module ELynx.Simulate.PointProcess-  ( PointProcess (..),-    TimeSpec,-    simulate,-    toReconstructedTree,-    simulateReconstructedTree,-    simulateNReconstructedTrees,-  )-where--import Control.Monad-import Control.Monad.Primitive-import Data.Function-import Data.List-import Data.Sequence (Seq)-import qualified Data.Sequence as S-import ELynx.Data.Tree.Measurable-import ELynx.Data.Tree.Phylogeny-import ELynx.Data.Tree.Rooted-import ELynx.Distribution.BirthDeath-import ELynx.Distribution.BirthDeathCritical-import ELynx.Distribution.BirthDeathCriticalNoTime-import ELynx.Distribution.BirthDeathNearlyCritical-import ELynx.Distribution.TimeOfOrigin-import ELynx.Distribution.TimeOfOriginNearCritical-import ELynx.Distribution.Types-import qualified Statistics.Distribution as D-  ( genContVar,-  )-import System.Random.MWC---- Require near critical process if birth and death rates are closer than this value.-epsNearCriticalPointProcess :: Double-epsNearCriticalPointProcess = 1e-5---- Also the distribution of origins needs a Tailor expansion for near critical values.------ TODO: Check why the two epsilons are chosen differently.-epsNearCriticalTimeOfOrigin :: Double-epsNearCriticalTimeOfOrigin = 1e-8---- Require critical process if birth and death rates are closer than this value.-eps :: Double-eps = 1e-12--(=~=) :: Double -> Double -> Bool-x =~= y = eps > abs (x - y)---- Sort a list and also return original indices.-sortListWithIndices :: Ord a => [a] -> [(a, Int)]-sortListWithIndices xs = sortBy (compare `on` fst) $ zip xs ([0 ..] :: [Int])---- Insert element into random position of list.-randomInsertList :: PrimMonad m => a -> [a] -> Gen (PrimState m) -> m [a]-randomInsertList e v g = do-  let l = length v-  i <- uniformR (0, l) g-  return $ take i v ++ [e] ++ drop i v---- | A __point process__ for \(n\) points and of age \(t_{or}\) is defined as--- follows. Draw $n$ points on the horizontal axis at \(1,2,\ldots,n\). Pick--- \(n-1\) points at locations \((i+1/2, s_i)\), \(i=1,2,\ldots,n-1\);--- \(0 < s_i < t_{or}\). There is a bijection between (ranked) oriented trees--- and the point process. Usually, a will be 'String' (or 'Int') and b will be--- 'Double'.-data PointProcess a b = PointProcess-  { points :: ![a],-    values :: ![b],-    origin :: !b-  }-  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)---- | Sample a point process using the 'BirthDeathDistribution'. The names of the--- points will be integers.-simulate ::-  (PrimMonad m) =>-  -- | Number of points (samples)-  Int ->-  -- | Time of origin or MRCA-  TimeSpec ->-  -- | Birth rate-  Rate ->-  -- | Death rate-  Rate ->-  -- | Generator (see 'System.Random.MWC')-  Gen (PrimState m) ->-  m (PointProcess Int Double)--- 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-  | -- 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."-  | -- For the critical process, we have no idea about the time of origin, but can-    -- use a specially derived distribution.-    m =~= l =-    do-      !vs <- replicateM (n - 1) (D.genContVar (BDCNTD l) g)-      -- XXX: The length of the root branch will be 0.-      let t = maximum vs-      return $ PointProcess [0 .. (n - 1)] vs t-  | -- For the near critical process, we use a special distribution.-    abs (m - l) <= epsNearCriticalTimeOfOrigin =-    do-      t <- D.genContVar (TONCD n l m) g-      simulate n (Just (t, False)) 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--- 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."-  | -- 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.-    -- -- | m < 0.0   = error "Death rate needs to be positive."-    -- Now, we have three different cases.-    -- 1. The critical branching process.-    -- 2. The near critical branching process.-    -- 3. Normal values :).-    (m =~= l) && not c = 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-    !vs <- replicateM (n - 1) (D.genContVar (BDNCD t l m) g)-    return $ PointProcess [0 .. (n - 1)] vs t-  | not c = do-    !vs <- replicateM (n - 1) (D.genContVar (BDD t l m) g)-    return $ PointProcess [0 .. (n - 1)] vs t-  | (m =~= l) && c = 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-    !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-    !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).-sortPP :: (Ord b) => PointProcess a b -> ([b], [Int])-sortPP (PointProcess _ vs _) = (vsSorted, isSorted)-  where-    vsIsSorted = sortListWithIndices vs-    vsSorted = map fst vsIsSorted-    isSorted = flattenIndices $ map snd vsIsSorted---- Decrement indices that are above the one that is merged.-flattenIndices :: [Int] -> [Int]-flattenIndices is = snd $ mapAccumL fAcc [] is---- TODO: This is the bottleneck for simulating large trees.------ The accumulating function. Count the number of indices which are before the--- current index and lower than the current index.-fAcc :: [Int] -> Int -> ([Int], Int)-fAcc is i = (i : is, i') where i' = i - length (filter (< i) is)---- | See 'simulateReconstructedTree', but n times.-simulateNReconstructedTrees ::-  (PrimMonad m) =>-  -- | Number of trees-  Int ->-  -- | Number of points (samples)-  Int ->-  -- | Time of origin or MRCA-  TimeSpec ->-  -- | Birth rate-  Rate ->-  -- | Death rate-  Rate ->-  -- | Generator (see 'System.Random.MWC')-  Gen (PrimState m) ->-  m (Forest Length Int)-simulateNReconstructedTrees nT nP t l m g-  | nT <= 0 = return []-  | otherwise = replicateM nT $ simulateReconstructedTree nP t l m g---- | Use the point process to simulate a reconstructed tree (see--- 'toReconstructedTree') possibly with specific height and a fixed number of--- leaves according to the birth and death process.-simulateReconstructedTree ::-  (PrimMonad m) =>-  -- | Number of points (samples)-  Int ->-  -- | Time of origin or MRCA-  TimeSpec ->-  -- | Birth rate-  Rate ->-  -- | Death rate-  Rate ->-  -- | Generator (see 'System.Random.MWC')-  Gen (PrimState m) ->-  m (Tree Length Int)-simulateReconstructedTree n t l m g =-  toReconstructedTree 0 <$> simulate n t l m g---- | Convert a point process to a reconstructed tree. See Lemma 2.2.---- Of course, I decided to only use one tree structure with extinct and extant--- leaves (actually a complete tree). So a tree created here just does not--- contain extinct leaves. A function 'isReconstructed' is provided to test if a--- tree is reconstructed (and not complete) in this sense. However, a complete--- tree might show up as "reconstructed", just because, by chance, it does not--- contain extinct leaves. I wanted to use a Monoid constraint to get the unit--- element, but this fails for classical 'Int's. So, I rather have another--- (useless) argument.-toReconstructedTree ::-  a -> -- Default node label.-  PointProcess a Double ->-  Tree Length a-toReconstructedTree l pp@(PointProcess ps vs o)-  | length ps /= length vs + 1 = error "Too few or too many points."-  | length vs <= 1 = error "Too few values."-  | -- -- Test is deactivated.-    -- -- | otherwise = if isReconstructed treeOrigin then treeOrigin else error "Error in algorithm."-    otherwise =-    treeOrigin-  where-    (vsSorted, isSorted) = sortPP pp-    !lvs = S.fromList [Node (Length 0) p [] | p <- ps]-    !heights = S.replicate (length ps) 0-    !treeRoot = toReconstructedTree' isSorted vsSorted l lvs heights-    !h = last vsSorted-    !treeOrigin = applyStem (+ (o - h)) treeRoot---- Move up the tree, connect nodes when they join according to the point process.-toReconstructedTree' ::-  [Int] -> -- Sorted indices, see 'sort'.-  [Double] -> -- Sorted merge values.-  a -> -- Default node label.-  Seq (Tree Length a) -> -- Leaves with accumulated root branch lengths.-  Seq Double -> -- Accumulated heights of the leaves.-  Tree Length a-toReconstructedTree' [] [] _ trs _ = trs `S.index` 0-toReconstructedTree' is vs l trs hs = toReconstructedTree' is' vs' l trs'' hs'-  where-    -- For the algorithm, see 'ELynx.Coalescent.simulate', but index starts-    -- at zero.--    !i = head is-    !is' = tail is-    !v = head vs-    !vs' = tail vs-    -- Left: l, right: r.-    !hl = hs `S.index` i-    !hr = hs `S.index` (i + 1)-    !dvl = v - hl-    !dvr = v - hr-    !tl = applyStem (+ dvl) $ trs `S.index` i-    !tr = applyStem (+ dvr) $ trs `S.index` (i + 1)-    !h' = hl + dvl -- Should be the same as 'hr + dvr'.-    !tm = Node (Length 0) l [tl, tr]-    !trs'' = (S.take i trs S.|> tm) S.>< S.drop (i + 2) trs-    !hs' = (S.take i hs S.|> h') S.>< S.drop (i + 2) hs

@@ -0,0 +1,115 @@+-- |+-- Module      :  ELynx.Topology.Phylogeny+-- Description :  Phylogenetic topologies+-- Copyright   :  (c) Dominik Schrempf, 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Sat Jul 18 13:15:49 2020.+--+-- A topology, as it is used in phylogenetics is a 'Topology' with unique leaf+-- labels, and the order of the topologies in the sub-forest is considered to be+-- meaningless.+--+-- Internally, however, the underlying 'Topology' data structure stores the+-- sub-forest as a (non-empty) list, which has a specific order. Hence, we have+-- to do some tricks when comparing topologies, and topology comparison is slow.+--+-- Also, the uniqueness of the leaves is not ensured by the data type, but has+-- to be checked at runtime. Functions relying on the tree to have unique leaves+-- do perform this check, and return 'Left' with an error message, if the tree+-- has duplicate leaves.+--+-- Note: Topologies are rooted.+--+-- Note: Topologies encoded in Newick format correspond to rooted topologies. By+-- convention only, a topology parsed from Newick format is usually thought to+-- be unrooted, when the root node is multifurcating and has three children.+-- This convention is not enforced here. Newick topologies are just parsed as+-- they are, and a rooted topology is returned.+--+-- The bifurcating root of a topology can be changed with 'roots', or 'rootAt'.+--+-- Topologies with multifurcating root nodes can be properly rooted using+-- 'outgroup'.+module ELynx.Topology.Phylogeny+  ( outgroup,+    roots,+    rootAt,+  )+where++import Data.Set (Set)+import ELynx.Topology.Rooted+import ELynx.Tree.Bipartition++-- TODO.++-- -- | Remove multifurcations.+-- --+-- -- A caterpillar like bifurcating tree is used to resolve all multifurcations on+-- -- a tree. The multifurcating nodes are copied.+-- --+-- -- Branch labels are not handled.+-- resolve :: Tree () a -> Tree () a+-- resolve t@(Node _ _ []) = t+-- resolve (Node _ l [x]) = Node () l [resolve x]+-- resolve (Node _ l [x, y]) = Node () l $ map resolve [x, y]+-- resolve (Node _ l (x : xs)) = Node () l $ map resolve [x, Node () l xs]++-- | Resolve a multifurcation at the root using an outgroup.+outgroup :: Ord a => Set a -> Topology a -> Either String (Topology a)+outgroup = undefined++-- -- | For a rooted tree with a bifurcating root node, get all possible rooted+-- -- trees.+-- --+-- -- The root node is moved.+-- --+-- -- For a tree with @l=2@ leaves, there is one rooted tree. For a bifurcating+-- -- tree with @l>2@ leaves, there are @(2l-3)@ rooted trees. For a general tree+-- -- with a bifurcating root node, and a total number of @n>2@ nodes, there are+-- -- (n-2) rooted trees.+-- --+-- -- Moving a multifurcating root node to another branch would change the+-- -- topology, and so, a bifurcating root is required. To resolve a multifurcating+-- -- root, please see and use TODO.+-- --+-- -- Branch labels are not handled, but see 'rootsBranch'.+-- --+-- -- 'rootAt' roots the tree at a specific position.+-- --+-- -- Return 'Left' if the root node is not 'bifurcating'.+-- roots :: Tree () a -> Either String (Forest () a)+-- roots (Node _ _ []) = Left "roots: Root node is a leaf."+-- roots (Node _ _ [_]) = Left "roots: Root node has degree two."+-- roots t@(Node _ c [tL, tR]) = Right $ t : descend id () c tR tL ++ descend id () c tL tR+-- roots _ = Left "roots: Root node is multifurcating."++-- | For a rooted topology with a bifurcating root node, get all possible rooted+-- topologies.+roots :: Topology a -> Either String (Forest a)+roots = undefined++-- -- | Root a tree at a specific position.+-- --+-- -- Root the tree at the branch defined by the given bipartition. The original+-- -- root node is moved to the new position.+-- --+-- -- The root node must be bifurcating (see 'roots').+-- --+-- -- Branch labels are not handled, but see 'rootAtBranch'.+-- --+-- -- Return 'Left', if:+-- -- - the root node is not bifurcating;+-- -- - the tree has duplicate leaves;+-- -- - the bipartition does not match the leaves of the tree.+-- rootAt :: Ord a => Bipartition a -> Tree () a -> Either String (Tree () a)+-- rootAt = rootAtBranch id++-- | Root a tree at a specific position.+rootAt :: Ord a => Bipartition a -> Topology a -> Either String (Forest a)+rootAt = undefined

@@ -0,0 +1,207 @@+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE DeriveGeneric #-}++-- |+-- Module      :  ELynx.Topology.Rooted+-- Description :  Topologies+-- Copyright   :  (c) Dominik Schrempf, 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Sat Jul 11 10:28:28 2020.+--+-- A 'Topology' differs from a classical rose 'Data.Tree.Tree' in that it does+-- not have internal node labels. The leaves have labels.+--+-- For rooted trees, please see 'ELynx.Tree.Rooted'.+--+-- In phylogenetics, the order of children of a topology node is arbitrary.+-- Internally, however, the underlying 'Topology' data structure stores the+-- sub-forest as a (non-empty) list, which has a specific order. Hence, we have+-- to do some tricks when comparing topologies, and topology comparison is slow.+module ELynx.Topology.Rooted+  ( -- * Data type+    Topology (..),+    Forest,+    fromTree,+    fromLabeledTree,++    -- * Functions+    degree,+    leaves,+    flatten,+    identify,+    prune,+    dropLeavesWith,+    zipTreesWith,+    zipTrees,+    duplicateLeaves,+  )+where++import Control.Applicative+import Control.DeepSeq+import Control.Monad+import Data.Aeson+import Data.Data+import Data.Foldable+import Data.List.NonEmpty (NonEmpty)+import qualified Data.List.NonEmpty as N+import Data.Maybe+import qualified Data.Set as S+import Data.Traversable+import qualified Data.Tree as T+import qualified ELynx.Tree.Rooted as R+import GHC.Generics++singleton :: NonEmpty a -> Bool+singleton xs = 1 == length (N.take 2 xs)++-- | Rooted topologies with leaf labels.+data Topology a+  = Node {forest :: Forest a}+  | Leaf {label :: a}+  deriving (Eq, Read, Show, Data, Generic)++-- | A shortcut.+type Forest a = NonEmpty (Topology a)++instance Functor Topology where+  fmap f (Node ts) = Node $ fmap (fmap f) ts+  fmap f (Leaf lb) = Leaf $ f lb++instance Foldable Topology where+  foldMap f (Node ts) = foldMap (foldMap f) ts+  foldMap f (Leaf lb) = f lb++  null _ = False+  {-# INLINE null #-}++  toList = flatten+  {-# INLINE toList #-}++instance Traversable Topology where+  traverse g (Node ts) = Node <$> traverse (traverse g) ts+  traverse g (Leaf lb) = Leaf <$> g lb++-- TODO: This type checks, but I doubt the implementation is bug-free.+instance Applicative Topology where+  pure = Leaf++  (Node tsF) <*> tx = Node $ fmap (<*> tx) tsF+  (Leaf lbF) <*> tx = lbF <$> tx++  liftA2 f (Node tsX) ty = Node $ fmap (\tx -> liftA2 f tx ty) tsX+  liftA2 f (Leaf lbX) (Node tsY) = Node $ fmap (f lbX <$>) tsY+  liftA2 f (Leaf lbX) (Leaf lbY) = Leaf $ f lbX lbY++  (Node tsX) *> ty@(Node tsY) = Node $ tsY <> fmap (*> ty) tsX+  (Leaf _) *> (Node tsY) = Node tsY+  _ *> (Leaf lbY) = Leaf lbY++  (Node tsX) <* ty = Node $ fmap (<* ty) tsX+  (Leaf lbX) <* _ = Leaf lbX++-- TODO: This type checks, but I doubt the implementation is bug-free.+instance Monad Topology where+  (Node ts) >>= f = Node $ fmap (>>= f) ts+  (Leaf lb) >>= f = case f lb of+    Node ts' -> Node ts'+    Leaf lb' -> Leaf lb'++instance NFData a => NFData (Topology a) where+  rnf (Node ts) = rnf ts+  rnf (Leaf lb) = rnf lb++instance ToJSON a => ToJSON (Topology a)++instance FromJSON a => FromJSON (Topology a)++-- | The degree of the root node.+degree :: Topology a -> Int+degree (Node ts) = (+ 1) $ length ts+degree (Leaf _) = 1++-- | Set of leaves.+leaves :: Ord a => Topology a -> [a]+leaves (Leaf lb) = [lb]+leaves (Node ts) = concatMap leaves ts++-- | Return leaf labels in pre-order.+flatten :: Topology a -> [a]+flatten t = squish t []+  where+    squish (Node ts) xs = foldr squish xs ts+    squish (Leaf lb) xs = lb : xs++-- TODO: Provide and fix tests, provide arbitrary instances.++-- | Convert a rooted rose tree to a rooted topology. Internal node labels are lost.+fromTree :: T.Tree a -> Topology a+fromTree (T.Node lb []) = Leaf lb+fromTree (T.Node _ xs) = Node $ fromTree <$> N.fromList xs++-- | Convert a rooted, labeled rose tree to a rooted topology. Branch labels and+-- internal node labels are lost.+fromLabeledTree :: R.Tree e a -> Topology a+fromLabeledTree (R.Node _ lb []) = Leaf lb+fromLabeledTree (R.Node _ _ xs) = Node $ fromLabeledTree <$> N.fromList xs++-- | Label the leaves with unique integers starting at 0.+identify :: Traversable t => t a -> t Int+identify = snd . mapAccumL (\i _ -> (i + 1, i)) (0 :: Int)++-- | Prune degree two nodes.+prune :: Topology a -> Topology a+prune (Node ts)+  | singleton ts = Node $ fmap prune $ forest $ N.head ts+  | otherwise = Node $ fmap prune ts+prune (Leaf lb) = Leaf lb++-- | Drop leaves satisfying predicate.+--+-- Degree two nodes may arise.+--+-- Return 'Nothing' if all leaves satisfy the predicate.+dropLeavesWith :: (a -> Bool) -> Topology a -> Maybe (Topology a)+dropLeavesWith p (Leaf lb)+  | p lb = Nothing+  | otherwise = Just $ Leaf lb+dropLeavesWith p (Node ts) =+  if null ts'+    then Nothing+    else -- XXX: May be slow, unnecessary conversion to and from list.+      Just $ Node $ N.fromList ts'+  where+    ts' = catMaybes $ N.toList $ fmap (dropLeavesWith p) ts++-- | Zip leaves of two equal topologies.+--+-- Return 'Nothing' if the topologies are different.+zipTreesWith :: (a1 -> a2 -> a) -> Topology a1 -> Topology a2 -> Maybe (Topology a)+zipTreesWith f (Node tsL) (Node tsR) =+  if N.length tsL == N.length tsR+    then -- XXX: May be slow, unnecessary conversion to and from list.+      zipWithM (zipTreesWith f) (N.toList tsL) (N.toList tsR) >>= Just . Node . N.fromList+    else Nothing+zipTreesWith f (Leaf lbL) (Leaf lbR) = Just $ Leaf $ f lbL lbR+zipTreesWith _ _ _ = Nothing++-- | Zip leaves of two equal topologies.+--+-- Return 'Nothing' if the topologies are different.+zipTrees :: Topology a1 -> Topology a2 -> Maybe (Topology (a1, a2))+zipTrees = zipTreesWith (,)++duplicates :: Ord a => [a] -> Bool+duplicates = go S.empty+  where+    go _ [] = False+    go seen (x : xs) = x `S.member` seen || go (S.insert x seen) xs++-- | Check if a topology has duplicate leaves.+duplicateLeaves :: Ord a => Topology a -> Bool+duplicateLeaves = duplicates . leaves

@@ -0,0 +1,82 @@+-- TODO: Topology data type.+-- data Topology a = Node (NonEmptySet (Topology a)) | Leaf a++-- |+-- Module      :  ELynx.Tree+-- Description :  Phylogenetic trees+-- Copyright   :  (c) Dominik Schrempf 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Sat Mar 21 16:27:20 2020.+module ELynx.Tree+  ( -- * Rooted trees+    module ELynx.Tree.Rooted,+    module ELynx.Tree.Zipper,++    -- * Branch label classes+    module ELynx.Tree.Measurable,+    module ELynx.Tree.Splittable,+    module ELynx.Tree.Supported,++    -- * Node label classes+    module ELynx.Tree.Named,++    -- * Phylogenies+    module ELynx.Tree.Phylogeny,++    -- * Partitions and distances+    module ELynx.Tree.Bipartition,+    module ELynx.Tree.Partition,+    module ELynx.Tree.Distance,++    -- * Import and Export+    module ELynx.Tree.Export.Newick,+    module ELynx.Tree.Export.Nexus,+    module ELynx.Tree.Import.Newick,+    module ELynx.Tree.Import.Nexus,+  )+where++import ELynx.Tree.Bipartition+import ELynx.Tree.Distance+import ELynx.Tree.Export.Newick+import ELynx.Tree.Export.Nexus+import ELynx.Tree.Import.Newick+import ELynx.Tree.Import.Nexus+import ELynx.Tree.Measurable+import ELynx.Tree.Named+import ELynx.Tree.Partition+import ELynx.Tree.Phylogeny+import ELynx.Tree.Rooted+import ELynx.Tree.Splittable+import ELynx.Tree.Supported+import ELynx.Tree.Zipper++-- -- | An evolutionary label has some information about where the corresponding+-- -- node is on the tree, and if the node is 'extant', 'extinct', 'internal', or+-- -- 'external'. The latter two could also be determined from the tree. This could+-- -- be species, genes or individuals; probably more.+-- class EvoLabel n where+--   extant          :: n -> Bool+--   extinct         :: n -> Bool++--   internal        :: n -> Bool+--   internal n = not $ extant n || extinct n++--   external        :: n -> Bool+--   external   = not . internal++-- -- -- | Glue branches together, so that one new tree emerges. It's root node is+-- -- -- new, the sub-forest has to be given (a list of trees).+-- -- glue :: (NodeType c)+-- --      => PhyloLabel a b c       -- ^ New root node.+-- --      -> [PhyloTree a b c]      -- ^ Sub-forest.+-- --      -> PhyloTree a b c+-- -- glue s@(PhyloLabel _ _ n) ts+-- --   | extant n  = error "Root node cannot be of type 'Exant'."+-- --   | extinct n = error "Root node cannot be of type 'Extinct'."+-- --   | otherwise = Node s ts

@@ -0,0 +1,198 @@+{-# LANGUAGE GeneralizedNewtypeDeriving #-}++-- |+-- Module      :  ELynx.Tree.Bipartition+-- Description :  Bipartitions on trees+-- Copyright   :  (c) Dominik Schrempf 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Fri Aug 30 15:28:17 2019.+--+-- 'Bipartition's are weird in that+-- > Bipartition x y == Bipartition y x+-- is True.+--+-- Also,+-- > Bipartition x y > Bipartition y x+-- is False, even when @x > y@.+--+-- That's why we have to make sure that for+-- > Bipartition x y+-- we always have @x >= y@.+module ELynx.Tree.Bipartition+  ( groups,++    -- * Data type+    Bipartition (fromBipartition),+    bp,+    bpUnsafe,+    toSet,+    bpHuman,++    -- * Work with 'Bipartition's+    bipartition,+    bipartitions,+    getComplementaryLeaves,+    bipartitionToBranch,+  )+where++import Control.Comonad+import Control.DeepSeq+import Data.List hiding (partition)+import Data.Map (Map)+import qualified Data.Map as M+import Data.Set (Set)+import qualified Data.Set as S+import ELynx.Tree.Rooted++-- | Each node of a tree is root of an induced subtree. Set the node labels to+-- the leaves of the induced subtrees.+groups :: Tree e a -> Tree e [a]+-- 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.+--+-- 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.+newtype Bipartition a = Bipartition+  { fromBipartition :: (Set a, Set a)+  }+  deriving (Eq, Ord, Show, Read, NFData)++-- | Create a bipartition from two sets.+--+-- Ensure that the smaller set comes first.+--+-- Return 'Left' if one set is empty.+bp :: Ord a => Set a -> Set a -> Either String (Bipartition a)+bp xs ys+  | S.null xs = Left "bp: Left set empty."+  | S.null ys = Left "bp: Right set empty."+  | otherwise = Right $ bpUnsafe xs ys++-- | Create a bipartition from two sets.+--+-- Ensure that the smaller set comes first.+bpUnsafe :: Ord a => Set a -> Set a -> Bipartition a+bpUnsafe xs ys = if xs >= ys then Bipartition (xs, ys) else Bipartition (ys, xs)++-- | Conversion to a set containing both partitions.+toSet :: Ord a => Bipartition a -> Set a+toSet (Bipartition (x, y)) = S.union x y++-- I decided not to provide a human readable show instance because I need the+-- following identity to hold:+--+-- > read . show = id+--+-- This identity is met by the derived instance anyways. A more human readable+-- instance would most likely violate the identity.++-- | Show a bipartition in a human readable format. Use a provided function to+-- extract information of interest.+bpHuman :: Show a => Bipartition a -> String+bpHuman (Bipartition (x, y)) = "(" ++ setShow x ++ "|" ++ setShow y ++ ")"++-- Show the elements of a set in a human readable format.+setShow :: Show a => Set a -> String+setShow = intercalate "," . map show . S.toList++-- -- | Map a function over all elements in the 'Bipartition'.+-- bpMap :: Ord b => (a -> b) -> Bipartition a -> Bipartition b+-- bpMap f (Bipartition (x, y)) = bp (S.map f x) (S.map f y)++-- | For a bifurcating root, get the bipartition induced by the root node.+--+-- Return 'Left' if+-- - the root node is not bifurcating;+-- - a leave set is empty.+bipartition :: Ord a => Tree e a -> Either String (Bipartition a)+bipartition (Node _ _ [x, y]) = bp (S.fromList $ leaves x) (S.fromList $ leaves y)+bipartition _ = Left "bipartition: Root node is not bifurcating."++-- | Get all bipartitions of the tree.+--+-- Return 'Left' if the tree contains duplicate leaves.+bipartitions :: Ord a => Tree e a -> Either String (Set (Bipartition a))+bipartitions t+  | duplicateLeaves t = Left "bipartitions: Tree contains duplicate leaves."+  | otherwise = Right $ bipartitions' S.empty $ S.fromList <$> groups t++-- | Report the complementary leaves for each child.+getComplementaryLeaves ::+  (Ord a) =>+  -- Complementary leaves.+  Set a ->+  -- Tree with node labels storing leaves.+  Tree e (Set a) ->+  [Set a]+getComplementaryLeaves p (Node _ _ ts) =+  [ S.unions $ p : take i lvsChildren ++ drop (i + 1) lvsChildren+    | i <- [0 .. (n -1)]+  ]+  where+    n = length ts+    lvsChildren = map label ts++-- See 'bipartitions', but do not check if leaves are unique, nor if+-- bipartitions are valid.+bipartitions' :: Ord a => Set a -> Tree e (Set a) -> Set (Bipartition a)+bipartitions' p (Node _ p' []) = either (const S.empty) S.singleton $ bp p p'+bipartitions' p t@(Node _ p' ts) =+  S.unions $+    either (const S.empty) S.singleton (bp p p') :+      [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'.+--+-- Since the induced bipartitions of the daughter branches of a bifurcating root+-- node are equal, the branches leading to the root have to be combined in this+-- case. See http://evolution.genetics.washington.edu/phylip/doc/treedist.html+-- and how unrooted trees should be handled.+--+-- Further, branches connected to degree two nodes also induce the same+-- bipartitions and have to be combined.+--+-- For combining branches, a binary function is required. This requirement is+-- encoded in the 'Semigroup' type class constraint (see 'prune').+--+-- Return 'Left' if the tree contains duplicate leaves.+bipartitionToBranch ::+  (Semigroup e, Ord a) =>+  Tree e a ->+  Either String (Map (Bipartition a) e)+bipartitionToBranch t+  | duplicateLeaves t = Left "bipartitionToBranch: Tree contains duplicate leaves."+  | otherwise = Right $ bipartitionToBranch' S.empty pTree+  where+    pTree = S.fromList <$> groups t++-- When calculating the map, branches separated by various degree two nodes have+-- to be combined. Hence, not only the complementary leaves, but also the branch+-- label itself have to be passed along.+bipartitionToBranch' ::+  (Semigroup e, Ord a) =>+  -- Complementary leaves.+  Set a ->+  -- Partition tree.+  Tree e (Set a) ->+  Map (Bipartition a) e+bipartitionToBranch' p t@(Node b p' ts) =+  M.unionsWith (<>) $+    either (const M.empty) (`M.singleton` b) (bp p p') :+      [bipartitionToBranch' c s | (c, s) <- zip cs ts]+  where+    cs = getComplementaryLeaves p t

@@ -0,0 +1,139 @@+-- |+-- Module      :  ELynx.Tree.Distance+-- Description :  Compute distances between trees+-- Copyright   :  (c) Dominik Schrempf 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Thu Jun 13 17:15:54 2019.+--+-- Various distance functions for phylogenetic trees (and trees with branch+-- lengths in general).+--+-- 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.+module ELynx.Tree.Distance+  ( symmetric,+    incompatibleSplits,+    branchScore,+  )+where++-- adjacent,++import Data.Bifunctor+import Data.List+import qualified Data.Map as M+import Data.Monoid+import Data.Set (Set)+import qualified Data.Set as S+import ELynx.Tree.Bipartition+import ELynx.Tree.Measurable+import ELynx.Tree.Partition+import ELynx.Tree.Rooted++-- Symmetric difference between two 'Set's.+symmetricDifference :: Ord a => Set a -> Set a -> Set a+symmetricDifference xs ys = S.difference xs ys `S.union` S.difference ys xs++-- | Symmetric (Robinson-Foulds) distance between two trees.+--+-- Although a rooted tree data type is used, the distance between the unrooted+-- trees is returned.+--+-- Return 'Nothing' if the trees contain different leaves.+--+-- XXX: Comparing a list of trees may recompute bipartitions.+symmetric :: Ord a => Tree e1 a -> Tree e2 a -> Either String Int+symmetric t1 t2+  | S.fromList (leaves t1) /= S.fromList (leaves t2) = Left "symmetric: Trees contain different leaves."+  | otherwise = do+    bps1 <- bipartitions t1+    bps2 <- bipartitions t2+    return $ length $ symmetricDifference bps1 bps2++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++-- | Number of incompatible splits.+--+-- Similar to 'symmetric' but all bipartitions induced by multifurcations are+-- considered. For a detailed description of how the distance is calculated, see+-- 'ELynx.Tree.Bipartition.bipartitionCompatible'.+--+-- A multifurcation on a tree may (but not necessarily does) represent missing+-- information about the order of bifurcations. In this case, it is interesting+-- to get a set of compatible bifurcations of the tree. For example, the star tree+--+-- > (A,B,C,D);+--+-- induces the following bipartitions:+--+-- > A|BCD+-- > B|ACD+-- > C|ABD+-- > D|ABC+--+-- However, the tree is additionally compatible with the following hidden+-- bipartitions:+--+-- > AB|CD+-- > AC|BD+-- > AD|BC+--+-- For an explanation of how compatibility of partitions is checked, see+-- 'compatible'. Before using 'compatible', bipartitions are simply converted to+-- partitions with two subsets.+--+-- A bipartition is incompatible with a tree if it is incompatible with all+-- induced multifurcations of the tree.+--+-- XXX: Comparing a list of trees with this function recomputes bipartitions.+incompatibleSplits :: (Show a, Ord a) => Tree e1 a -> Tree e2 a -> Either String Int+incompatibleSplits t1 t2+  | S.fromList (leaves t1) /= S.fromList (leaves t2) =+    Left "incompatibleSplits: Trees do not have equal leaf sets."+  | otherwise = do+    -- Bipartitions.+    bs1 <- bipartitions t1+    bs2 <- bipartitions t2+    -- traceShowM $ "bs1" ++ show (S.map bpHuman bs1)+    -- traceShowM $ "bs2" ++ show (S.map bpHuman bs2)+    let -- Putative incompatible bipartitions of trees one and two, respectively.+        putIncBs1 = bs1 S.\\ bs2+        putIncBs2 = bs2 S.\\ bs1+    -- Partitions.+    ms1 <- partitions t1+    ms2 <- partitions t2+    -- traceShowM $ "putIncBs1 " ++ show (S.map bpHuman putIncBs1)+    -- traceShowM $ "putIncBs2 " ++ show (S.map bpHuman putIncBs2)+    return $ countIncompatibilities putIncBs1 ms2 + countIncompatibilities putIncBs2 ms1++-- | Compute branch score distance between two trees.+--+-- Although a rooted tree data type is used, the distance between the unrooted+-- trees is returned.+--+-- XXX: Comparing a list of trees with this function recomputes bipartitions.+branchScore :: (Measurable e1, Measurable e2, Ord a) => Tree e1 a -> Tree e2 a -> Either String Double+branchScore t1 t2+  | S.fromList (leaves t1) /= S.fromList (leaves t2) = Left "branchScoreWith: Trees do not have equal leaf sets."+  | otherwise = do+    bpToBr1 <- bipartitionToBranch $ first (Sum . getLen) t1+    bpToBr2 <- bipartitionToBranch $ first (Sum . getLen) t2+    let dBs = M.unionWith (-) bpToBr1 bpToBr2+        dsSquared = foldl' (\acc e -> acc + e * e) 0 dBs+    return $ sqrt $ getSum dsSquared

@@ -0,0 +1,95 @@+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE DeriveGeneric #-}++-- |+--   Module      :  ELynx.Tree.Distribution.BirthDeath+--   Description :  Birth and death distribution+--   Copyright   :  (c) Dominik Schrempf 2018+--   License     :  GPL-3.0-or-later+--+--   Maintainer  :  dominik.schrempf@gmail.com+--   Stability   :  unstable+--   Portability :  portable+--+-- Creation date: Tue Feb 13 13:16:18 2018.+--+-- See Gernhard, T. (2008). The conditioned reconstructed process. Journal of+-- Theoretical Biology, 253(4), 769–778. http://doi.org/10.1016/j.jtbi.2008.04.005.+--+-- Distribution of the values of the point process such that it corresponds to+-- reconstructed trees under the birth and death process.+module ELynx.Tree.Distribution.BirthDeath+  ( BirthDeathDistribution (..),+    cumulative,+    density,+    quantile,+  )+where++import Data.Data+  ( Data,+    Typeable,+  )+import ELynx.Tree.Distribution.Types+import GHC.Generics (Generic)+import qualified Statistics.Distribution as D++-- | Distribution of the values of the point process such that it corresponds to+-- a reconstructed tree of the birth and death process.+data BirthDeathDistribution = BDD+  { -- | Time to origin of the tree.+    bddTOr :: Time,+    -- | Birth rate.+    bddLa :: Rate,+    -- | Death rate.+    bddMu :: Rate+  }+  deriving (Eq, Typeable, Data, Generic)++instance D.Distribution BirthDeathDistribution where+  cumulative = cumulative++-- | Cumulative distribution function Eq. (3).+cumulative :: BirthDeathDistribution -> Time -> Double+cumulative (BDD t l m) x+  | x <= 0 = 0+  | x > t = 1+  | otherwise = t1 * t2+  where+    d = l - m+    t1 = (1.0 - exp (- d * x)) / (l - m * exp (- d * x))+    t2 = (l - m * exp (- d * t)) / (1.0 - exp (- d * t))++instance D.ContDistr BirthDeathDistribution where+  density = density+  quantile = quantile++-- | Density function Eq. (2).+density :: BirthDeathDistribution -> Time -> Double+density (BDD t l m) x+  | x < 0 = 0+  | x > t = 0+  | otherwise = d ** 2 * t1 * t2+  where+    d = l - m+    t1 = exp (- d * x) / ((l - m * exp (- d * x)) ** 2)+    t2 = (l - m * exp (- d * t)) / (1.0 - exp (- d * t))++-- | Inverted cumulative probability distribution 'cumulative'. See also+-- 'D.ContDistr'.+quantile :: BirthDeathDistribution -> Double -> Time+quantile (BDD t l m) p+  | p >= 0 && p <= 1 =+    res+  | otherwise =+    error $+      "PointProcess.quantile: p must be in range [0,1] but got "+        ++ show p+        ++ "."+  where+    d = l - m+    t2 = (l - m * exp (- d * t)) / (1.0 - exp (- d * t))+    res = (-1.0 / d) * log ((1.0 - p * l / t2) / (1.0 - p * m / t2))++instance D.ContGen BirthDeathDistribution where+  genContVar = D.genContinuous

@@ -0,0 +1,84 @@+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE DeriveGeneric #-}++-- |+--   Module      :  ELynx.Tree.Distribution.BirthDeathCritical+--   Description :  Birth and death distribution+--   Copyright   :  (c) Dominik Schrempf 2018+--   License     :  GPL-3.0-or-later+--+--   Maintainer  :  dominik.schrempf@gmail.com+--   Stability   :  unstable+--   Portability :  portable+--+-- Creation date: Tue Feb 13 13:16:18 2018.+--+-- See Gernhard, T. (2008). The conditioned reconstructed process. Journal of+-- Theoretical Biology, 253(4), 769–778. http://doi.org/10.1016/j.jtbi.2008.04.005.+--+-- Distribution of the values of the point process such that it corresponds to+-- reconstructed trees under the birth and death process; critical birth and death+-- process with lambda=mu.+module ELynx.Tree.Distribution.BirthDeathCritical+  ( BirthDeathCriticalDistribution (..),+    cumulative,+    density,+    quantile,+  )+where++import Data.Data+  ( Data,+    Typeable,+  )+import ELynx.Tree.Distribution.Types+import GHC.Generics (Generic)+import qualified Statistics.Distribution as D++-- | Distribution of the values of the point process such that it corresponds to+-- a reconstructed tree of the birth and death process.+data BirthDeathCriticalDistribution = BDCD+  { -- | Time to origin of the tree.+    bdcdTOr :: Time,+    -- | Birth and death rate.+    bdcdLa :: Rate+  }+  deriving (Eq, Typeable, Data, Generic)++instance D.Distribution BirthDeathCriticalDistribution where+  cumulative = cumulative++-- | Cumulative distribution function section 2.1.2, second formula.+cumulative :: BirthDeathCriticalDistribution -> Time -> Double+cumulative (BDCD t l) x+  | x <= 0 = 0+  | x > t = 1+  | otherwise = x / (1.0 + l * x) * (1.0 + l * t) / t++instance D.ContDistr BirthDeathCriticalDistribution where+  density = density+  quantile = quantile++-- | Density function section 2.1.2, first formula.+density :: BirthDeathCriticalDistribution -> Time -> Double+density (BDCD t l) x+  | x < 0 = 0+  | x > t = 0+  | otherwise = (1.0 + l * t) / (t * (1.0 + l * x) ** 2)++-- | Inverted cumulative probability distribution 'cumulative'. See also+-- 'D.ContDistr'.+quantile :: BirthDeathCriticalDistribution -> Double -> Time+quantile (BDCD t l) p+  | p >= 0 && p <= 1 =+    res+  | otherwise =+    error $+      "PointProcess.quantile: p must be in [0,1] range. Got: "+        ++ show p+        ++ "."+  where+    res = p * t / (1 + l * t - l * p * t)++instance D.ContGen BirthDeathCriticalDistribution where+  genContVar = D.genContinuous

@@ -0,0 +1,76 @@+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE DeriveGeneric #-}++-- |+--   Module      :  ELynx.Tree.Distribution.BirthDeathCriticalNoTime+--   Description :  Birth and death distribution+--   Copyright   :  (c) Dominik Schrempf 2018+--   License     :  GPL-3.0-or-later+--+--   Maintainer  :  dominik.schrempf@gmail.com+--   Stability   :  unstable+--   Portability :  portable+--+-- Creation date: Tue Feb 13 13:16:18 2018.+--+-- See Gernhard, T. (2008). The conditioned reconstructed process. Journal of+-- Theoretical Biology, 253(4), 769–778. http://doi.org/10.1016/j.jtbi.2008.04.005.+--+-- Distribution of the values of the point process such that it corresponds to+-- reconstructed trees under the birth and death process; critical birth and death+-- process with lambda=mu.+module ELynx.Tree.Distribution.BirthDeathCriticalNoTime+  ( BirthDeathCriticalNoTimeDistribution (..),+    cumulative,+    density,+    quantile,+  )+where++import Data.Data+  ( Data,+    Typeable,+  )+import ELynx.Tree.Distribution.Types+import GHC.Generics (Generic)+import qualified Statistics.Distribution as D++-- | Distribution of the values of the point process such that it corresponds to+-- a reconstructed tree of the birth and death process.+newtype BirthDeathCriticalNoTimeDistribution = BDCNTD+  { -- | Birth and death rate.+    bdcntdLa :: Rate+  }+  deriving (Eq, Typeable, Data, Generic)++instance D.Distribution BirthDeathCriticalNoTimeDistribution where+  cumulative = cumulative++-- | Cumulative distribution function section 2.1.2, second formula.+cumulative :: BirthDeathCriticalNoTimeDistribution -> Time -> Double+cumulative (BDCNTD l) x+  | x <= 0 = 0+  | otherwise = x * l / (1.0 + x * l)++instance D.ContDistr BirthDeathCriticalNoTimeDistribution where+  density = density+  quantile = quantile++-- | Density function section 2.1.2, first formula; t cancels out because it is+-- expected to be much larger than 1.0; because t \in [0, \infty].+density :: BirthDeathCriticalNoTimeDistribution -> Time -> Double+density (BDCNTD l) x+  | x < 0 = 0+  | otherwise = l / ((1.0 + x * l) ** 2)++-- | Inverted cumulative probability distribution 'cumulative'. See also+-- 'D.ContDistr'.+quantile :: BirthDeathCriticalNoTimeDistribution -> Double -> Time+quantile (BDCNTD l) p+  | p >= 0 && p <= 1 =+    p / (l - l * p)+  | otherwise =+    error $ "PointProcess.quantile: p must be in [0,1]. Got: " ++ show p ++ "."++instance D.ContGen BirthDeathCriticalNoTimeDistribution where+  genContVar = D.genContinuous

@@ -0,0 +1,109 @@+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE DeriveGeneric #-}++-- |+--   Module      :  ELynx.Tree.Distribution.BirthDeathNearlyCritical+--   Description :  Birth and death distribution+--   Copyright   :  (c) Dominik Schrempf 2018+--   License     :  GPL-3.0-or-later+--+--   Maintainer  :  dominik.schrempf@gmail.com+--   Stability   :  unstable+--   Portability :  portable+--+-- Creation date: Tue Feb 13 13:16:18 2018.+--+-- See Gernhard, T. (2008). The conditioned reconstructed process. Journal of+-- Theoretical Biology, 253(4), 769–778. http://doi.org/10.1016/j.jtbi.2008.04.005.+--+-- Distribution of the values of the point process such that it corresponds to+-- reconstructed trees under the birth and death process; nearly critical birth and+-- death process with lambda~mu.+--+-- Basically, this is a Taylor expansion of Eq. (2) and Eq. (3).+module ELynx.Tree.Distribution.BirthDeathNearlyCritical+  ( BirthDeathNearlyCriticalDistribution (..),+    cumulative,+    density,+    quantile,+  )+where++import Data.Data+  ( Data,+    Typeable,+  )+import ELynx.Tree.Distribution.Types+import GHC.Generics (Generic)+import qualified Statistics.Distribution as D++-- | Distribution of the values of the point process such that it corresponds to+-- a reconstructed tree of the birth and death process.+data BirthDeathNearlyCriticalDistribution = BDNCD+  { -- | Time to origin of the tree.+    bdncdTOr :: Time,+    -- | Birth and death rate.+    bdncdLa :: Rate,+    -- | Birth and death rate.+    bdncdMu :: Rate+  }+  deriving (Eq, Typeable, Data, Generic)++instance D.Distribution BirthDeathNearlyCriticalDistribution where+  cumulative = cumulative++-- | Cumulative distribution function section 2.1.2, second formula.+cumulative :: BirthDeathNearlyCriticalDistribution -> Time -> Double+cumulative (BDNCD t l m) s+  | s <= 0 = 0+  | s > t = 1+  | otherwise = o0 + o1+  where+    o0 = s * (1.0 + t * l) / t / (1.0 + s * l)+    o1 = (- s * s + s * t) * (m - l) / (2.0 * t * (1.0 + s * l) ** 2)++instance D.ContDistr BirthDeathNearlyCriticalDistribution where+  density = density+  quantile = quantile++-- | Density function section 2.1.2, first formula.+density :: BirthDeathNearlyCriticalDistribution -> Time -> Double+density (BDNCD t l m) s+  | s < 0 = 0+  | s > t = 0+  | otherwise = o0 + o1+  where+    o0 = (1.0 + t * l) / (t * (1.0 + s * l) ** 2)+    o1 = (-2.0 * s + t - s * t * l) * (m - l) / (2.0 * t * (1.0 + s * l) ** 3)++-- | Inverted cumulative probability distribution 'cumulative'. See also+-- 'D.ContDistr'.+quantile :: BirthDeathNearlyCriticalDistribution -> Double -> Time+quantile (BDNCD t l m) p+  | p >= 0 && p <= 1 =+    res+  | otherwise =+    error $+      "PointProcess.quantile: p must be in [0,1] range. Got: "+        ++ show p+        ++ "."+  where+    den = l * (-3.0 + 2.0 * t * (-1.0 + p) * l) + m+    t1 = (2.0 + t * (l - 4.0 * p * l + m)) / den+    t2Nom =+      4.0+        + t+        * ( l+              * (4.0 + t * l + 8.0 * p * (1.0 + t * l))+              + 2.0+              * (2.0 + t * l - 4.0 * p * (1.0 + t * l))+              * m+              + t+              * m+              * m+          )+    t2 = t2Nom / (den ** 2)+    res = 0.5 * (t1 + sqrt t2)++instance D.ContGen BirthDeathNearlyCriticalDistribution where+  genContVar = D.genContinuous

@@ -0,0 +1,26 @@+-- |+-- Module      :  ELynx.Tree.Distribution.CoalescentContinuous+-- Description :  Distribution of coalescent times+-- Copyright   :  (c) Dominik Schrempf 2018+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Wed May 16 12:40:45 2018.+module ELynx.Tree.Distribution.CoalescentContinuous+  ( coalescentDistributionCont,+  )+where++import Numeric.SpecFunctions (choose)+import Statistics.Distribution.Exponential++-- | Distribution of the next coalescent event for a number of samples @n@. The+-- time is measured in units of effective number of population size.+coalescentDistributionCont ::+  -- | Sample size.+  Int ->+  ExponentialDistribution+coalescentDistributionCont n = exponential (choose n 2)

@@ -0,0 +1,95 @@+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE DeriveGeneric #-}++-- |+--   Module      :  ELynx.Tree.Distribution.TimeOfOrigin+--   Description :  Distribution of time of origin for birth and death trees+--   Copyright   :  (c) Dominik Schrempf 2018+--   License     :  GPL-3.0-or-later+--+--   Maintainer  :  dominik.schrempf@gmail.com+--   Stability   :  unstable+--   Portability :  portable+--+-- Creation date: Tue Feb 13 13:16:18 2018.+--+-- See Gernhard, T. (2008). The conditioned reconstructed process. Journal of+-- Theoretical Biology, 253(4), 769–778. http://doi.org/10.1016/j.jtbi.2008.04.005.+--+-- Distribution of the time of origin for birth and death trees. See corollary 3.3+-- in the paper cited above.+module ELynx.Tree.Distribution.TimeOfOrigin+  ( TimeOfOriginDistribution (..),+    cumulative,+    density,+    quantile,+  )+where++import Data.Data+  ( Data,+    Typeable,+  )+import ELynx.Tree.Distribution.Types+import GHC.Generics (Generic)+import qualified Statistics.Distribution as D++-- | Distribution of the time of origin for a phylogenetic tree evolving under+-- the birth and death process and conditioned on observing n leaves today.+data TimeOfOriginDistribution = TOD+  { -- | Number of leaves of the tree.+    todTN :: Int,+    -- | Birth rate.+    todLa :: Rate,+    -- | Death rate.+    todMu :: Rate+  }+  deriving (Eq, Typeable, Data, Generic)++instance D.Distribution TimeOfOriginDistribution where+  cumulative = cumulative++-- | Cumulative distribution function Corollary 3.3.+cumulative :: TimeOfOriginDistribution -> Time -> Double+cumulative (TOD n l m) x+  | x <= 0 = 0+  | otherwise = te ** fromIntegral n+  where+    d = l - m+    te = l * (1.0 - exp (- d * x)) / (l - m * exp (- d * x))++instance D.ContDistr TimeOfOriginDistribution where+  density = density+  quantile = quantile++-- | The density function Eq. (5).+density :: TimeOfOriginDistribution -> Time -> Double+density (TOD nn l m) x+  | x < 0 = 0+  | otherwise = n * l ** n * d ** 2 * t1 ** (n - 1.0) * ex / (t2 ** (n + 1.0))+  where+    d = l - m+    n = fromIntegral nn+    ex = exp (- d * x)+    t1 = 1.0 - ex+    t2 = l - m * ex++-- | The inverted cumulative probability distribution 'cumulative'. See also+-- 'D.ContDistr'.+quantile :: TimeOfOriginDistribution -> Double -> Time+quantile (TOD n' l m) p+  | p >= 0 && p <= 1 =+    -1.0 / d * log (t1 / t2)+  | otherwise =+    error $+      "PointProcess.quantile: p must be in [0,1] range. Got: "+        ++ show p+        ++ "."+  where+    d = l - m+    n = fromIntegral n'+    t1 = l * (1.0 - p ** (1.0 / n))+    t2 = l - p ** (1.0 / n) * m++instance D.ContGen TimeOfOriginDistribution where+  genContVar = D.genContinuous

@@ -0,0 +1,96 @@+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE DeriveGeneric #-}++-- |+--   Module      :  ELynx.Tree.Distribution.TimeOfOriginNearCritical+--   Description :  Distribution of time of origin for birth and death trees+--   Copyright   :  (c) Dominik Schrempf 2018+--   License     :  GPL-3.0-or-later+--+--   Maintainer  :  dominik.schrempf@gmail.com+--   Stability   :  unstable+--   Portability :  portable+--+-- Creation date: Tue Feb 13 13:16:18 2018.+--+-- See Gernhard, T. (2008). The conditioned reconstructed process. Journal of+-- Theoretical Biology, 253(4), 769–778. http://doi.org/10.1016/j.jtbi.2008.04.005.+--+-- Distribution of the time of origin for birth and death trees. See corollary 3.3+-- in the paper cited above.+module ELynx.Tree.Distribution.TimeOfOriginNearCritical+  ( TimeOfOriginNearCriticalDistribution (..),+    cumulative,+    density,+    quantile,+  )+where++import Data.Data+  ( Data,+    Typeable,+  )+import ELynx.Tree.Distribution.Types+import GHC.Generics (Generic)+import qualified Statistics.Distribution as D++-- | Distribution of the time of origin for a phylogenetic tree evolving under+-- the birth and death process and conditioned on observing n leaves today.+data TimeOfOriginNearCriticalDistribution = TONCD+  { -- | Number of leaves of the tree.+    todTN :: Int,+    -- | Birth rate.+    todLa :: Rate,+    -- | Death rate.+    todMu :: Rate+  }+  deriving (Eq, Typeable, Data, Generic)++instance D.Distribution TimeOfOriginNearCriticalDistribution where+  cumulative = cumulative++-- | Cumulative distribution function; see Mathematica notebook.+cumulative :: TimeOfOriginNearCriticalDistribution -> Time -> Double+cumulative (TONCD n' l m) t+  | t <= 0 = 0+  | otherwise = t1 + t2+  where+    d = l - m+    n = fromIntegral n'+    t1 = (t * l / (1.0 + t * l)) ** n+    t2 = (n * t * t1) * d / (2.0 * (1.0 + t * l))++instance D.ContDistr TimeOfOriginNearCriticalDistribution where+  density = density+  quantile = quantile++-- | The density function Eq. (5).+density :: TimeOfOriginNearCriticalDistribution -> Time -> Double+density (TONCD n' l m) t+  | t < 0 = 0+  | otherwise = nom / den+  where+    n = fromIntegral n'+    nom =+      n * (t * l / (1 + t * l)) ** n * (2 + (3 + n) * t * l - (1 + n) * t * m)+    den = 2 * t * (1 + t * l) ** 2++-- | The inverted cumulative probability distribution 'cumulative'. See also+-- 'D.ContDistr'.+quantile :: TimeOfOriginNearCriticalDistribution -> Double -> Time+quantile (TONCD n' l m) p+  | p >= 0 && p <= 1 =+    t1 + t2nom / t2den+  | otherwise =+    error $+      "PointProcess.quantile: p must be in [0,1] range. Got: "+        ++ show p+        ++ "."+  where+    n = fromIntegral n'+    t1 = - p ** (1 / n) / ((-1 + p ** (1 / n)) * l)+    t2nom = p ** (2 / n) * (m - l)+    t2den = 2 * (-1 + p ** (1 / n)) ** 2 * l ** 2++instance D.ContGen TimeOfOriginNearCriticalDistribution where+  genContVar = D.genContinuous

@@ -0,0 +1,22 @@+-- |+-- Module      :  ELynx.Tree.Distribution.Types+-- Description :  Data types for distributions on trees+-- Copyright   :  (c) Dominik Schrempf 2018+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Wed May 16 12:21:57 2018.+module ELynx.Tree.Distribution.Types+  ( Time,+    Rate,+  )+where++-- | Branch lengths are measured in time.+type Time = Double++-- | Birth or death rates.+type Rate = Double

@@ -0,0 +1,53 @@+-- |+-- Module      :  ELynx.Tree.Export.Newick+-- Description :  Export tree objects to Newick format+-- Copyright   :  (c) Dominik Schrempf 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Thu Jan 17 13:51:47 2019.+--+-- Some functions are inspired by+-- [Biobase.Newick.Import](https://hackage.haskell.org/package/BiobaseNewick).+--+-- See nomenclature in 'ELynx.Tree'.+module ELynx.Tree.Export.Newick+  ( toNewick,+    toNewickBuilder,+  )+where++import qualified Data.ByteString.Builder as BB+import qualified Data.ByteString.Lazy.Char8 as BL+import Data.List (intersperse)+import ELynx.Tree.Named+import ELynx.Tree.Phylogeny+import ELynx.Tree.Rooted++-- | See 'toNewick'.+toNewickBuilder :: Named a => Tree Phylo a -> BB.Builder+toNewickBuilder t = go t <> BB.char8 ';'+  where+    go (Node b l []) = lbl b l+    go (Node b l ts) =+      BB.char8 '('+        <> 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)+    lbl x y =+      BB.lazyByteString (getName y)+        <> mBrLenBuilder x+        -- After reading several discussion, I go for the "more semantical+        -- form" with branch support values in square brackets.+        <> mBrSupBuilder x++-- | General conversion of a tree into a Newick 'BL.Bytestring'. Use provided+-- functions to extract node labels and branch lengths builder objects. See also+-- Biobase.Newick.Export.+toNewick :: Named a => Tree Phylo a -> BL.ByteString+toNewick = BB.toLazyByteString . toNewickBuilder

@@ -0,0 +1,31 @@+{-# LANGUAGE OverloadedStrings #-}++-- |+-- Module      :  ELynx.Tree.Export.Nexus+-- Description :  Export trees to Nexus files+-- Copyright   :  (c) Dominik Schrempf 2020+-- License     :  GPL-3+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Tue Apr 28 20:24:19 2020.+module ELynx.Tree.Export.Nexus+  ( toNexusTrees,+  )+where++import qualified Data.ByteString.Lazy.Char8 as BL+import ELynx.Export.Nexus+import ELynx.Tree.Export.Newick+import ELynx.Tree.Named+import ELynx.Tree.Phylogeny+import ELynx.Tree.Rooted++-- | Export a list of (NAME, TREE) to a Nexus file.+toNexusTrees :: Named a => [(BL.ByteString, Tree Phylo a)] -> BL.ByteString+toNexusTrees ts = toNexus "TREES" (map tree ts)++tree :: Named a => (BL.ByteString, Tree Phylo a) -> BL.ByteString+tree (n, t) = "  TREE " <> n <> " = " <> toNewick t

@@ -0,0 +1,235 @@+{-# LANGUAGE DeriveGeneric #-}++-- Module      :  ELynx.Tree.Import.Newick+-- Description :  Import Newick trees+-- Copyright   :  (c) Dominik Schrempf 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Thu Jan 17 14:56:27 2019.+--+-- Some functions are inspired by+-- [Biobase.Newick.Import](https://hackage.haskell.org/package/BiobaseNewick).+--+-- [Specifications](http://evolution.genetics.washington.edu/phylip/newicktree.html)+--+-- 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.+module ELynx.Tree.Import.Newick+  ( NewickFormat (..),+    description,+    newick,+    oneNewick,+    someNewick,+  )+where++import Control.Applicative+import Data.Aeson (FromJSON, ToJSON)+import Data.Attoparsec.ByteString.Char8+import qualified Data.ByteString.Char8 as BS+import ELynx.Tree.Measurable+import ELynx.Tree.Phylogeny+import ELynx.Tree.Rooted hiding (forest, label)+import ELynx.Tree.Supported+import GHC.Generics+import Prelude hiding (takeWhile)++-- | Newick tree format.+--+-- >>> unlines $ map (("- " <>) . description) (allValues :: [NewickFormat])+-- - Standard: Branch support values are stored in square brackets after branch lengths.+-- - IqTree:   Branch support values are stored as node names after the closing bracket of forests.+-- - 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.+data NewickFormat = Standard | IqTree | RevBayes+  deriving (Eq, Show, Read, Bounded, Enum, Generic)++instance FromJSON NewickFormat++instance ToJSON NewickFormat++-- | Short description of the supported Newick formats.+description :: NewickFormat -> String+description Standard =+  "Standard: Branch support values are stored in square brackets after branch lengths."+description IqTree =+  "IqTree:   Branch support values are stored as node names after the closing bracket of forests."+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 :: 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.+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.+someNewick :: NewickFormat -> Parser (Forest Phylo BS.ByteString)+someNewick Standard = someNewickStandard+someNewick IqTree = someNewickIqTree+someNewick RevBayes = someNewickRevBayes++-- 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"++-- Parse a single Newick tree. Fails when end of file is not reached.+oneNewickStandard :: Parser (Tree Phylo BS.ByteString)+oneNewickStandard = newickStandard <* endOfInput <?> "oneNewickStandard"++-- Parse one ore more Newick trees until end of file.+someNewickStandard :: Parser (Forest Phylo BS.ByteString)+someNewickStandard = some newickStandard <* endOfInput <?> "someNewickStandard"++tree :: Parser (Tree Phylo BS.ByteString)+tree = branched <|> leaf <?> "tree"++branched :: Parser (Tree Phylo BS.ByteString)+branched = (<?> "branched") $ do+  f <- forest+  n <- name+  p <- phylo+  return $ Node p n f++-- A 'forest' is a set of trees separated by @,@ and enclosed by parentheses.+forest :: Parser (Forest Phylo BS.ByteString)+forest = char '(' *> (tree `sepBy1` char ',') <* char ')' <?> "forest"++-- A 'leaf' has a 'name' and a 'phylo' branch.+leaf :: Parser (Tree Phylo BS.ByteString)+leaf = (<?> "leaf") $ do+  n <- name+  p <- phylo+  return $ Node p n []++nameChar :: Char -> Bool+nameChar c = c `notElem` " :;()[],"++-- A name can be any string of printable characters except blanks, colons,+-- semicolons, parentheses, and square brackets (and commas).+name :: Parser BS.ByteString+name = takeWhile nameChar <?> "name"++phylo :: Parser Phylo+phylo = Phylo <$> optional branchLength <*> optional branchSupport <?> "phylo"++-- Branch length.+branchLength :: Parser BranchLength+branchLength = char ':' *> double <?> "branchLength"++branchSupport :: Parser BranchSupport+branchSupport = (<?> "branchSupport") $+  do+    _ <- char '['+    s <- double+    _ <- char ']'+    return s++--------------------------------------------------------------------------------+-- IQ-TREE.++-- IQ-TREE stores the branch support as node names after the closing bracket of+-- a forest. Parse a single Newick tree. Also succeeds when more trees follow.+newickIqTree :: Parser (Tree Phylo BS.ByteString)+newickIqTree = skipWhile isSpace *> treeIqTree <* char ';' <* skipWhile isSpace <?> "newickIqTree"++-- See 'newickIqTree'. Parse a single Newick tree. Fails when end of file is not+-- reached.+oneNewickIqTree :: Parser (Tree Phylo BS.ByteString)+oneNewickIqTree = newickIqTree <* endOfInput <?> "oneNewickIqTree"++-- See 'newickIqTree'. Parse one ore more Newick trees until end of file.+someNewickIqTree :: Parser (Forest Phylo BS.ByteString)+someNewickIqTree = some newickIqTree <* endOfInput <?> "someNewickIqTree"++-- IQ-TREE stores the branch support as node names after the closing bracket of a forest.+treeIqTree :: Parser (Tree Phylo BS.ByteString)+treeIqTree = branchedIqTree <|> leaf <?> "treeIqTree"++-- IQ-TREE stores the branch support as node names after the closing bracket of a forest.+branchedIqTree :: Parser (Tree Phylo BS.ByteString)+branchedIqTree = (<?> "branchedIqTree") $ do+  f <- forestIqTree+  s <- optional double+  n <- name+  b <- optional branchLength+  return $ Node (Phylo b s) n f++-- IQ-TREE stores the branch support as node names after the closing bracket of a forest.+forestIqTree :: Parser (Forest Phylo BS.ByteString)+forestIqTree = (<?> "forestIqTree") $ do+  _ <- char '('+  f <- treeIqTree `sepBy1` char ','+  _ <- char ')'+  return f++--------------------------------------------------------------------------------+-- RevBayes.++-- RevBayes uses square brackets and key-value pairs to define information+-- about nodes and branches. Parse a single Newick tree. Also succeeds when more+-- trees follow.+--+-- XXX: Key value pairs are ignored at the moment.+newickRevBayes :: Parser (Tree Phylo BS.ByteString)+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.+oneNewickRevBayes :: Parser (Tree Phylo BS.ByteString)+oneNewickRevBayes = newickRevBayes <* endOfInput <?> "oneNewickRevBayes"++-- See 'newickRevBayes'. Parse one ore more Newick trees until end of file.+someNewickRevBayes :: Parser (Forest Phylo BS.ByteString)+someNewickRevBayes = some newickRevBayes <* endOfInput <?> "someNewickRevBayes"++treeRevBayes :: Parser (Tree Phylo BS.ByteString)+treeRevBayes = branchedRevBayes <|> leafRevBayes <?> "treeRevBayes"++branchedRevBayes :: Parser (Tree Phylo BS.ByteString)+branchedRevBayes = (<?> "branchedRevgBayes") $ do+  f <- forestRevBayes+  n <- nameRevBayes+  b <- optional branchLengthRevBayes+  return $ Node (Phylo b Nothing) n f++forestRevBayes :: Parser (Forest Phylo BS.ByteString)+forestRevBayes = (<?> "forestRevBayes") $ do+  _ <- char '('+  f <- treeRevBayes `sepBy1` char ','+  _ <- char ')'+  return f++nameRevBayes :: Parser BS.ByteString+nameRevBayes = name <* optional brackets <?> "nameRevBayes"++branchLengthRevBayes :: Parser BranchLength+branchLengthRevBayes = branchLength <* optional brackets <?> "branchLengthRevBayes"++leafRevBayes :: Parser (Tree Phylo BS.ByteString)+leafRevBayes = (<?> "leafRevBayes") $ do+  n <- nameRevBayes+  b <- optional branchLengthRevBayes+  return $ Node (Phylo b Nothing) n []++-- Drop anything between brackets.+brackets :: Parser ()+brackets = (<?> "brackets") $ do+  _ <- char '['+  _ <- takeWhile (/= ']')+  _ <- char ']'+  return ()

@@ -0,0 +1,45 @@+{-# LANGUAGE OverloadedStrings #-}++-- |+-- Module      :  ELynx.Tree.Import.Nexus+-- Description :  Import trees from Nexus files+-- Copyright   :  (c) Dominik Schrempf 2020+-- License     :  GPL-3+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Tue Apr 28 17:44:13 2020.+module ELynx.Tree.Import.Nexus+  ( nexusTrees,+  )+where++import Control.Applicative+import Data.Attoparsec.ByteString.Char8+import qualified Data.ByteString.Char8 as BS+import ELynx.Import.Nexus+import ELynx.Tree.Import.Newick+import ELynx.Tree.Phylogeny+import ELynx.Tree.Rooted+import Prelude hiding (takeWhile)++-- | Parse a Nexus files with a TREES block.+nexusTrees :: NewickFormat -> Parser [(BS.ByteString, Tree Phylo BS.ByteString)]+nexusTrees = nexus . trees++trees :: NewickFormat -> Block [(BS.ByteString, Tree Phylo BS.ByteString)]+trees f = Block "TREES" (some $ namedNewick f)++namedNewick :: NewickFormat -> Parser (BS.ByteString, Tree Phylo BS.ByteString)+namedNewick f = do+  _ <- skipWhile isSpace+  _ <- string "TREE"+  _ <- skipWhile isSpace+  n <- takeWhile1 (\x -> isAlpha_ascii x || isDigit x)+  _ <- skipWhile isSpace+  _ <- char '='+  _ <- skipWhile isSpace+  t <- newick f+  return (n, t)

@@ -0,0 +1,124 @@+-- |+-- Module      :  ELynx.Tree.Measurable+-- Description :  Measurable branch labels+-- Copyright   :  (c) Dominik Schrempf 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Thu Jan 17 14:16:34 2019.+--+-- Non-negativity of branch lengths is not (yet) ensured. To ensure+-- non-negativity, a newtype wrapper could be used, but this would be a major+-- refactor.+module ELynx.Tree.Measurable+  ( BranchLength,+    Measurable (..),+    applyStem,+    getStem,+    setStem,+    height,+    rootHeight,+    distancesOriginLeaves,+    totalBranchLength,+    normalizeBranchLengths,+    normalizeHeight,+    ultrametric,+    makeUltrametric,+  )+where++import Data.Bifoldable+import Data.Bifunctor+import ELynx.Tree.Rooted++-- | Branch length.+type BranchLength = Double++-- | A branch label with measurable and modifiable branch length.+class Measurable e where+  -- | Length of attached branch.+  getLen :: e -> BranchLength++  -- | Set attached branch length.+  setLen :: BranchLength -> e -> e++instance Measurable Double where+  getLen = id+  setLen = const++-- Apply a function to a branch support label.+apply :: Measurable e => (BranchLength -> BranchLength) -> e -> e+apply f l = setLen (f s) l where s = getLen l++-- | Lengthen the stem of a tree.+applyStem :: Measurable e => (BranchLength -> BranchLength) -> Tree e a -> Tree e a+applyStem f t = t {branch = apply f b}+  where+    b = branch t++-- | Get the length of the stem of a tree.+getStem :: Measurable e => Tree e a -> BranchLength+getStem (Node br _ _) = getLen br++-- | Set the length of the stem of a tree.+setStem :: Measurable e => BranchLength -> Tree e a -> Tree e a+setStem x = applyStem (const x)++-- | The maximum distance between origin and leaves.+--+-- The height includes the length of the stem.+height :: Measurable e => Tree e a -> BranchLength+height = maximum . distancesOriginLeaves++-- | The maximum distance between root node and leaves.+rootHeight :: Measurable e => Tree e a -> BranchLength+rootHeight (Node _ _ []) = 0+rootHeight t = maximum $ concatMap distancesOriginLeaves (forest t)++-- | Distances from the origin of a tree to the leaves.+--+-- The distances include the length of the stem.+distancesOriginLeaves :: Measurable e => Tree e a -> [BranchLength]+distancesOriginLeaves (Node br _ []) = [getLen br]+distancesOriginLeaves (Node br _ ts) = map (getLen br +) (concatMap distancesOriginLeaves ts)++-- | Total branch length of a tree.+totalBranchLength :: Measurable e => Tree e a -> BranchLength+totalBranchLength = bifoldl' (+) const 0 . first getLen++-- | Normalize branch lengths so that the sum is 1.0.+normalizeBranchLengths :: Measurable e => Tree e a -> Tree e a+normalizeBranchLengths t = first (apply (/ s)) t+  where+    s = totalBranchLength t++-- | Normalize height of tree to 1.0.+normalizeHeight :: Measurable e => Tree e a -> Tree e a+normalizeHeight t = first (apply (/ h)) t+  where+    h = height t++eps :: Double+eps = 1e-12++allNearlyEqual :: [Double] -> Bool+allNearlyEqual [] = True+allNearlyEqual xs = all (\y -> eps > abs (x - y)) (tail xs)+  where+    x = head xs++-- | Check if a tree is ultrametric.+ultrametric :: Measurable e => Tree e a -> Bool+ultrametric = allNearlyEqual . distancesOriginLeaves++-- | Elongate terminal branches such that the tree becomes ultrametric.+makeUltrametric :: Measurable e => Tree e a -> Tree e a+makeUltrametric t = go 0 t+  where+    h = height t+    go :: Measurable e => BranchLength -> Tree e a -> Tree e a+    go h' (Node br lb []) = let dh = h - h' - getLen br in Node (apply (+ dh) br) lb []+    go h' (Node br lb ts) = let h'' = h' + getLen br in Node br lb $ map (go h'') ts

@@ -0,0 +1,43 @@+-- |+-- Module      :  ELynx.Tree.Named+-- Description :  Trees with named nodes+-- Copyright   :  (c) Dominik Schrempf 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Thu Jan 24 20:09:20 2019.+module ELynx.Tree.Named+  ( Named (..),+  )+where++import qualified Data.ByteString.Builder as BB+import qualified Data.ByteString.Char8 as BS+import qualified Data.ByteString.Lazy.Char8 as BL+import Data.Double.Conversion.ByteString as BC++-- | Data types with names.+class Named a where+  -- Use lazy byte strings because Newick strings are built using chunks.+  getName :: a -> BL.ByteString++instance Named () where+  getName = const BL.empty++instance Named Int where+  getName = BB.toLazyByteString . BB.intDec++instance Named Double where+  getName = BL.fromStrict . toShortest++instance Named Char where+  getName = BB.toLazyByteString . BB.char8++instance Named BL.ByteString where+  getName = id++instance Named BS.ByteString where+  getName = BL.fromStrict

@@ -0,0 +1,137 @@+-- |+-- Module      :  ELynx.Tree.Partition+-- Description :  Partitions on rose trees+-- Copyright   :  (c) Dominik Schrempf 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Thu Dec 12 12:58:49 2019.+--+-- A multifurcation induces a 'Partition', similar to branches inducing+-- 'ELynx.Tree.Bipartition's.+module ELynx.Tree.Partition+  ( -- * Data type+    Partition (fromPartition),+    mp,+    mpUnsafe,+    bpToMp,+    mpHuman,++    -- * Work with 'Partition's+    partition,+    partitions,+    compatible,+  )+where++import Data.List hiding (partition)+import Data.Set (Set)+import qualified Data.Set as S+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'.+--+-- 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.+newtype Partition a = Partition+  { fromPartition :: Set (Set a)+  }+  deriving (Eq, Ord, Show, Read)++-- 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+  [] -> Left "mp: Empty list."+  xs' -> Right $ mpUnsafe xs'++-- | Create a partition.+mpUnsafe :: Ord a => [Set a] -> Partition a+mpUnsafe xs = Partition (S.fromList xs)++-- | Convert a bipartition to a partition.+bpToMp :: Ord a => Bipartition a -> Partition a+bpToMp = mpUnsafe . tupleToList . fromBipartition+  where+    -- Be careful with tuples, because 'toList' does something very weird. It only+    -- takes the second element of the tuple!+    --+    -- toList :: Foldable t => t a -> [a]+    tupleToList (x, y) = [x, y]++-- | 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) =+  "(" ++ intercalate "|" (map setShow (S.toList xs)) ++ ")"++-- Show the elements of a set in a human readable format.+setShow :: Show a => Set a -> String+setShow = intercalate "," . map show . S.toList++-- | Get partition defined by the root of the tree.+--+-- Return 'Left' if:+-- - the tree is a leaf;+-- - the tree contains duplicate leaves.+partition :: Ord a => Tree e a -> Either String (Partition a)+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++-- | Get all 'Partition's of a tree.+--+-- Return 'Left' if tree contains duplicate leaves.+partitions :: Ord a => Tree e a -> Either String (Set (Partition a))+partitions t+  | duplicateLeaves t = Left "partitions: Tree contains duplicate leaves."+  | otherwise = Right $ partitions' S.empty $ S.fromList <$> groups t++-- See 'partitions', but do not check if leaves are unique.+partitions' :: Ord a => Set a -> Tree e (Set a) -> Set (Partition a)+partitions' _ (Node _ _ []) = S.empty+partitions' p t@(Node _ _ ts) =+  S.unions $+    either (const S.empty) S.singleton (mp (p : map label ts)) :+    zipWith partitions' cs ts+  where+    cs = getComplementaryLeaves p t++-- | 'Partition's are compatible if they do not contain conflicting+-- information. This function checks if two partitions are compatible with+-- each other. Thereby, a variation of the following algorithm is used:+--+-- @+-- mp1 `compatible` mp2+-- for set1 in mp1:+--   for set2 in mp2:+--     if set1 `S.isSubSetOf` set2:+--       remove set1 from mp1+--     if set2 `S.isSubSetOf` set1:+--       remove set2 from mp2+-- if either mp2 or mp2 is empty, they are compatible+-- @+compatible :: (Show a, Ord a) => Partition a -> Partition a -> Bool+compatible l r = S.null (S.filter (`remove` rs) ls) || S.null (S.filter (`remove` ls) rs)+  where+    ls = fromPartition l+    rs = fromPartition r++remove :: Ord a => Set a -> Set (Set a) -> Bool+remove s = not . any (s `S.isSubsetOf`)

@@ -0,0 +1,487 @@+{-# LANGUAGE DeriveAnyClass #-}+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE DerivingVia #-}++-- |+-- Module      :  ELynx.Tree.Phylogeny+-- Description :  Phylogenetic trees+-- Copyright   :  (c) Dominik Schrempf 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Thu Jan 17 16:08:54 2019.+--+-- A phylogeny is a 'Tree' with unique leaf labels, and the order of the trees+-- in the sub-forest is considered to be meaningless.+--+-- Internally, however, the underlying 'Tree' data structure stores the+-- sub-forest as a list, which has a specific order. Hence, we have to do some+-- tricks when comparing trees, and tree comparison is slow.+--+-- Also, the uniqueness of the leaves is not ensured by the data type, but has+-- to be checked at runtime. Functions relying on the tree to have unique leaves+-- do perform this check, and return 'Left' with an error message, if the tree+-- has duplicate leaves.+--+-- Note: 'Tree's are rooted.+--+-- Note: 'Tree's encoded in Newick format correspond to rooted trees. By+-- convention only, a tree parsed from Newick format is usually thought to be+-- unrooted, when the root node is multifurcating and has three children. This+-- convention is not used here. Newick trees are just parsed as they are, and a+-- rooted tree is returned.+--+-- The bifurcating root of a tree can be changed with 'rootAt' or 'midpoint'; a+-- list of all rooted trees is returned by 'roots'.+--+-- Trees with multifurcating root nodes can be rooted using 'outgroup'.+module ELynx.Tree.Phylogeny+  ( -- * Functions+    equal,+    intersect,+    bifurcating,+    outgroup,+    midpoint,+    roots,+    rootAt,++    -- * Branch labels+    Phylo (..),+    toPhyloTree,+    measurableToPhyloTree,+    supportedToPhyloTree,+    Length (..),+    phyloToLengthTree,+    Support (..),+    phyloToSupportTree,+    phyloToSupportTreeUnsafe,+    PhyloExplicit (..),+    toExplicitTree,+  )+where++import Control.DeepSeq+import Data.Aeson+import Data.Bifoldable+import Data.Bifunctor+import Data.Bitraversable+import Data.List hiding (intersect)+import Data.Maybe+import Data.Monoid+import Data.Semigroup+import Data.Set (Set)+import qualified Data.Set as S+import ELynx.Tree.Bipartition+import ELynx.Tree.Measurable+import ELynx.Tree.Rooted+import ELynx.Tree.Splittable+import ELynx.Tree.Supported+import GHC.Generics++-- | The equality check is slow because the order of children is considered to+-- be arbitrary.+equal :: (Eq e, Eq a) => Tree e a -> Tree e a -> Bool+equal ~(Node brL lbL tsL) ~(Node brR lbR tsR) =+  (brL == brR)+    && (lbL == lbR)+    && (length tsL == length tsR)+    && all (`elem` tsR) tsL++-- | Compute the intersection of trees.+--+-- The intersections are the largest subtrees sharing the same leaf set.+--+-- Degree two nodes are pruned with 'prune'.+--+-- Return 'Left' if:+-- - the intersection of leaves is empty.+intersect ::+  (Semigroup e, Eq e, Ord a) => Forest e a -> Either String (Forest e a)+intersect ts+  | S.null lvsCommon = Left "intersect: Intersection of leaves is empty."+  | otherwise = case sequence [dropLeavesWith (predicate ls) t | (ls, t) <- zip leavesToDrop ts] of+    Nothing -> Left "intersect: A tree is empty."+    Just ts' -> Right ts'+  where+    -- Leaf sets.+    lvss = map (S.fromList . leaves) ts+    -- Common leaf set.+    lvsCommon = foldl1' S.intersection lvss+    -- Leaves to drop for each tree in the forest.+    leavesToDrop = map (S.\\ lvsCommon) lvss+    -- Predicate.+    predicate lvsToDr l = l `S.member` lvsToDr++-- | Check if a tree is bifurcating.+--+-- A Bifurcating tree only contains degree one (leaves) and degree three nodes+-- (internal bifurcating nodes).+bifurcating :: Tree e a -> Bool+bifurcating (Node _ _ []) = True+bifurcating (Node _ _ [x, y]) = bifurcating x && bifurcating y+bifurcating _ = False++-- I believe that manual treatment with 'outgroup' is preferable.++-- -- | Remove multifurcations.+-- --+-- -- A caterpillar like bifurcating structure is used to resolve all+-- -- multifurcations on a tree.+-- --+-- -- Multifurcating nodes are copied and branches are 'split'.+-- resolve :: Splittable e => Tree e a -> Tree e a+-- resolve t@(Node _ _ []) = t+-- resolve (Node br lb [x]) = Node br lb [resolve x]+-- resolve (Node br lb [x, y]) = Node br lb $ map resolve [x, y]+-- resolve (Node br lb (Node brL lbL xsL : xs)) = Node br lb [Node brL' lbL (map resolve xsL), Node brL' lb (map resolve xs)]+--   where brL' = split brL++-- | Resolve a multifurcating root using an outgroup.+--+-- A bifurcating root node with the provided label is introduced. The affected+-- branch is 'split'.+--+-- Note, the degree of the former root node is decreased by one.+--+-- If the root node is bifurcating, use 'rootAt'.+--+-- Return 'Left' if+-- - the tree has duplicate leaves;+-- - the root node is not multifurcating;+-- - the provided outgroup is not found on the tree or is polyphyletic.+outgroup :: (Semigroup e, Splittable e, Ord a) => Set a -> a -> Tree e a -> Either String (Tree e a)+outgroup _ _ (Node _ _ []) = Left "outgroup: Root node is a leaf."+outgroup _ _ (Node _ _ [_]) = Left "outgroup: Root node has degree two."+outgroup _ _ (Node _ _ [_, _]) = Left "outgroup: Root node is bifurcating."+outgroup o r t@(Node b l ts)+  | duplicateLeaves t = Left "outgroup: Tree has duplicate leaves."+  | otherwise = do+    bip <- bp o (S.fromList lvs S.\\ o)+    rootAt bip t'+  where+    lvs = leaves t+    (Node brO lbO tsO) = head ts+    -- Introduce a bifurcating root node.+    t' = Node b r [Node (split brO) lbO tsO, Node (split brO) l (tail ts)]++-- The 'midpoint' algorithm is pretty stupid because it calculates all rooted+-- trees and then finds the one minimizing the difference between the heights of+-- the left and right sub tree. Actually, one just needs to move left or right,+-- with the aim to minimize the height difference between the left and right sub+-- tree.++-- | Root tree at the midpoint.+--+-- Return 'Left' if+-- - the root node is not bifurcating.+midpoint :: (Semigroup e, Splittable e, Measurable e) => Tree e a -> Either String (Tree e a)+midpoint (Node _ _ []) = Left "midpoint: Root node is a leaf."+midpoint (Node _ _ [_]) = Left "midpoint: Root node has degree two."+midpoint t@(Node _ _ [_, _]) = getMidpoint <$> roots t+midpoint _ = Left "midpoint: Root node is multifurcating."++findMinIndex :: Ord a => [a] -> Int+findMinIndex (x : xs) = go (0, x) 1 xs+  where+    go (i, _) _ [] = i+    go (i, z) j (y : ys) = if z < y then go (i, z) (j + 1) ys else go (j, y) (j + 1) ys+findMinIndex [] = error "findMinIndex: Empty list."++getMidpoint :: Measurable e => [Tree e a] -> Tree e a+getMidpoint ts = case t of+  (Node br lb [l, r]) ->+    let hl = height l+        hr = height r+        dh = (hl - hr) / 2+     in Node br lb [applyStem (subtract dh) l, applyStem (+ dh) r]+  -- Explicitly use 'error' here, because roots is supposed to return trees with+  -- bifurcating root nodes.+  _ -> error "getMidpoint: Root node is not bifurcating."+  where+    dhs = map getDeltaHeight ts+    i = findMinIndex dhs+    t = ts !! i++-- find index of minimum; take this tree and move root to the midpoint of the branch++-- Get delta height of left and right sub tree.+getDeltaHeight :: Measurable e => Tree e a -> Double+getDeltaHeight (Node _ _ [l, r]) = abs $ height l - height r+-- Explicitly use 'error' here, because roots is supposed to return trees with+-- bifurcating root nodes.+getDeltaHeight _ = error "getDeltaHeight: Root node is not bifurcating."++-- | For a rooted tree with a bifurcating root node, get all possible rooted+-- trees.+--+-- The root node is moved.+--+-- For a tree with @l=2@ leaves, there is one rooted tree. For a bifurcating+-- tree with @l>2@ leaves, there are @(2l-3)@ rooted trees. For a general tree+-- with a bifurcating root node, and a total number of @n>2@ nodes, there are+-- (n-2) rooted trees.+--+-- Moving a multifurcating root node to another branch would change the degree+-- of the root node. Hence, a bifurcating root is required. To resolve a+-- multifurcating root, please use 'outgroup'.+--+-- Connect branches according to the provided 'Semigroup' instance.+--+-- Upon insertion of the root, split the affected branch into one out of two+-- equal entities according to a given function.+--+-- Return 'Left' if the root node is not 'bifurcating'.+roots :: (Semigroup e, Splittable e) => Tree e a -> Either String (Forest e a)+roots (Node _ _ []) = Left "roots: Root node is a leaf."+roots (Node _ _ [_]) = Left "roots: Root node has degree two."+roots t@(Node b c [tL, tR]) = Right $ t : descend b c tR tL ++ descend b c tL tR+roots _ = Left "roots: Root node is multifurcating."++complementaryForests :: Tree e a -> Forest e a -> [Forest e a]+complementaryForests t ts = [t : take i ts ++ drop (i + 1) ts | i <- [0 .. (n -1)]]+  where+    n = length ts++-- From the bifurcating root, descend into one of the two pits.+--+-- descend splitFunction rootBranch rootLabel complementaryTree downwardsTree+descend :: (Semigroup e, Splittable e) => e -> a -> Tree e a -> Tree e a -> Forest e a+descend _ _ _ (Node _ _ []) = []+descend brR lbR tC (Node brD lbD tsD) =+  [ Node brR lbR [Node (split brDd) lbD f, Node (split brDd) lbDd tsDd]+    | (Node brDd lbDd tsDd, f) <- zip tsD cfs+  ]+    ++ concat+      [ descend brR lbR (Node (split brDd) lbD f) (Node (split brDd) lbDd tsDd)+        | (Node brDd lbDd tsDd, f) <- zip tsD cfs+      ]+  where+    brC' = branch tC <> brD+    tC' = tC {branch = brC'}+    cfs = complementaryForests tC' tsD++-- | Root a tree at a specific position.+--+-- Root the tree at the branch defined by the given bipartition. The original+-- root node is moved to the new position.+--+-- The root node must be bifurcating (see 'roots' and 'outgroup').+--+-- Connect branches according to the provided 'Semigroup' instance.+--+-- Upon insertion of the root, split the affected branch according to the+-- provided 'Splittable' instance.+--+-- Return 'Left', if:+-- - the root node is not bifurcating;+-- - the tree has duplicate leaves;+-- - the bipartition does not match the leaves of the tree.+rootAt ::+  (Semigroup e, Splittable e, Eq a, Ord a) =>+  Bipartition a ->+  Tree e a ->+  Either String (Tree e a)+rootAt b t+  -- Tree is checked for being bifurcating in 'roots'.+  --+  -- Do not use 'duplicateLeaves' here, because we also need to compare the leaf+  -- set with the bipartition.+  | length lvLst /= S.size lvSet = Left "rootAt: Tree has duplicate leaves."+  | toSet b /= lvSet = Left "rootAt: Bipartition does not match leaves of tree."+  | otherwise = rootAt' b t+  where+    lvLst = leaves t+    lvSet = S.fromList $ leaves t++-- Assume the leaves of the tree are unique.+rootAt' ::+  (Semigroup e, Splittable e, Ord a) =>+  Bipartition a ->+  Tree e a ->+  Either String (Tree e a)+rootAt' b t = do+  ts <- roots t+  case find (\x -> Right b == bipartition x) ts of+    Nothing -> Left "rootAt': Bipartition not found on tree."+    Just t' -> Right t'++-- | Branch label for phylogenetic trees.+--+-- Branches may have a length and a support value.+data Phylo = Phylo+  { brLen :: Maybe BranchLength,+    brSup :: Maybe BranchSupport+  }+  deriving (Read, Show, Eq, Ord, Generic, NFData)++instance Semigroup Phylo where+  Phylo mBL mSL <> Phylo mBR mSR =+    Phylo+      (getSum <$> (Sum <$> mBL) <> (Sum <$> mBR))+      (getMin <$> (Min <$> mSL) <> (Min <$> mSR))++instance ToJSON Phylo++instance FromJSON Phylo++-- | Set all branch length and support values to 'Just' the value.+--+-- Useful to export a tree with branch lengths in Newick format.+toPhyloTree :: (Measurable e, Supported e) => Tree e a -> Tree Phylo a+toPhyloTree = first toPhyloLabel++toPhyloLabel :: (Measurable e, Supported e) => e -> Phylo+toPhyloLabel x = Phylo (Just $ getLen x) (Just $ getSup x)++-- | Set all branch support values to 'Nothing'.+--+-- Useful to export a tree with branch lengths to Newick format.+measurableToPhyloTree :: Measurable e => Tree e a -> Tree Phylo a+measurableToPhyloTree = first measurableToPhyloLabel++measurableToPhyloLabel :: Measurable e => e -> Phylo+measurableToPhyloLabel x = Phylo (Just $ getLen x) Nothing++-- | Set all branch lengths to 'Nothing'.+--+-- Useful to export a tree with branch support to Newick format.+supportedToPhyloTree :: Supported e => Tree e a -> Tree Phylo a+supportedToPhyloTree = first supportedToPhyloLabel++supportedToPhyloLabel :: Supported e => e -> Phylo+supportedToPhyloLabel x = Phylo Nothing (Just $ getSup x)++-- | Branch length label.+--+-- For conversion, see 'phyloToLengthTree'.+newtype Length = Length {fromLength :: BranchLength}+  deriving (Read, Show, Eq, Ord, Generic, NFData)+  deriving (Num, Fractional, Floating) via Double+  deriving (Semigroup, Monoid) via Sum Double++instance Measurable Length where+  getLen = fromLength+  setLen b _ = Length b++instance Splittable Length where+  split = Length . (/ 2.0) . fromLength++instance ToJSON Length++instance FromJSON Length++-- | If root branch length is not available, set it to 0.+--+-- Return 'Left' if any other branch length is unavailable.+phyloToLengthTree :: Tree Phylo a -> Either String (Tree Length a)+phyloToLengthTree =+  maybe (Left "phyloToLengthTree: Length unavailable for some branches.") Right+    . bitraverse toLength pure+    . cleanRootLength++cleanRootLength :: Tree Phylo a -> Tree Phylo a+cleanRootLength (Node (Phylo Nothing s) l f) = Node (Phylo (Just 0) s) l f+cleanRootLength t = t++toLength :: Phylo -> Maybe Length+toLength p = Length <$> brLen p++-- | Branch support label.+--+-- For conversion, see 'phyloToSupportTree'.+newtype Support = Support {fromSupport :: BranchSupport}+  deriving (Read, Show, Eq, Ord, Generic, NFData)+  deriving (Num, Fractional, Floating) via Double+  deriving (Semigroup) via Min Double++instance Supported Support where+  getSup = fromSupport+  setSup s _ = Support s++instance Splittable Support where+  split = id++instance ToJSON Support++instance FromJSON Support++-- | Set branch support values of branches leading to the leaves and of the root+-- branch to maximum support.+--+-- Return 'Left' if any other branch has no available support value.+phyloToSupportTree :: Tree Phylo a -> Either String (Tree Support a)+phyloToSupportTree t =+  maybe+    (Left "phyloToSupportTree: Support unavailable for some branches.")+    Right+    $ bitraverse toSupport pure $+      cleanLeafSupport m $+        cleanRootSupport m t+  where+    m = getMaxSupport t++-- | Set all unavailable branch support values to maximum support.+phyloToSupportTreeUnsafe :: Tree Phylo a -> Tree Support a+phyloToSupportTreeUnsafe t = cleanSupport m t+  where+    m = getMaxSupport t++-- If all branch support values are below 1.0, set the max support to 1.0.+getMaxSupport :: Tree Phylo a -> BranchSupport+getMaxSupport = fromJust . max (Just 1.0) . bimaximum . bimap brSup (const Nothing)++cleanRootSupport :: BranchSupport -> Tree Phylo a -> Tree Phylo a+cleanRootSupport maxSup (Node (Phylo b Nothing) l xs) = Node (Phylo b (Just maxSup)) l xs+cleanRootSupport _ t = t++cleanLeafSupport :: BranchSupport -> Tree Phylo a -> Tree Phylo a+cleanLeafSupport s (Node (Phylo b Nothing) l []) = Node (Phylo b (Just s)) l []+cleanLeafSupport s (Node b l xs) = Node b l $ map (cleanLeafSupport s) xs++toSupport :: Phylo -> Maybe Support+toSupport (Phylo _ Nothing) = Nothing+toSupport (Phylo _ (Just s)) = Just $ Support s++cleanSupport :: BranchSupport -> Tree Phylo a -> Tree Support a+cleanSupport maxSup (Node (Phylo _ s) l xs) = Node (Support $ fromMaybe maxSup s) l $ map (cleanSupport maxSup) xs++-- | Explicit branch label for phylogenetic trees.+data PhyloExplicit = PhyloExplicit+  { sBrLen :: BranchLength,+    sBrSup :: BranchSupport+  }+  deriving (Read, Show, Eq, Ord, Generic)++instance Semigroup PhyloExplicit where+  PhyloExplicit bL sL <> PhyloExplicit bR sR = PhyloExplicit (bL + bR) (min sL sR)++instance Measurable PhyloExplicit where+  getLen = sBrLen+  setLen b l = l {sBrLen = b}++instance Splittable PhyloExplicit where+  split l = l {sBrLen = b'}+    where+      b' = sBrLen l / 2.0++instance Supported PhyloExplicit where+  getSup = sBrSup+  setSup s l = l {sBrSup = s}++instance ToJSON PhyloExplicit++instance FromJSON PhyloExplicit++-- | Conversion to a 'PhyloExplicit' tree.+--+-- See 'phyloToLengthTree' and 'phyloToSupportTree'.+toExplicitTree :: Tree Phylo a -> Either String (Tree PhyloExplicit a)+toExplicitTree t = do+  lt <- first fromLength <$> phyloToLengthTree t+  st <- first fromSupport <$> phyloToSupportTree t+  case zipTreesWith PhyloExplicit const lt st of+    Nothing -> error "toExplicitTree: This is a bug. Can not zip two trees with the same topology."+    Just zt -> return zt

@@ -0,0 +1,355 @@+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE DeriveGeneric #-}++-- |+-- Module      :  ELynx.Tree.Rooted+-- Description :  Rooted trees with labeled branches+-- Copyright   :  (c) Dominik Schrempf 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Thu Jan 17 09:57:29 2019.+--+-- Rooted 'Tree's differes from a classical rose 'Data.Tree.Tree' in that it has+-- labeled branches.+--+-- For rooted topologies, please see 'ELynx.Topology.Rooted'.+--+-- A 'Tree' is defined as:+--+-- @+-- data Tree e a = Node+--   { branch :: e,+--     label :: a,+--     forest :: Forest e a+--   }+-- @+--+-- where+--+-- @+-- type Forest e a = [Tree e a]+-- @+--+-- This means, that the word 'Node' is reserved for the constructor of a tree,+-- and that a 'Node' has an attached 'branch', a 'label', and a sub-'forest'.+-- The value constructor /Node/ and the record function /label/ are not to be+-- confused. The elements of the sub-forest are often called /children/.+--+-- With respect to phylogenetic analyses, using the 'Tree' data type has some+-- disadvantages:+--+-- 1. All trees are rooted. Unrooted trees can be treated with a rooted data+-- structure, as it is used here. However, some functions may be meaningless.+--+-- 2. Changing branch labels, node labels, or the topology of the tree are slow+-- operations, especially, when the changes are close to the leaves of the tree.+--+-- In mathematical terms: A 'Tree' is a directed acyclic graph without loops,+-- with vertex labels, with edge labels. Let me know if this definition is+-- incomplete.+module ELynx.Tree.Rooted+  ( -- * Data type+    Tree (..),+    Forest,+    toTreeBranchLabels,+    toTreeNodeLabels,++    -- * Access leaves, branches and labels+    leaves,+    duplicateLeaves,+    branches,+    setBranches,+    labels,+    setLabels,+    identify,++    -- * Change structure+    degree,+    prune,+    dropNodesWith,+    dropLeavesWith,+    zipTreesWith,+    zipTrees,+  )+where++import Control.Applicative+import Control.Comonad+import Control.DeepSeq+import Control.Monad+import Control.Monad.Fix+import Data.Aeson+import Data.Bifoldable+import Data.Bifunctor+import Data.Bitraversable+import Data.Data+import Data.Foldable+import Data.List+import Data.Maybe+import qualified Data.Set as S+import qualified Data.Tree as T+import GHC.Generics++-- | 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.+--+-- Lifted instances are not provided.+data Tree e a = Node+  { branch :: e,+    label :: a,+    forest :: Forest e a+  }+  deriving (Eq, Read, Show, Data, Generic)++-- | A shorthand.+type Forest e a = [Tree e a]++-- | Map over node labels.+instance Functor (Tree e) where+  fmap f ~(Node br lb ts) = Node br (f lb) $ map (fmap f) ts+  x <$ ~(Node br _ ts) = Node br x (map (x <$) ts)++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+  second g ~(Node br lb ts) = Node br (g lb) $ map (second g) ts++-- | Combine node labels in pre-order.+instance Foldable (Tree e) where+  foldMap f ~(Node _ lb ts) = f lb <> foldMap (foldMap f) ts+  null _ = False+  {-# INLINE null #-}+  toList = labels+  {-# INLINE toList #-}++instance Bifoldable Tree where+  bifoldMap f g ~(Node br lb ts) = f br <> g lb <> foldMap (bifoldMap f g) ts++instance Traversable (Tree e) where+  traverse g ~(Node br lb ts) = Node br <$> g lb <*> traverse (traverse g) ts++instance Bitraversable Tree where+  bitraverse f g ~(Node br lb ts) = Node <$> f br <*> g lb <*> traverse (bitraverse f g) ts++-- The following code provides a zip-like applicative instance. However,+-- the zip-like instance makes the Monad instance meaningless. So, either we+-- provide only 'Applicative' in zip-like form, or we use the classic instance+-- for 'Applicative' and 'Monad'.++-- -- | Note: The 'Applicative' instance of 'Tree' is similar to the one of+-- -- 'Control.Applicative.ZipList', and differs from the instance of+-- -- 'Data.Tree.Tree'!+-- --+-- -- >>> let t = Node "" 0 [Node "" 1 [], Node "" 2 []] :: Tree String Int+-- -- >>> let f = Node "+3" (+3) [Node "*5" (*5) [], Node "+10" (+10) []] :: Tree String (Int -> Int)+-- -- >>> f <*> t+-- -- Node {branch = "+3", label = 3, forest = [Node {branch = "*5", label = 5, forest = []},Node {branch = "+10", label = 12, forest = []}]}+-- --+-- -- Note: The 'Monoid' instance of the branch labels determines how the branches+-- -- are combined. For example, distances can be summed using the+-- -- 'Data.Monoid.Sum' monoid.+-- instance Monoid e => Applicative (Tree e) where+--   pure lb = Node mempty lb []+--   ~(Node brF lbF tsF) <*> ~(Node brX lbX tsX) =+--     Node (brF <> brX) (lbF lbX) (zipWith (<*>) tsF tsX)+--   liftA2 f ~(Node brX lbX tsX) ~(Node brY lbY tsY) =+--     Node (brX <> brY) (f lbX lbY) (zipWith (liftA2 f) tsX tsY)+--   ~(Node brX _ tsX) *> ~(Node brY lbY tsY) =+--     Node (brX <> brY) lbY (zipWith (*>) tsX tsY)+--   ~(Node brX lbX tsX) <* ~(Node brY _ tsY) =+--     Node (brX <> brY) lbX (zipWith (<*) tsX tsY)++-- | The 'Semigroup' instance of the branch labels determines how the+-- branches are combined. For example, distances can be summed using+-- 'Data.Semigroup.Sum'.+--+-- The 'Monoid' instance of the branch labels determines the default branch+-- label when using 'pure'.+instance Monoid e => Applicative (Tree e) where+  pure lb = Node mempty lb []+  ~(Node brF lbF tsF) <*> ~tx@(Node brX lbX tsX) =+    Node (brF <> brX) (lbF lbX) (map (lbF <$>) tsX ++ map (<*> tx) tsF)+  liftA2 f ~(Node brX lbX tsX) ~ty@(Node brY lbY tsY) =+    Node (brX <> brY) (f lbX lbY) (map (f lbX <$>) tsY ++ map (\tx -> liftA2 f tx ty) tsX)+  ~(Node brX _ tsX) *> ~ty@(Node brY lbY tsY) =+    Node (brX <> brY) lbY (tsY ++ map (*> ty) tsX)+  ~(Node brX lbX tsX) <* ~ty@(Node brY _ tsY) =+    Node (brX <> brY) lbX (map (lbX <$) tsY ++ map (<* ty) tsX)++-- | The 'Semigroup' instance of the branch labels determines how the branches+-- are combined. For example, distances can be summed using+-- 'Data.Semigroup.Sum'.+instance Monoid e => Monad (Tree e) where+  ~(Node br lb ts) >>= f = case f lb of+    Node br' lb' ts' -> Node (br <> br') lb' (ts' ++ map (>>= f) ts)++-- -- Cannot provide MonadZip instance because branch labels cannot be+-- -- recovered from combined label.+-- instance Monoid e => MonadZip (Tree e) where+--   mzipWith f (Node brL lbL tsL) (Node brR lbR tsR) =+--     Node (brL <> brR) (f lbL lbR) (mzipWith (mzipWith f) tsL tsR)+--+--   munzip (Node br (lbL, lbR) ts) = (Node ? lbL tsL, Node ? lbR tsR)+--     where+--       (tsL, tsR) = munzip (map munzip ts)++instance Monoid e => MonadFix (Tree e) where+  mfix = mfixTree++mfixTree :: (a -> Tree e a) -> Tree e a+mfixTree f+  | Node br lb ts <- fix (f . label) =+    Node+      br+      lb+      ( zipWith+          (\i _ -> mfixTree ((!! i) . forest . f))+          [0 ..]+          ts+      )++instance Comonad (Tree e) where+  duplicate t@(Node br _ ts) = Node br t (map duplicate ts)+  extract (Node _ lb _) = lb+  {-# INLINE extract #-}++instance (NFData e, NFData a) => NFData (Tree e a) where+  rnf (Node br lb ts) = rnf br `seq` rnf lb `seq` rnf ts++instance (ToJSON e, ToJSON a) => ToJSON (Tree e a)++instance (FromJSON e, FromJSON a) => FromJSON (Tree e a)++-- | Conversion to 'T.Tree' using branch labels.+toTreeBranchLabels :: Tree e a -> T.Tree e+toTreeBranchLabels (Node br _ ts) = T.Node br (map toTreeBranchLabels ts)++-- | Conversion to 'T.Tree' using node labels.+toTreeNodeLabels :: Tree e a -> T.Tree a+toTreeNodeLabels (Node _ lb ts) = T.Node lb (map toTreeNodeLabels ts)++-- | Get leaves.+leaves :: Tree e a -> [a]+leaves (Node _ lb []) = [lb]+leaves (Node _ _ ts) = concatMap leaves ts++-- | Check if a tree has duplicate leaves.+duplicateLeaves :: Ord a => Tree e a -> Bool+duplicateLeaves = duplicates . leaves++-- | Get branch labels in pre-order.+branches :: Tree e a -> [e]+branches t = squish t []+  where+    squish (Node br _ ts) xs = br : foldr squish xs ts++-- | Set branch labels in pre-order.+--+-- Return 'Nothing' if the provided list of branch labels is too short.+setBranches :: Bitraversable t => [f] -> t e a -> Maybe (t f a)+setBranches xs = bisequenceA . snd . bimapAccumL setBranch noChange xs+  where+    setBranch [] _ = ([], Nothing)+    setBranch (y : ys) _ = (ys, Just y)+    noChange ys z = (ys, Just z)++-- | Return node labels in pre-order.+labels :: Tree e a -> [a]+labels t = squish t []+  where+    squish (Node _ lb ts) xs = lb : foldr squish xs ts++-- | Set node labels in pre-order.+--+-- Return 'Nothing' if the provided list of node labels is too short.+setLabels :: Traversable t => [b] -> t a -> Maybe (t b)+setLabels xs = sequenceA . snd . mapAccumL setLabel xs+  where+    setLabel [] _ = ([], Nothing)+    setLabel (y : ys) _ = (ys, Just y)++-- | Label the nodes with unique integers starting at the root with 0.+identify :: Traversable t => t a -> t Int+identify = snd . mapAccumL (\i _ -> (i + 1, i)) (0 :: Int)++-- | The degree of the root node.+degree :: Tree e a -> Int+degree = (+ 1) . length . forest++-- | Prune degree two nodes.+--+-- The information stored in a pruned node is lost. The branches are combined+-- according to their 'Semigroup' instance of the form @\daughterBranch+-- parentBranch -> combinedBranch@.+prune :: Semigroup e => Tree e a -> Tree e a+prune t@(Node _ _ []) = t+prune (Node paBr _ [Node daBr daLb daTs]) = Node (daBr <> paBr) daLb daTs+prune (Node paBr paLb paTs) = Node paBr paLb $ map prune paTs++-- | Drop nodes satisfying predicate.+--+-- Degree two nodes may arise.+--+-- Also drop parent nodes of which all daughter nodes are dropped.+--+-- Return 'Nothing' if the root node satisfies the predicate.+dropNodesWith :: (a -> Bool) -> Tree e a -> Maybe (Tree e a)+dropNodesWith p (Node br lb ts)+  | p lb = Nothing+  | otherwise =+    if null ts'+      then Nothing+      else Just $ Node br lb ts'+  where+    ts' = mapMaybe (dropNodesWith p) ts++-- | Drop leaves satisfying predicate.+--+-- Degree two nodes may arise.+--+-- Also drop parent nodes of which all leaves are dropped.+--+-- Return 'Nothing' if all leaves satisfy the predicate.+dropLeavesWith :: (a -> Bool) -> Tree e a -> Maybe (Tree e a)+dropLeavesWith p (Node br lb [])+  | p lb = Nothing+  | otherwise = Just $ Node br lb []+dropLeavesWith p (Node br lb ts) =+  if null ts'+    then Nothing+    else Just $ Node br lb ts'+  where+    ts' = mapMaybe (dropLeavesWith p) ts++-- | Zip two trees with the same topology.+--+-- Return 'Nothing' if the topologies are different.+zipTreesWith ::+  (e1 -> e2 -> e) ->+  (a1 -> a2 -> a) ->+  Tree e1 a1 ->+  Tree e2 a2 ->+  Maybe (Tree e a)+zipTreesWith f g (Node brL lbL tsL) (Node brR lbR tsR) =+  if length tsL == length tsR+    then -- I am proud of that :)).+      zipWithM (zipTreesWith f g) tsL tsR >>= Just . Node (f brL brR) (g lbL lbR)+    else Nothing++-- | Zip two trees with the same topology.+--+-- Return 'Nothing' if the topologies are different.+zipTrees :: Tree e1 a1 -> Tree e2 a2 -> Maybe (Tree (e1, e2) (a1, a2))+zipTrees = zipTreesWith (,) (,)++duplicates :: Ord a => [a] -> Bool+duplicates = go S.empty+  where+    go _ [] = False+    go seen (x : xs) = x `S.member` seen || go (S.insert x seen) xs

@@ -0,0 +1,60 @@+-- |+-- Module      :  ELynx.Tree.Simulate.Coalescent+-- Description :  Generate coalescent trees+-- Copyright   :  (c) Dominik Schrempf 2018+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Wed May 16 13:13:11 2018.+module ELynx.Tree.Simulate.Coalescent+  ( simulate,+  )+where++import Control.Monad.Primitive+import ELynx.Tree.Measurable+import ELynx.Tree.Phylogeny+import ELynx.Tree.Rooted+import ELynx.Tree.Distribution.CoalescentContinuous+import Statistics.Distribution+import System.Random.MWC++-- | Simulate a coalescent tree with @n@ leaves. The branch lengths are in units+-- of effective population size.+simulate ::+  (PrimMonad m) =>+  -- | Number of leaves.+  Int ->+  Gen (PrimState m) ->+  m (Tree Length Int)+simulate n = simulate' n 0 trs+  where+    trs = [Node (Length 0) i [] | i <- [0 .. n - 1]]++simulate' ::+  (PrimMonad m) =>+  Int ->+  Int ->+  Forest Length Int ->+  Gen (PrimState m) ->+  m (Tree Length Int)+simulate' n a trs g+  | n <= 0 = error "Cannot construct trees without leaves."+  | n == 1 && length trs /= 1 = error "Too many trees provided."+  | n == 1 && length trs == 1 = return $ head trs+  | otherwise = do+    -- Indices of the leaves to join will be i-1 and i.+    i <- uniformR (1, n - 1) g+    -- The time of the coalescent event.+    t <- genContVar (coalescentDistributionCont n) g+    let trs' = map (applyStem (+ t)) trs -- Move time 't' up on the tree.+        tl = trs' !! (i - 1)+        tr = trs' !! i+        -- Join the two chosen trees.+        tm = Node (Length 0) a [tl, tr]+        -- Take the trees on the left, the merged tree, and the trees on the right.+        trs'' = take (i - 1) trs' ++ [tm] ++ drop (i + 1) trs'+    simulate' (n - 1) a trs'' g

@@ -0,0 +1,288 @@+{-# LANGUAGE BangPatterns #-}++-- |+--   Module      :  ELynx.Tree.Simulate.PointProcess+--   Description :  Point process and functions+--   Copyright   :  (c) Dominik Schrempf 2018+--   License     :  GPL-3.0-or-later+--+--   Maintainer  :  dominik.schrempf@gmail.com+--   Stability   :  unstable+--   Portability :  portable+--+-- Creation date: Tue Feb 13 13:16:18 2018.+--+-- See Gernhard, T. (2008). The conditioned reconstructed process. Journal of+-- Theoretical Biology, 253(4), 769–778. http://doi.org/10.1016/j.jtbi.2008.04.005.+--+-- The point process can be used to simulate reconstructed trees under the birth+-- and death process.+module ELynx.Tree.Simulate.PointProcess+  ( PointProcess (..),+    TimeSpec,+    simulate,+    toReconstructedTree,+    simulateReconstructedTree,+    simulateNReconstructedTrees,+  )+where++import Control.Monad+import Control.Monad.Primitive+import Data.Function+import Data.List+import Data.Sequence (Seq)+import qualified Data.Sequence as S+import ELynx.Tree.Distribution.BirthDeath+import ELynx.Tree.Distribution.BirthDeathCritical+import ELynx.Tree.Distribution.BirthDeathCriticalNoTime+import ELynx.Tree.Distribution.BirthDeathNearlyCritical+import ELynx.Tree.Distribution.TimeOfOrigin+import ELynx.Tree.Distribution.TimeOfOriginNearCritical+import ELynx.Tree.Distribution.Types+import ELynx.Tree.Measurable+import ELynx.Tree.Phylogeny+import ELynx.Tree.Rooted+import qualified Statistics.Distribution as D+  ( genContVar,+  )+import System.Random.MWC++-- Require near critical process if birth and death rates are closer than this value.+epsNearCriticalPointProcess :: Double+epsNearCriticalPointProcess = 1e-5++-- Also the distribution of origins needs a Tailor expansion for near critical values.+--+-- TODO: Check why the two epsilons are chosen differently.+epsNearCriticalTimeOfOrigin :: Double+epsNearCriticalTimeOfOrigin = 1e-8++-- Require critical process if birth and death rates are closer than this value.+eps :: Double+eps = 1e-12++(=~=) :: Double -> Double -> Bool+x =~= y = eps > abs (x - y)++-- Sort a list and also return original indices.+sortListWithIndices :: Ord a => [a] -> [(a, Int)]+sortListWithIndices xs = sortBy (compare `on` fst) $ zip xs ([0 ..] :: [Int])++-- Insert element into random position of list.+randomInsertList :: PrimMonad m => a -> [a] -> Gen (PrimState m) -> m [a]+randomInsertList e v g = do+  let l = length v+  i <- uniformR (0, l) g+  return $ take i v ++ [e] ++ drop i v++-- | A __point process__ for \(n\) points and of age \(t_{or}\) is defined as+-- follows. Draw $n$ points on the horizontal axis at \(1,2,\ldots,n\). Pick+-- \(n-1\) points at locations \((i+1/2, s_i)\), \(i=1,2,\ldots,n-1\);+-- \(0 < s_i < t_{or}\). There is a bijection between (ranked) oriented trees+-- and the point process. Usually, a will be 'String' (or 'Int') and b will be+-- 'Double'.+data PointProcess a b = PointProcess+  { points :: ![a],+    values :: ![b],+    origin :: !b+  }+  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)++-- | Sample a point process using the 'BirthDeathDistribution'. The names of the+-- points will be integers.+simulate ::+  (PrimMonad m) =>+  -- | Number of points (samples)+  Int ->+  -- | Time of origin or MRCA+  TimeSpec ->+  -- | Birth rate+  Rate ->+  -- | Death rate+  Rate ->+  -- | Generator (see 'System.Random.MWC')+  Gen (PrimState m) ->+  m (PointProcess Int Double)+-- 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+  | -- 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."+  | -- For the critical process, we have no idea about the time of origin, but can+    -- use a specially derived distribution.+    m =~= l =+    do+      !vs <- replicateM (n - 1) (D.genContVar (BDCNTD l) g)+      -- XXX: The length of the root branch will be 0.+      let t = maximum vs+      return $ PointProcess [0 .. (n - 1)] vs t+  | -- For the near critical process, we use a special distribution.+    abs (m - l) <= epsNearCriticalTimeOfOrigin =+    do+      t <- D.genContVar (TONCD n l m) g+      simulate n (Just (t, False)) 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+-- 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."+  | -- 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.+    -- -- | m < 0.0   = error "Death rate needs to be positive."+    -- Now, we have three different cases.+    -- 1. The critical branching process.+    -- 2. The near critical branching process.+    -- 3. Normal values :).+    (m =~= l) && not c = 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+    !vs <- replicateM (n - 1) (D.genContVar (BDNCD t l m) g)+    return $ PointProcess [0 .. (n - 1)] vs t+  | not c = do+    !vs <- replicateM (n - 1) (D.genContVar (BDD t l m) g)+    return $ PointProcess [0 .. (n - 1)] vs t+  | (m =~= l) && c = 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+    !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+    !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).+sortPP :: (Ord b) => PointProcess a b -> ([b], [Int])+sortPP (PointProcess _ vs _) = (vsSorted, isSorted)+  where+    vsIsSorted = sortListWithIndices vs+    vsSorted = map fst vsIsSorted+    isSorted = flattenIndices $ map snd vsIsSorted++-- Decrement indices that are above the one that is merged.+flattenIndices :: [Int] -> [Int]+flattenIndices is = snd $ mapAccumL fAcc [] is++-- TODO: This is the bottleneck for simulating large trees.+--+-- The accumulating function. Count the number of indices which are before the+-- current index and lower than the current index.+fAcc :: [Int] -> Int -> ([Int], Int)+fAcc is i = (i : is, i') where i' = i - length (filter (< i) is)++-- | See 'simulateReconstructedTree', but n times.+simulateNReconstructedTrees ::+  (PrimMonad m) =>+  -- | Number of trees+  Int ->+  -- | Number of points (samples)+  Int ->+  -- | Time of origin or MRCA+  TimeSpec ->+  -- | Birth rate+  Rate ->+  -- | Death rate+  Rate ->+  -- | Generator (see 'System.Random.MWC')+  Gen (PrimState m) ->+  m (Forest Length Int)+simulateNReconstructedTrees nT nP t l m g+  | nT <= 0 = return []+  | otherwise = replicateM nT $ simulateReconstructedTree nP t l m g++-- | Use the point process to simulate a reconstructed tree (see+-- 'toReconstructedTree') possibly with specific height and a fixed number of+-- leaves according to the birth and death process.+simulateReconstructedTree ::+  (PrimMonad m) =>+  -- | Number of points (samples)+  Int ->+  -- | Time of origin or MRCA+  TimeSpec ->+  -- | Birth rate+  Rate ->+  -- | Death rate+  Rate ->+  -- | Generator (see 'System.Random.MWC')+  Gen (PrimState m) ->+  m (Tree Length Int)+simulateReconstructedTree n t l m g =+  toReconstructedTree 0 <$> simulate n t l m g++-- | Convert a point process to a reconstructed tree. See Lemma 2.2.++-- Of course, I decided to only use one tree structure with extinct and extant+-- leaves (actually a complete tree). So a tree created here just does not+-- contain extinct leaves. A function 'isReconstructed' is provided to test if a+-- tree is reconstructed (and not complete) in this sense. However, a complete+-- tree might show up as "reconstructed", just because, by chance, it does not+-- contain extinct leaves. I wanted to use a Monoid constraint to get the unit+-- element, but this fails for classical 'Int's. So, I rather have another+-- (useless) argument.+toReconstructedTree ::+  a -> -- Default node label.+  PointProcess a Double ->+  Tree Length a+toReconstructedTree l pp@(PointProcess ps vs o)+  | length ps /= length vs + 1 = error "Too few or too many points."+  | length vs <= 1 = error "Too few values."+  | -- -- Test is deactivated.+    -- -- | otherwise = if isReconstructed treeOrigin then treeOrigin else error "Error in algorithm."+    otherwise =+    treeOrigin+  where+    (vsSorted, isSorted) = sortPP pp+    !lvs = S.fromList [Node (Length 0) p [] | p <- ps]+    !heights = S.replicate (length ps) 0+    !treeRoot = toReconstructedTree' isSorted vsSorted l lvs heights+    !h = last vsSorted+    !treeOrigin = applyStem (+ (o - h)) treeRoot++-- Move up the tree, connect nodes when they join according to the point process.+toReconstructedTree' ::+  [Int] -> -- Sorted indices, see 'sort'.+  [Double] -> -- Sorted merge values.+  a -> -- Default node label.+  Seq (Tree Length a) -> -- Leaves with accumulated root branch lengths.+  Seq Double -> -- Accumulated heights of the leaves.+  Tree Length a+toReconstructedTree' [] [] _ trs _ = trs `S.index` 0+toReconstructedTree' is vs l trs hs = toReconstructedTree' is' vs' l trs'' hs'+  where+    -- For the algorithm, see 'simulate' but index starts at zero.++    !i = head is+    !is' = tail is+    !v = head vs+    !vs' = tail vs+    -- Left: l, right: r.+    !hl = hs `S.index` i+    !hr = hs `S.index` (i + 1)+    !dvl = v - hl+    !dvr = v - hr+    !tl = applyStem (+ dvl) $ trs `S.index` i+    !tr = applyStem (+ dvr) $ trs `S.index` (i + 1)+    !h' = hl + dvl -- Should be the same as 'hr + dvr'.+    !tm = Node (Length 0) l [tl, tr]+    !trs'' = (S.take i trs S.|> tm) S.>< S.drop (i + 2) trs+    !hs' = (S.take i hs S.|> h') S.>< S.drop (i + 2) hs