diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -1,11 +1,18 @@ -# Changelog for ELynx +# Revision history for ELynx ## Unreleased changes +- Improve rooting functions. +- Improve `Topology` data type (but still a lot to do). +- Various additions to the documentation. +- Rename `Measurable` to `HasLength`, `Supported` to `HasSupport`, and `Named` + to `HasLength`. +- Cabal and stack file changes. -## Version 0.5.0 + +## Version 0.5.0.1 - `modLen`, `modSup`. - Newtype wrappers for branch length, branch support, and node name. Those data diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -2,7 +2,7 @@ # The ELynx Suite -Version: 0.5.0. +Version: 0.5.0.1. Reproducible evolution made easy.

@@ -19,18 +19,20 @@ The library packages are: -- **elynx-nexus:** Nexus file support. -- **elynx-markov:** Simulate multi sequence alignments along phylogenetic trees. -- **elynx-seq:** Handle evolutionary sequences and multi sequence alignments. -- **elynx-tools:** Tools for the provided executables. -- **elynx-tree:** Handle phylogenetic trees. +- **[elynx-nexus](https://hackage.haskell.org/package/elynx-nexus):** Nexus file support. +- **[elynx-markov](https://hackage.haskell.org/package/elynx-markov):** Simulate multi sequence alignments along phylogenetic trees. +- **[elynx-seq](https://hackage.haskell.org/package/elynx-seq):** Handle evolutionary sequences and multi sequence alignments. +- **[elynx-tools](https://hackage.haskell.org/package/elynx-tools):** Tools for the provided executables. +- **[elynx-tree](https://hackage.haskell.org/package/elynx-tree):** Handle phylogenetic trees. The executables are: -- **slynx:** Analyze, modify, and simulate evolutionary sequences. -- **tlynx:** Analyze, modify, and simulate phylogenetic trees. -- **elynx:** Validate and redo past analyses. +- **[slynx](https://hackage.haskell.org/package/slynx):** Analyze, modify, and simulate evolutionary sequences. +- **[tlynx](https://hackage.haskell.org/package/tlynx):** Analyze, modify, and simulate phylogenetic trees. +- **[elynx](https://hackage.haskell.org/package/elynx):** Validate and redo past analyses. +Documentation is available on [Hackage](https://hackage.haskell.org/) (use direct links above). + **ELynx is actively developed. We happily receive comments, ideas, feature requests, and pull requests!** @@ -65,34 +67,15 @@ [PATH](https://en.wikipedia.org/wiki/PATH_(variable)) environment variable. Then, they can be used directly. -# Documentation - -Documentation is available on [Hackage](https://hackage.haskell.org/). - -Libraries: - -- [elynx-nexus](https://hackage.haskell.org/package/elynx-nexus) -- [elynx-markov](https://hackage.haskell.org/package/elynx-markov) -- [elynx-seq](https://hackage.haskell.org/package/elynx-seq) -- [elynx-tools](https://hackage.haskell.org/package/elynx-tools) -- [elynx-tree](https://hackage.haskell.org/package/elynx-tree) - -Executables: - -- [elynx](https://hackage.haskell.org/package/elynx) -- [slynx](https://hackage.haskell.org/package/slynx) -- [tlynx](https://hackage.haskell.org/package/tlynx) - - # SLynx Handle evolutionary sequences. slynx --help | head -n -16 - ELynx Suite version 0.5.0. + ELynx Suite version 0.5.0.1. Developed by Dominik Schrempf. - Compiled on November 10, 2020, at 14:29 pm, UTC. + Compiled on December 18, 2020, at 10:14 am, UTC. Usage: slynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME] [-f|--force] [--no-elynx-file] COMMAND @@ -137,9 +120,9 @@ slynx concatenate --help - ELynx Suite version 0.5.0. + ELynx Suite version 0.5.0.1. Developed by Dominik Schrempf. - Compiled on November 10, 2020, at 14:29 pm, UTC. + Compiled on December 18, 2020, at 10:14 am, UTC. Usage: slynx concatenate (-a|--alphabet NAME) INPUT-FILE Concatenate sequences found in input files. @@ -158,9 +141,9 @@ slynx examine --help - ELynx Suite version 0.5.0. + ELynx Suite version 0.5.0.1. Developed by Dominik Schrempf. - Compiled on November 10, 2020, at 14:29 pm, UTC. + Compiled on December 18, 2020, at 10:14 am, UTC. Usage: slynx examine (-a|--alphabet NAME) INPUT-FILE [--per-site] Examine sequences. If data is a multi sequence alignment, additionally analyze columns. @@ -180,9 +163,9 @@ slynx filter-rows --help - ELynx Suite version 0.5.0. + ELynx Suite version 0.5.0.1. Developed by Dominik Schrempf. - Compiled on November 10, 2020, at 14:29 pm, UTC. + Compiled on December 18, 2020, at 10:14 am, UTC. Usage: slynx filter-rows (-a|--alphabet NAME) INPUT-FILE [--longer-than LENGTH] [--shorter-than LENGTH] [--standard-characters] @@ -203,9 +186,9 @@ slynx filter-columns --help - ELynx Suite version 0.5.0. + ELynx Suite version 0.5.0.1. Developed by Dominik Schrempf. - Compiled on November 10, 2020, at 14:29 pm, UTC. + Compiled on December 18, 2020, at 10:14 am, UTC. Usage: slynx filter-columns (-a|--alphabet NAME) INPUT-FILE [--standard-chars DOUBLE] @@ -227,9 +210,9 @@ slynx simulate --help - ELynx Suite version 0.5.0. + ELynx Suite version 0.5.0.1. Developed by Dominik Schrempf. - Compiled on November 10, 2020, at 14:29 pm, UTC. + Compiled on December 18, 2020, at 10:14 am, UTC. Usage: slynx simulate (-t|--tree-file Name) [-s|--substitution-model MODEL] [-m|--mixture-model MODEL] [-e|--edm-file NAME] @@ -306,9 +289,9 @@ slynx sub-sample --help - ELynx Suite version 0.5.0. + ELynx Suite version 0.5.0.1. Developed by Dominik Schrempf. - Compiled on November 10, 2020, at 14:29 pm, UTC. + Compiled on December 18, 2020, at 10:14 am, UTC. Usage: slynx sub-sample (-a|--alphabet NAME) INPUT-FILE (-n|--number-of-sites INT) @@ -336,9 +319,9 @@ slynx translate --help - ELynx Suite version 0.5.0. + ELynx Suite version 0.5.0.1. Developed by Dominik Schrempf. - Compiled on November 10, 2020, at 14:29 pm, UTC. + Compiled on December 18, 2020, at 10:14 am, UTC. Usage: slynx translate (-a|--alphabet NAME) INPUT-FILE (-r|--reading-frame INT) (-u|--universal-code CODE) @@ -361,9 +344,9 @@ tlynx --help | head -n -16 - ELynx Suite version 0.5.0. + ELynx Suite version 0.5.0.1. Developed by Dominik Schrempf. - Compiled on November 10, 2020, at 14:29 pm, UTC. + Compiled on December 18, 2020, at 10:14 am, UTC. Usage: tlynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME] [-f|--force] [--no-elynx-file] COMMAND @@ -401,9 +384,9 @@ tlynx compare --help - ELynx Suite version 0.5.0. + ELynx Suite version 0.5.0.1. Developed by Dominik Schrempf. - Compiled on November 10, 2020, at 14:29 pm, UTC. + Compiled on December 18, 2020, at 10:14 am, UTC. Usage: tlynx compare [-n|--normalize] [-b|--bipartitions] [-t|--intersect] [-f|--newick-format FORMAT] NAMES @@ -430,9 +413,9 @@ tlynx examine --help - ELynx Suite version 0.5.0. + ELynx Suite version 0.5.0.1. Developed by Dominik Schrempf. - Compiled on November 10, 2020, at 14:29 pm, UTC. + Compiled on December 18, 2020, at 10:14 am, UTC. Usage: tlynx examine INPUT-FILE [-f|--newick-format FORMAT] Compute summary statistics of phylogenetic trees. @@ -454,9 +437,9 @@ tlynx simulate --help - ELynx Suite version 0.5.0. + ELynx Suite version 0.5.0.1. Developed by Dominik Schrempf. - Compiled on November 10, 2020, at 14:29 pm, UTC. + Compiled on December 18, 2020, at 10:14 am, UTC. Usage: tlynx simulate (-t|--nTrees INT) (-n|--nLeaves INT) PROCESS [-u|--sub-sample DOUBLE] [-s|--summary-statistics] @@ -491,9 +474,9 @@ elynx --help | head -n -16 - ELynx Suite version 0.5.0. + ELynx Suite version 0.5.0.1. Developed by Dominik Schrempf. - Compiled on November 10, 2020, at 14:29 pm, UTC. + Compiled on December 18, 2020, at 10:14 am, UTC. Usage: elynx COMMAND Validate and redo past ELynx analyses diff --git a/bench/Length.hs b/bench/Length.hs --- a/bench/Length.hs +++ b/bench/Length.hs @@ -19,10 +19,10 @@ where import Data.Foldable -import ELynx.Tree.Measurable +import ELynx.Tree.Length lengthSumFoldl' :: [Length] -> Length -lengthSumFoldl' = foldl' (\x y -> either error id $ toLength $ fromLength x + fromLength y) 0 +lengthSumFoldl' = foldl' (\x y -> toLength "lengthSumFoldl'" $ fromLength x + fromLength y) 0 lengthSumFoldl'Unsafe :: [Length] -> Length lengthSumFoldl'Unsafe = foldl' (\x y -> toLengthUnsafe $ fromLength x + fromLength y) 0 diff --git a/bench/Lens.hs b/bench/Lens.hs --- a/bench/Lens.hs +++ b/bench/Lens.hs @@ -18,20 +18,20 @@ where import Data.Foldable -import ELynx.Tree.Measurable +import ELynx.Tree.Length import Lens.Micro -len :: Measurable a => Lens' a Length +len :: HasLength a => Lens' a Length len = lens getLen (flip setLen) -sumWithGetter :: Measurable a => [a] -> Length +sumWithGetter :: HasLength a => [a] -> Length sumWithGetter = foldl' (\x y -> x ^. len + y ^. len) 0 -sumWithSetter :: Measurable a => [a] -> Length +sumWithSetter :: HasLength a => [a] -> Length sumWithSetter = sumWithGetter . map (\x -> x & len %~ (+ 10)) -sumWithAccessorFunction :: Measurable a => [a] -> Length +sumWithAccessorFunction :: HasLength a => [a] -> Length sumWithAccessorFunction = foldl' (\x y -> getLen x + getLen y) 0 -sumWithModifyFunction :: Measurable a => [a] -> Length +sumWithModifyFunction :: HasLength a => [a] -> Length sumWithModifyFunction = sumWithAccessorFunction . map (modLen (+ 10)) diff --git a/elynx-tree.cabal b/elynx-tree.cabal --- a/elynx-tree.cabal +++ b/elynx-tree.cabal @@ -1,18 +1,19 @@ -cabal-version: 2.2 -name: elynx-tree -version: 0.5.0 -synopsis: Handle phylogenetic trees -description: Examine, compare, and simulate phylogenetic trees in a reproducible way. Please see the README on GitHub at . -category: Bioinformatics -homepage: https://github.com/dschrempf/elynx#readme -bug-reports: https://github.com/dschrempf/elynx/issues -author: Dominik Schrempf -maintainer: dominik.schrempf@gmail.com -copyright: Dominik Schrempf (2020) -license: GPL-3.0-or-later -license-file: LICENSE -build-type: Simple +cabal-version: 2.2 +name: elynx-tree +version: 0.5.0.1 +license: GPL-3.0-or-later +license-file: LICENSE +copyright: Dominik Schrempf (2020) +maintainer: dominik.schrempf@gmail.com +author: Dominik Schrempf +homepage: https://github.com/dschrempf/elynx#readme +bug-reports: https://github.com/dschrempf/elynx/issues +synopsis: Handle phylogenetic trees +description: + Examine, compare, and simulate phylogenetic trees in a reproducible way. Please see the README on GitHub at . +category: Bioinformatics +build-type: Simple extra-source-files: README.md ChangeLog.md @@ -34,114 +35,111 @@ data/UltraMetric.tree source-repository head - type: git - location: https://github.com/dschrempf/elynx + type: git + location: https://github.com/dschrempf/elynx library - exposed-modules: - ELynx.Topology - ELynx.Topology.Phylogeny - ELynx.Topology.Rooted - ELynx.Tree - ELynx.Tree.Bipartition - ELynx.Tree.Distance - ELynx.Tree.Measurable - ELynx.Tree.Named - ELynx.Tree.Parallel - 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: - Paths_elynx_tree - hs-source-dirs: - src - ghc-options: -Wall - build-depends: - aeson - , attoparsec - , base >=4.7 && <5 - , bytestring - , comonad - , containers - , deepseq - , double-conversion - , elynx-nexus - , math-functions - , mwc-random - , parallel - , primitive - , statistics - default-language: Haskell2010 + exposed-modules: + ELynx.Topology + ELynx.Topology.Phylogeny + ELynx.Topology.Rooted + ELynx.Tree + ELynx.Tree.Bipartition + ELynx.Tree.Distance + ELynx.Tree.Length + ELynx.Tree.Name + ELynx.Tree.Parallel + ELynx.Tree.Partition + ELynx.Tree.Phylogeny + ELynx.Tree.Rooted + ELynx.Tree.Splittable + ELynx.Tree.Support + 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 + hs-source-dirs: src + other-modules: Paths_elynx_tree + autogen-modules: Paths_elynx_tree + default-language: Haskell2010 + ghc-options: -Wall -Wunused-packages + build-depends: + aeson >=1.5.4.1 && <1.6, + attoparsec >=0.13.2.4 && <0.14, + base >=4.14.1.0 && <4.15, + bytestring >=0.10.10.0 && <0.11, + comonad >=5.0.6 && <5.1, + containers >=0.6.2.1 && <0.7, + deepseq >=1.4.4.0 && <1.5, + double-conversion >=2.0.2.0 && <2.1, + elynx-nexus >=0.5.0.1 && <0.6, + math-functions >=0.3.4.1 && <0.4, + mwc-random >=0.14.0.0 && <0.15, + parallel >=3.2.2.0 && <3.3, + primitive >=0.7.1.0 && <0.8, + statistics >=0.15.2.0 && <0.16 + test-suite tree-test - type: exitcode-stdio-1.0 - main-is: Spec.hs - other-modules: - 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 - ghc-options: -Wall - build-depends: - QuickCheck - , attoparsec - , base >=4.7 && <5 - , bytestring - , containers - , elynx-tools - , elynx-tree - , hspec - default-language: Haskell2010 + type: exitcode-stdio-1.0 + main-is: Spec.hs + hs-source-dirs: test + other-modules: + ELynx.Topology.RootedSpec + ELynx.Tree.Arbitrary + ELynx.Tree.BipartitionSpec + ELynx.Tree.DistanceSpec + ELynx.Tree.PartitionSpec + ELynx.Tree.PhylogenySpec + ELynx.Tree.RootedSpec + ELynx.Tree.SupportSpec + ELynx.Tree.Export.NewickSpec + ELynx.Tree.Export.NexusSpec + ELynx.Tree.Import.NewickSpec + ELynx.Tree.Import.NexusSpec + Paths_elynx_tree + default-language: Haskell2010 + ghc-options: -Wall -Wunused-packages + build-depends: + QuickCheck >=2.13.2 && <2.14, + attoparsec >=0.13.2.4 && <0.14, + base >=4.14.1.0 && <4.15, + bytestring >=0.10.10.0 && <0.11, + containers >=0.6.2.1 && <0.7, + elynx-tools >=0.5.0.1 && <0.6, + elynx-tree -any, + hspec >=2.7.4 && <2.8 + benchmark tree-bench - type: exitcode-stdio-1.0 - main-is: Bench.hs - other-modules: - Length - Lens - Paths_elynx_tree - hs-source-dirs: - bench - ghc-options: -Wall -threaded -rtsopts -with-rtsopts=-N - -- ghc-options: -Wall - build-depends: - base >=4.7 && <5 - , bytestring - , criterion - , deepseq - , elynx-tools - , elynx-tree - , microlens - , mwc-random - , parallel - default-language: Haskell2010 + type: exitcode-stdio-1.0 + main-is: Bench.hs + hs-source-dirs: bench + other-modules: + Length + Lens + Paths_elynx_tree + + default-language: Haskell2010 + ghc-options: + -Wall -Wunused-packages -threaded -rtsopts -with-rtsopts=-N + + build-depends: + base >=4.7 && <5, + criterion >=1.5.7.0 && <1.6, + elynx-tools >=0.5.0.1 && <0.6, + elynx-tree -any, + microlens >=0.4.11.2 && <0.5, + mwc-random >=0.14.0.0 && <0.15, + parallel >=3.2.2.0 && <3.3 diff --git a/src/ELynx/Topology/Phylogeny.hs b/src/ELynx/Topology/Phylogeny.hs --- a/src/ELynx/Topology/Phylogeny.hs +++ b/src/ELynx/Topology/Phylogeny.hs @@ -10,6 +10,8 @@ -- -- Creation date: Sat Jul 18 13:15:49 2020. -- +-- THIS MODULE IS INCOMPLETE. +-- -- 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. @@ -23,93 +25,54 @@ -- 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. +-- Note: 'Topology's are rooted. -- --- The bifurcating root of a topology can be changed with 'roots', or 'rootAt'. +-- Note: 'Topology's 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 or more +-- children. This convention is not enforced here. Newick topologies are just +-- parsed as they are, and a rooted topology is returned. -- --- Topologies with multifurcating root nodes can be properly rooted using --- 'outgroup'. +-- THIS MODULE IS INCOMPLETE. module ELynx.Topology.Phylogeny - ( outgroup, - roots, - rootAt, + ( equal, + equal', ) where -import Data.Set (Set) +import Data.List +import Data.Maybe 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." +-- | The equality check is slow because the order of children is considered to +-- be arbitrary. +-- +-- NOTE: The equality check is only meaningful if the topologies have unique +-- leaves. +-- +-- Return 'Left' if a topology does not have unique leaves. +equal :: (Eq a, Ord a) => Topology a -> Topology a -> Either String Bool +equal tL tR + | duplicateLeaves tL = Left "equal: Left topology has duplicate leaves." + | duplicateLeaves tR = Left "equal: Right topology has duplicate leaves." + | otherwise = Right $ equal' tL tR --- | For a rooted topology with a bifurcating root node, get all possible rooted --- topologies. -roots :: Topology a -> Either String (Forest a) -roots = undefined +-- | Same as 'equal', but assume that leaves are unique. +equal' :: Eq a => Topology a -> Topology a -> Bool +equal' (Leaf lbL) (Leaf lbR) = lbL == lbR +equal' (Node tsL) (Node tsR) = + (length tsL == length tsR) + && all (`elem'` tsR) tsL + where + elem' t ts = isJust $ find (equal' t) ts +equal' _ _ = False --- -- | 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 +-- TODO. --- | Root a tree at a specific position. -rootAt :: Ord a => Bipartition a -> Topology a -> Either String (Forest a) -rootAt = undefined +-- A multifurcating root node can be resolved to a bifurcating root node with +-- 'outgroup'. +-- +-- The bifurcating root node can be changed with 'outgroup' or 'midpoint'. +-- +-- For a given topology with bifurcating root node, a list of all rooted +-- topologies is returned by 'roots'. diff --git a/src/ELynx/Topology/Rooted.hs b/src/ELynx/Topology/Rooted.hs --- a/src/ELynx/Topology/Rooted.hs +++ b/src/ELynx/Topology/Rooted.hs @@ -13,15 +13,15 @@ -- -- 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. +-- THIS MODULE IS INCOMPLETE. -- --- For rooted trees, please see 'ELynx.Tree.Rooted'. +-- A rooted 'Topology' differs from a classical rooted rose 'Data.Tree.Tree' in +-- that it does not have internal node labels. The leaves have labels. -- --- 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. +-- For rooted trees with branch labels, please see "ELynx.Tree.Rooted". Please +-- also see the note about tree traversals therein. +-- +-- THIS MODULE IS INCOMPLETE. module ELynx.Topology.Rooted ( -- * Data type Topology (..), @@ -29,16 +29,17 @@ fromTree, fromLabeledTree, - -- * Functions - degree, + -- * Access leaves, branches and labels leaves, - flatten, + duplicateLeaves, identify, + + -- * Structure + degree, prune, dropLeavesWith, zipTreesWith, zipTrees, - duplicateLeaves, ) where @@ -80,7 +81,7 @@ null _ = False {-# INLINE null #-} - toList = flatten + toList = leaves {-# INLINE toList #-} instance Traversable Topology where @@ -120,23 +121,6 @@ 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. @@ -150,10 +134,40 @@ fromLabeledTree (R.Node _ lb []) = Leaf lb fromLabeledTree (R.Node _ _ xs) = Node $ fromLabeledTree <$> N.fromList xs +-- TODO: Maybe use foldr similar to 'flatten'. +-- | Set of leaves. +leaves :: Topology a -> [a] +leaves (Leaf lb) = [lb] +leaves (Node ts) = concatMap leaves ts + +-- -- TODO: Check if this implementation of 'leaves' is faster. +-- -- | 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 + +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 + +-- TODO: This is the same as in ELynx.Tree.Rooted. -- | 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) +-- | The degree of the root node. +degree :: Topology a -> Int +degree (Node ts) = (+ 1) $ length ts +degree (Leaf _) = 1 + -- | Prune degree two nodes. prune :: Topology a -> Topology a prune (Node ts) @@ -195,13 +209,3 @@ -- 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 diff --git a/src/ELynx/Tree.hs b/src/ELynx/Tree.hs --- a/src/ELynx/Tree.hs +++ b/src/ELynx/Tree.hs @@ -22,12 +22,12 @@ module ELynx.Tree.Parallel, -- * Branch label classes - module ELynx.Tree.Measurable, + module ELynx.Tree.Length, module ELynx.Tree.Splittable, - module ELynx.Tree.Supported, + module ELynx.Tree.Support, -- * Node label classes - module ELynx.Tree.Named, + module ELynx.Tree.Name, -- * Phylogenies module ELynx.Tree.Phylogeny, @@ -51,14 +51,14 @@ 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.Length +import ELynx.Tree.Name import ELynx.Tree.Parallel import ELynx.Tree.Partition import ELynx.Tree.Phylogeny import ELynx.Tree.Rooted import ELynx.Tree.Splittable -import ELynx.Tree.Supported +import ELynx.Tree.Support import ELynx.Tree.Zipper -- -- | An evolutionary label has some information about where the corresponding diff --git a/src/ELynx/Tree/Distance.hs b/src/ELynx/Tree/Distance.hs --- a/src/ELynx/Tree/Distance.hs +++ b/src/ELynx/Tree/Distance.hs @@ -36,7 +36,7 @@ import Data.Set (Set) import qualified Data.Set as S import ELynx.Tree.Bipartition -import ELynx.Tree.Measurable +import ELynx.Tree.Length import ELynx.Tree.Partition import ELynx.Tree.Rooted @@ -124,7 +124,7 @@ -- trees is returned. -- -- XXX: Comparing a list of trees recomputes bipartitions. -branchScore :: (Measurable e1, Measurable e2, Ord a) => Tree e1 a -> Tree e2 a -> Either String Double +branchScore :: (HasLength e1, HasLength 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 diff --git a/src/ELynx/Tree/Export/Newick.hs b/src/ELynx/Tree/Export/Newick.hs --- a/src/ELynx/Tree/Export/Newick.hs +++ b/src/ELynx/Tree/Export/Newick.hs @@ -23,11 +23,11 @@ import qualified Data.ByteString.Builder as BB import qualified Data.ByteString.Lazy.Char8 as BL import Data.List (intersperse) -import ELynx.Tree.Measurable -import ELynx.Tree.Named +import ELynx.Tree.Length +import ELynx.Tree.Name import ELynx.Tree.Phylogeny import ELynx.Tree.Rooted -import ELynx.Tree.Supported +import ELynx.Tree.Support buildBrLen :: Length -> BB.Builder buildBrLen bl = BB.char8 ':' <> BB.doubleDec (fromLength bl) @@ -36,7 +36,7 @@ buildBrSup bs = BB.char8 '[' <> BB.doubleDec (fromSupport bs) <> BB.char8 ']' -- | See 'toNewick'. -toNewickBuilder :: Named a => Tree Phylo a -> BB.Builder +toNewickBuilder :: HasName a => Tree Phylo a -> BB.Builder toNewickBuilder t = go t <> BB.char8 ';' where go (Node b l []) = lbl b l @@ -65,5 +65,5 @@ -- @ -- "ACTUALNAME[posterior=-2839.2,age_95%_HPD={4.80804,31.6041}]" -- @ -toNewick :: Named a => Tree Phylo a -> BL.ByteString +toNewick :: HasName a => Tree Phylo a -> BL.ByteString toNewick = BB.toLazyByteString . toNewickBuilder diff --git a/src/ELynx/Tree/Export/Nexus.hs b/src/ELynx/Tree/Export/Nexus.hs --- a/src/ELynx/Tree/Export/Nexus.hs +++ b/src/ELynx/Tree/Export/Nexus.hs @@ -19,13 +19,13 @@ import qualified Data.ByteString.Lazy.Char8 as BL import ELynx.Export.Nexus import ELynx.Tree.Export.Newick -import ELynx.Tree.Named +import ELynx.Tree.Name 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 :: HasName 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 :: HasName a => (BL.ByteString, Tree Phylo a) -> BL.ByteString tree (n, t) = " TREE " <> n <> " = " <> toNewick t diff --git a/src/ELynx/Tree/Import/Newick.hs b/src/ELynx/Tree/Import/Newick.hs --- a/src/ELynx/Tree/Import/Newick.hs +++ b/src/ELynx/Tree/Import/Newick.hs @@ -31,7 +31,7 @@ -- @ module ELynx.Tree.Import.Newick ( NewickFormat (..), - description, + describeNewickFormat, newick, parseNewick, oneNewick, @@ -46,11 +46,11 @@ import Data.Attoparsec.ByteString.Char8 import qualified Data.ByteString.Char8 as BS import qualified Data.ByteString.Lazy.Char8 as BL -import ELynx.Tree.Measurable -import ELynx.Tree.Named +import ELynx.Tree.Length +import ELynx.Tree.Name import ELynx.Tree.Phylogeny import ELynx.Tree.Rooted hiding (forest, label) -import ELynx.Tree.Supported +import ELynx.Tree.Support import GHC.Generics import Prelude hiding (takeWhile) @@ -68,12 +68,12 @@ instance ToJSON NewickFormat -- | Short description of the supported Newick formats. -description :: NewickFormat -> String -description Standard = +describeNewickFormat :: NewickFormat -> String +describeNewickFormat Standard = "Standard: Branch support values are stored in square brackets after branch lengths." -description IqTree = +describeNewickFormat IqTree = "IqTree: Branch support values are stored as node names after the closing bracket of forests." -description RevBayes = +describeNewickFormat 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." -- | Newick tree parser. Also succeeds when more trees follow. @@ -155,7 +155,7 @@ branchLength = do _ <- char ':' "branchLengthDelimiter" l <- double "branchLength" - return $ either error id $ toLength l + return $ toLength "branchLength" l branchSupport :: Parser Support branchSupport = @@ -163,7 +163,7 @@ _ <- char '[' "branchSupportBegin" s <- double "branchSupport" _ <- char ']' "branchSupportEnd" - return $ either error id $ toSupport s + return $ toSupport "branchSupport" s -------------------------------------------------------------------------------- -- IQ-TREE. @@ -191,7 +191,7 @@ branchedIqTree = ( "branchedIqTree") $ do f <- forestIqTree ms <- optional double - let s = either error id . toSupport <$> ms + let s = toSupport "branchedIqTree" <$> ms n <- name b <- optional branchLength return $ Node (Phylo b s) n f diff --git a/src/ELynx/Tree/Import/Nexus.hs b/src/ELynx/Tree/Import/Nexus.hs --- a/src/ELynx/Tree/Import/Nexus.hs +++ b/src/ELynx/Tree/Import/Nexus.hs @@ -21,7 +21,7 @@ import qualified Data.ByteString.Char8 as BS import ELynx.Import.Nexus import ELynx.Tree.Import.Newick -import ELynx.Tree.Named +import ELynx.Tree.Name import ELynx.Tree.Phylogeny import ELynx.Tree.Rooted import Prelude hiding (takeWhile) diff --git a/src/ELynx/Tree/Length.hs b/src/ELynx/Tree/Length.hs new file mode 100644 --- /dev/null +++ b/src/ELynx/Tree/Length.hs @@ -0,0 +1,172 @@ +{-# LANGUAGE DeriveAnyClass #-} +{-# LANGUAGE DeriveGeneric #-} +{-# LANGUAGE DerivingVia #-} + +-- | +-- Module : ELynx.Tree.Length +-- Description : Measurable 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 lengths is not completely ensured. See the documentation of +-- 'Length'. +module ELynx.Tree.Length + ( -- * Non-negative length + Length (fromLength), + toLength, + toLengthUnsafe, + checkLength, + HasLength (..), + height, + rootHeight, + + -- * Functions on trees + distancesOriginLeaves, + totalBranchLength, + normalizeBranchLengths, + normalizeHeight, + ultrametric, + makeUltrametric, + ) +where + +import Control.DeepSeq +import Data.Aeson +import Data.Bifoldable +import Data.Bifunctor +import Data.Semigroup +import ELynx.Tree.Rooted +import ELynx.Tree.Splittable +import GHC.Generics + +-- | Non-negative length. +-- +-- However, non-negativity is only checked with 'toLength', and negative values +-- can be obtained using the 'Num' and related instances. +-- +-- Safe conversion is roughly 50 percent slower. +-- +-- @ +-- benchmarking length/length sum foldl' with safe conversion +-- time 110.4 ms (109.8 ms .. 111.0 ms) +-- 1.000 R² (1.000 R² .. 1.000 R²) +-- mean 110.2 ms (110.0 ms .. 110.6 ms) +-- std dev 501.8 μs (359.1 μs .. 730.0 μs) +-- +-- benchmarking length/length sum foldl' num instance +-- time 89.37 ms (85.13 ms .. 94.27 ms) +-- 0.996 R² (0.992 R² .. 1.000 R²) +-- mean 86.53 ms (85.63 ms .. 88.52 ms) +-- std dev 2.239 ms (1.069 ms .. 3.421 ms) +-- +-- benchmarking length/double sum foldl' +-- time 85.47 ms (84.88 ms .. 86.42 ms) +-- 1.000 R² (0.999 R² .. 1.000 R²) +-- mean 85.56 ms (85.26 ms .. 86.02 ms) +-- std dev 611.9 μs (101.5 μs .. 851.7 μs) +-- @ +newtype Length = Length {fromLength :: Double} + deriving (Read, Show, Generic, NFData) + deriving (Enum, Eq, Floating, Fractional, Num, Ord, Real, RealFloat, RealFrac) via Double + deriving (Semigroup, Monoid) via Sum Double + +instance Splittable Length where + split = (/ 2.0) + +instance ToJSON Length + +instance FromJSON Length + +instance HasLength Length where + getLen = id + setLen = const + modLen f = f + +-- | If negative, call 'error' indicating the calling function name. +toLength :: String -> Double -> Length +toLength s x + | x < 0 = error $ s ++ ": Length is negative: " ++ show x ++ "." + | otherwise = Length x + +-- | Do not check if value is negative. +toLengthUnsafe :: Double -> Length +toLengthUnsafe = Length + +-- | If negative, call 'error' with given calling function name. +checkLength :: String -> Length -> Length +checkLength s = toLength s . fromLength + +-- | A data type with measurable and modifiable values. +class HasLength e where + -- | Get length. + getLen :: e -> Length + + -- | Set length. + setLen :: Length -> e -> e + + -- For computational efficiency. + + -- | Modify length. + modLen :: (Length -> Length) -> e -> e + +-- | The maximum distance between origin and leaves. +-- +-- The height includes the branch length of the stem. +height :: HasLength e => Tree e a -> Length +height = maximum . distancesOriginLeaves + +-- | The maximum distance between root node and leaves. +rootHeight :: HasLength e => Tree e a -> Length +rootHeight (Node _ _ []) = 0 +rootHeight t = maximum $ concatMap distancesOriginLeaves (forest t) + +-- | Distances from the origin of a tree to the leaves. +-- +-- The distances include the branch length of the stem. +distancesOriginLeaves :: HasLength e => Tree e a -> [Length] +distancesOriginLeaves (Node br _ []) = [getLen br] +distancesOriginLeaves (Node br _ ts) = map (getLen br +) (concatMap distancesOriginLeaves ts) + +-- | Total branch length of a tree. +totalBranchLength :: HasLength e => Tree e a -> Length +totalBranchLength = bifoldl' (+) const 0 . first getLen + +-- | Normalize branch lengths so that the sum is 1.0. +normalizeBranchLengths :: HasLength e => Tree e a -> Tree e a +normalizeBranchLengths t = first (modLen (/ s)) t + where + s = totalBranchLength t + +-- | Normalize height of tree to 1.0. +normalizeHeight :: HasLength e => Tree e a -> Tree e a +normalizeHeight t = first (modLen (/ h)) t + where + h = height t + +eps :: Double +eps = 1e-12 + +allNearlyEqual :: [Length] -> Bool +allNearlyEqual [] = True +allNearlyEqual xs = all (\y -> eps > abs (fromLength $ x - y)) (tail xs) + where + x = head xs + +-- | Check if a tree is ultrametric. +ultrametric :: HasLength e => Tree e a -> Bool +ultrametric = allNearlyEqual . distancesOriginLeaves + +-- | Elongate terminal branches such that the tree becomes ultrametric. +makeUltrametric :: HasLength e => Tree e a -> Tree e a +makeUltrametric t = go 0 t + where + h = height t + go :: HasLength e => Length -> Tree e a -> Tree e a + go h' (Node br lb []) = let dh = h - h' - getLen br in Node (modLen (+ dh) br) lb [] + go h' (Node br lb ts) = let h'' = h' + getLen br in Node br lb $ map (go h'') ts diff --git a/src/ELynx/Tree/Measurable.hs b/src/ELynx/Tree/Measurable.hs deleted file mode 100644 --- a/src/ELynx/Tree/Measurable.hs +++ /dev/null @@ -1,173 +0,0 @@ -{-# LANGUAGE DeriveAnyClass #-} -{-# LANGUAGE DeriveGeneric #-} -{-# LANGUAGE DerivingVia #-} - --- | --- Module : ELynx.Tree.Measurable --- Description : Measurable 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 lengths is not completely ensured. See the documentation of --- 'Length'. -module ELynx.Tree.Measurable - ( -- * Non-negative length - Length (fromLength), - toLength, - toLengthUnsafe, - checkLength, - Measurable (..), - height, - rootHeight, - - -- * Functions on trees - distancesOriginLeaves, - totalBranchLength, - normalizeBranchLengths, - normalizeHeight, - ultrametric, - makeUltrametric, - ) -where - -import Control.DeepSeq -import Data.Aeson -import Data.Bifoldable -import Data.Bifunctor -import Data.Semigroup -import ELynx.Tree.Rooted -import ELynx.Tree.Splittable -import GHC.Generics - --- | Non-negative length. --- --- However, non-negativity is only checked with 'toLength', and negative values --- can be obtained using the 'Num' and related instances. --- --- Safe operations with conversion from and to length are roughly 50 percent --- slower. --- --- @ --- benchmarking length/length sum foldl' with safe conversion --- time 110.4 ms (109.8 ms .. 111.0 ms) --- 1.000 R² (1.000 R² .. 1.000 R²) --- mean 110.2 ms (110.0 ms .. 110.6 ms) --- std dev 501.8 μs (359.1 μs .. 730.0 μs) --- --- benchmarking length/length sum foldl' num instance --- time 89.37 ms (85.13 ms .. 94.27 ms) --- 0.996 R² (0.992 R² .. 1.000 R²) --- mean 86.53 ms (85.63 ms .. 88.52 ms) --- std dev 2.239 ms (1.069 ms .. 3.421 ms) --- --- benchmarking length/double sum foldl' --- time 85.47 ms (84.88 ms .. 86.42 ms) --- 1.000 R² (0.999 R² .. 1.000 R²) --- mean 85.56 ms (85.26 ms .. 86.02 ms) --- std dev 611.9 μs (101.5 μs .. 851.7 μs) --- @ -newtype Length = Length {fromLength :: Double} - deriving (Read, Show, Generic, NFData) - deriving (Enum, Eq, Floating, Fractional, Num, Ord, Real, RealFloat, RealFrac) via Double - deriving (Semigroup, Monoid) via Sum Double - -instance Splittable Length where - split = (/ 2.0) - -instance ToJSON Length - -instance FromJSON Length - -instance Measurable Length where - getLen = id - setLen = const - modLen f = f - --- | 'Nothing' if length is negative. -toLength :: Double -> Either String Length -toLength x - | x < 0 = Left $ "length: Length is negative: " ++ show x ++ "." - | otherwise = Right $ Length x - --- | Do not check if support value is negative. -toLengthUnsafe :: Double -> Length -toLengthUnsafe = Length - --- | 'Nothing' if length is negative. -checkLength :: Length -> Either String Length -checkLength = toLength . fromLength - --- | A data type with measurable and modifiable length. -class Measurable e where - -- | Get length. - getLen :: e -> Length - - -- | Set length. - setLen :: Length -> e -> e - - -- For computational efficiency. - - -- | Modify length. - modLen :: (Length -> Length) -> e -> e - --- | The maximum distance between origin and leaves. --- --- The height includes the branch length of the stem. -height :: Measurable e => Tree e a -> Length -height = maximum . distancesOriginLeaves - --- | The maximum distance between root node and leaves. -rootHeight :: Measurable e => Tree e a -> Length -rootHeight (Node _ _ []) = 0 -rootHeight t = maximum $ concatMap distancesOriginLeaves (forest t) - --- | Distances from the origin of a tree to the leaves. --- --- The distances include the branch length of the stem. -distancesOriginLeaves :: Measurable e => Tree e a -> [Length] -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 -> Length -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 (modLen (/ 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 (modLen (/ h)) t - where - h = height t - -eps :: Double -eps = 1e-12 - -allNearlyEqual :: [Length] -> Bool -allNearlyEqual [] = True -allNearlyEqual xs = all (\y -> eps > abs (fromLength $ 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 => Length -> Tree e a -> Tree e a - go h' (Node br lb []) = let dh = h - h' - getLen br in Node (modLen (+ dh) br) lb [] - go h' (Node br lb ts) = let h'' = h' + getLen br in Node br lb $ map (go h'') ts diff --git a/src/ELynx/Tree/Name.hs b/src/ELynx/Tree/Name.hs new file mode 100644 --- /dev/null +++ b/src/ELynx/Tree/Name.hs @@ -0,0 +1,72 @@ +{-# LANGUAGE DerivingVia #-} +{-# LANGUAGE TypeSynonymInstances #-} + +-- | +-- Module : ELynx.Tree.Name +-- 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.Name + ( Name (..), + HasName (..), + ) +where + +import Control.DeepSeq +import Data.Aeson +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 +import Data.String + +-- | Node name. +-- +-- Use lazy byte strings because Newick strings are built using chunks. +newtype Name = Name {fromName :: BL.ByteString} + deriving (Show, Read, Eq) + deriving (Ord, Monoid, Semigroup, IsString, NFData) via BL.ByteString + +-- XXX: This is pretty lame, but I need those instances. At the moment, I just +-- go via 'String', but this is certainly not the best solution. + +instance ToJSON Name where + toJSON = toJSON . BL.unpack . fromName + toEncoding = toEncoding . BL.unpack . fromName + +instance FromJSON Name where + parseJSON = fmap (Name . BL.pack) . parseJSON + +-- | Class of types having a name. +class HasName a where + getName :: a -> Name + +instance HasName Name where + getName = id + +instance HasName () where + getName = const (Name BL.empty) + +instance HasName Int where + getName = Name . BB.toLazyByteString . BB.intDec + +instance HasName Double where + getName = Name . BL.fromStrict . toShortest + +instance HasName Char where + getName = Name . BB.toLazyByteString . BB.char8 + +instance (HasName a) => HasName [a] where + getName = Name . BL.concat . map (fromName . getName) + +instance HasName BL.ByteString where + getName = Name + +instance HasName BS.ByteString where + getName = Name . BL.fromStrict diff --git a/src/ELynx/Tree/Named.hs b/src/ELynx/Tree/Named.hs deleted file mode 100644 --- a/src/ELynx/Tree/Named.hs +++ /dev/null @@ -1,72 +0,0 @@ -{-# LANGUAGE DerivingVia #-} -{-# LANGUAGE TypeSynonymInstances #-} - --- | --- 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 - ( Name (..), - Named (..), - ) -where - -import Control.DeepSeq -import Data.Aeson -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 -import Data.String - --- | Node name. --- --- Use lazy byte strings because Newick strings are built using chunks. -newtype Name = Name {fromName :: BL.ByteString} - deriving (Show, Eq) - deriving (Ord, Monoid, Semigroup, IsString, NFData) via BL.ByteString - --- XXX: This is pretty lame, but I need those instances. At the moment, I just --- go via 'String', but this is certainly not the best solution. - -instance ToJSON Name where - toJSON = toJSON . BL.unpack . fromName - toEncoding = toEncoding . BL.unpack . fromName - -instance FromJSON Name where - parseJSON = fmap (Name . BL.pack) . parseJSON - -instance Named Name where - getName = id - --- | Data types with names. -class Named a where - getName :: a -> Name - -instance Named () where - getName = const (Name BL.empty) - -instance Named Int where - getName = Name . BB.toLazyByteString . BB.intDec - -instance Named Double where - getName = Name . BL.fromStrict . toShortest - -instance Named Char where - getName = Name . BB.toLazyByteString . BB.char8 - -instance (Named a) => Named [a] where - getName = Name . BL.concat . map (fromName . getName) - -instance Named BL.ByteString where - getName = Name - -instance Named BS.ByteString where - getName = Name . BL.fromStrict diff --git a/src/ELynx/Tree/Phylogeny.hs b/src/ELynx/Tree/Phylogeny.hs --- a/src/ELynx/Tree/Phylogeny.hs +++ b/src/ELynx/Tree/Phylogeny.hs @@ -39,23 +39,26 @@ -- -- 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. +-- unrooted, when the root node is multifurcating and has three or more +-- children. This convention is not used here. Newick trees are just parsed as +-- they are, and a rooted tree is returned. -- --- Trees with multifurcating root nodes can be rooted using 'outgroup'. +-- A multifurcating root node can be resolved to a bifurcating root node with +-- 'outgroup'. -- --- Trees with bifurcating root nodes can be changed with 'rootAt' or 'midpoint'; --- a list of all rooted trees is returned by 'roots'. +-- The bifurcating root node can be changed with 'outgroup' or 'midpoint'. +-- +-- For a given tree with bifurcating root node, a list of all rooted trees is +-- returned by 'roots'. module ELynx.Tree.Phylogeny ( -- * Functions equal, + equal', intersect, bifurcating, outgroup, midpoint, roots, - rootAt, -- * Branch labels Phylo (..), @@ -82,22 +85,33 @@ import Data.Set (Set) import qualified Data.Set as S import ELynx.Tree.Bipartition -import ELynx.Tree.Measurable +import ELynx.Tree.Length import ELynx.Tree.Rooted import ELynx.Tree.Splittable -import ELynx.Tree.Supported +import ELynx.Tree.Support import GHC.Generics --- A faster check could probably be done using 'Ord' and sets. The leave set --- could be precomputed. -- | 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) = +-- +-- NOTE: The equality check is only meaningful if the trees have unique leaves. +-- +-- Return 'Left' if a tree does not have unique leaves. +equal :: (Eq e, Eq a, Ord a) => Tree e a -> Tree e a -> Either String Bool +equal tL tR + | duplicateLeaves tL = Left "equal: Left tree has duplicate leaves." + | duplicateLeaves tR = Left "equal: Right tree has duplicate leaves." + | otherwise = Right $ equal' tL tR + +-- | Same as 'equal', but assume that leaves are unique. +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 + && all (elem' tsR) tsL + where + elem' ts t = isJust $ find (equal' t) ts -- | Compute the intersection of trees. -- @@ -106,6 +120,7 @@ -- 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) @@ -148,33 +163,40 @@ -- 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. +-- | Root the tree using an outgroup. -- --- A bifurcating root node with the provided label is introduced. The affected --- branch is 'split'. +-- If the current root node is multifurcating, a bifurcating root node with the +-- empty label is introduced by 'split'ting the leftmost branch. The 'Monoid' +-- instance of the node label and the 'Splittable' instance of the branch length +-- are used. Note that in this case, the degree of the former root node is +-- decreased by one! -- --- Note, the degree of the former root node is decreased by one. +-- Given that the root note is bifurcating, the root node is moved to the +-- required position specified by the outgroup. -- --- If the root node is bifurcating, use 'rootAt'. +-- Branches are connected according to the provided 'Semigroup' instance. -- +-- Upon insertion of the root node at the required position, the affected branch +-- is 'split' according to the provided 'Splittable' instance. +-- -- Return 'Left' if +-- -- - the root node is not multifurcating; +-- -- - the tree has duplicate leaves; +-- -- - 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' +outgroup :: (Semigroup e, Splittable e, Monoid a, Ord a) => Set 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 o t@(Node _ _ [_, _]) = do + bip <- bp o (S.fromList (leaves t) S.\\ o) + rootAt bip t +outgroup o (Node b l ts) = outgroup o 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)] + t' = Node b mempty [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 @@ -185,8 +207,9 @@ -- | 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 :: (Semigroup e, Splittable e, HasLength 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 @@ -199,7 +222,7 @@ 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 :: HasLength e => [Tree e a] -> Tree e a getMidpoint ts = case t of (Node br lb [l, r]) -> let hl = height l @@ -222,7 +245,7 @@ -- 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 -> Length +getDeltaHeight :: HasLength e => Tree e a -> Length getDeltaHeight (Node _ _ [l, r]) = abs $ height l - height r -- Explicitly use 'error' here, because roots is supposed to return trees with -- bifurcating root nodes. @@ -238,14 +261,14 @@ -- 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 +-- A bifurcating root is required because moving a multifurcating root node to +-- another branch would change the degree of the root node. 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. +-- Split the affected branch into one out of two equal entities according the +-- provided 'Splittable' instance. -- -- Return 'Left' if the root node is not 'bifurcating'. roots :: (Semigroup e, Splittable e) => Tree e a -> Either String (Forest e a) @@ -277,7 +300,7 @@ tC' = tC {branch = brC'} cfs = complementaryForests tC' tsD --- | Root a tree at a specific position. +-- 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. @@ -290,8 +313,11 @@ -- 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) => @@ -318,7 +344,7 @@ Either String (Tree e a) rootAt' b t = do ts <- roots t - case find (\x -> Right b == bipartition x) ts of + case find (\x -> bipartition x == Right b) ts of Nothing -> Left "rootAt': Bipartition not found on tree." Just t' -> Right t' @@ -344,10 +370,10 @@ -- | Set all branch lengths 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 :: (HasLength e, HasSupport e) => Tree e a -> Tree Phylo a toPhyloTree = first toPhyloLabel -toPhyloLabel :: (Measurable e, Supported e) => e -> Phylo +toPhyloLabel :: (HasLength e, HasSupport e) => e -> Phylo toPhyloLabel x = Phylo (Just $ getLen x) (Just $ getSup x) -- | Set all branch lengths to 'Just' the values, and all support values to @@ -355,10 +381,10 @@ -- -- Useful to export a tree with branch lengths but without branch support values -- to Newick format. -measurableToPhyloTree :: Measurable e => Tree e a -> Tree Phylo a +measurableToPhyloTree :: HasLength e => Tree e a -> Tree Phylo a measurableToPhyloTree = first measurableToPhyloLabel -measurableToPhyloLabel :: Measurable e => e -> Phylo +measurableToPhyloLabel :: HasLength e => e -> Phylo measurableToPhyloLabel x = Phylo (Just $ getLen x) Nothing -- | Set all branch lengths to 'Nothing', and all support values to 'Just' the @@ -366,10 +392,10 @@ -- -- Useful to export a tree with branch support values but without branch lengths -- to Newick format. -supportedToPhyloTree :: Supported e => Tree e a -> Tree Phylo a +supportedToPhyloTree :: HasSupport e => Tree e a -> Tree Phylo a supportedToPhyloTree = first supportedToPhyloLabel -supportedToPhyloLabel :: Supported e => e -> Phylo +supportedToPhyloLabel :: HasSupport e => e -> Phylo supportedToPhyloLabel x = Phylo Nothing (Just $ getSup x) -- | If root branch length is not available, set it to 0. @@ -431,7 +457,7 @@ instance Semigroup PhyloExplicit where PhyloExplicit bL sL <> PhyloExplicit bR sR = PhyloExplicit (bL + bR) (min sL sR) -instance Measurable PhyloExplicit where +instance HasLength PhyloExplicit where getLen = sBrLen setLen b pl = pl {sBrLen = b} modLen f (PhyloExplicit l s) = PhyloExplicit (f l) s @@ -441,7 +467,7 @@ where b' = sBrLen l / 2.0 -instance Supported PhyloExplicit where +instance HasSupport PhyloExplicit where getSup = sBrSup setSup s pl = pl {sBrSup = s} modSup f (PhyloExplicit l s) = PhyloExplicit l (f s) diff --git a/src/ELynx/Tree/Rooted.hs b/src/ELynx/Tree/Rooted.hs --- a/src/ELynx/Tree/Rooted.hs +++ b/src/ELynx/Tree/Rooted.hs @@ -41,6 +41,21 @@ -- -- In mathematical terms: A 'Tree' is a directed acyclic graph without loops, -- with vertex labels, and with edge labels. +-- +-- A short recap of recursive tree traversals: +-- +-- - Pre-order: Root first, then sub trees from left to right. Also called depth +-- first. +-- +-- - In-order: Only valid for bifurcating trees. Left sub tree first, then root, +-- then right sub tree. +-- +-- - Post-order: Sub trees from left to right, then the root. Also called +-- breadth first. +-- +-- Here, pre-order traversals are used exclusively, for example, by accessor +-- functions such as 'branches', or 'labels' which is the same as 'toList'. +-- Please let me know, if post-order algorithms are required. module ELynx.Tree.Rooted ( -- * Data type Tree (..), @@ -168,7 +183,7 @@ -- -- The 'Monoid' instance of the branch labels determines the default branch -- label when using 'pure'. -instance Monoid e => Applicative (Tree e) where +instance (Semigroup e, 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) @@ -182,7 +197,10 @@ -- | 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 +-- +-- The 'Monoid' instance of the branch labels determines the default branch +-- label when using 'return'. +instance (Semigroup e, 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) @@ -231,11 +249,18 @@ toTreeNodeLabels :: Tree e a -> T.Tree a toTreeNodeLabels (Node _ lb ts) = T.Node lb (map toTreeNodeLabels ts) +-- TODO: Maybe use foldr similar to 'labels'. -- | Get leaves. leaves :: Tree e a -> [a] leaves (Node _ lb []) = [lb] leaves (Node _ _ ts) = concatMap leaves ts +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 tree has duplicate leaves. duplicateLeaves :: Ord a => Tree e a -> Bool duplicateLeaves = duplicates . leaves @@ -364,9 +389,3 @@ -- 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 diff --git a/src/ELynx/Tree/Simulate/Coalescent.hs b/src/ELynx/Tree/Simulate/Coalescent.hs --- a/src/ELynx/Tree/Simulate/Coalescent.hs +++ b/src/ELynx/Tree/Simulate/Coalescent.hs @@ -16,7 +16,7 @@ import Control.Monad.Primitive import ELynx.Tree.Distribution.CoalescentContinuous -import ELynx.Tree.Measurable +import ELynx.Tree.Length import ELynx.Tree.Rooted import Statistics.Distribution import System.Random.MWC diff --git a/src/ELynx/Tree/Simulate/PointProcess.hs b/src/ELynx/Tree/Simulate/PointProcess.hs --- a/src/ELynx/Tree/Simulate/PointProcess.hs +++ b/src/ELynx/Tree/Simulate/PointProcess.hs @@ -40,7 +40,7 @@ import ELynx.Tree.Distribution.TimeOfOrigin import ELynx.Tree.Distribution.TimeOfOriginNearCritical import ELynx.Tree.Distribution.Types -import ELynx.Tree.Measurable +import ELynx.Tree.Length import ELynx.Tree.Rooted import qualified Statistics.Distribution as D ( genContVar, diff --git a/src/ELynx/Tree/Support.hs b/src/ELynx/Tree/Support.hs new file mode 100644 --- /dev/null +++ b/src/ELynx/Tree/Support.hs @@ -0,0 +1,105 @@ +{-# LANGUAGE DeriveAnyClass #-} +{-# LANGUAGE DeriveGeneric #-} +{-# LANGUAGE DerivingVia #-} + +-- | +-- Module : ELynx.Tree.Support +-- Description : Labels with support values +-- 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. +module ELynx.Tree.Support + ( -- * Non-negative support value + Support (fromSupport), + toSupport, + toSupportUnsafe, + HasSupport (..), + + -- * Functions on trees + normalizeBranchSupport, + collapse, + ) +where + +import Control.DeepSeq +import Data.Aeson +import Data.Bifoldable +import Data.Bifunctor +import Data.List +import Data.Semigroup +import ELynx.Tree.Rooted +import ELynx.Tree.Splittable +import GHC.Generics + +-- | Non-negative support value. +-- +-- However, non-negativity is only checked with 'toSupport', and negative values +-- can be obtained using the 'Num' and related instances. +-- +-- See also the documentation of 'ELynx.Tree.Length.Length'. +newtype Support = Support {fromSupport :: Double} + deriving (Read, Show, Generic, NFData) + deriving (Enum, Eq, Floating, Fractional, Num, Ord, Real, RealFloat, RealFrac) via Double + deriving (Semigroup) via Min Double + +instance Splittable Support where + split = id + +instance ToJSON Support + +instance FromJSON Support + +instance HasSupport Support where + getSup = id + setSup = const + modSup f = f + +-- | If negative, call 'error' indicating the calling function name. +toSupport :: String -> Double -> Support +toSupport s x + | x < 0 = error $ s ++ ": Support is negative: " ++ show x ++ "." + | otherwise = Support x + +-- | Do not check if value is negative. +toSupportUnsafe :: Double -> Support +toSupportUnsafe = Support + +-- | A data type with measurable and modifiable values. +class HasSupport e where + getSup :: e -> Support + setSup :: Support -> e -> e + + -- For computational efficiency. + modSup :: (Support -> Support) -> e -> e + +-- | Normalize branch support values. The maximum branch support value will be +-- set to 1.0. +normalizeBranchSupport :: HasSupport e => Tree e a -> Tree e a +normalizeBranchSupport t = first (modSup (/ 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, HasSupport e) => Support -> 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 :: HasSupport e => Support -> Tree e a -> Bool +highP _ (Node _ _ []) = True +highP th (Node br _ _) = getSup br >= th + +-- See 'collapse'. +collapse' :: HasSupport e => Support -> 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 diff --git a/src/ELynx/Tree/Supported.hs b/src/ELynx/Tree/Supported.hs deleted file mode 100644 --- a/src/ELynx/Tree/Supported.hs +++ /dev/null @@ -1,106 +0,0 @@ -{-# LANGUAGE DeriveAnyClass #-} -{-# LANGUAGE DeriveGeneric #-} -{-# LANGUAGE DerivingVia #-} - --- | --- Module : ELynx.Tree.Supported --- Description : Labels with support values --- 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. -module ELynx.Tree.Supported - ( -- * Non-negative support value - Support (fromSupport), - toSupport, - toSupportUnsafe, - Supported (..), - - -- * Functions on trees - normalizeBranchSupport, - collapse, - ) -where - -import Control.DeepSeq -import Data.Aeson -import Data.Bifoldable -import Data.Bifunctor -import Data.List -import Data.Semigroup -import ELynx.Tree.Rooted -import ELynx.Tree.Splittable -import GHC.Generics - --- | Non-negative support value. --- --- However, non-negativity is only checked with 'toSupport', and negative values --- can be obtained using the 'Num' and related instances. --- --- See also the documentation of 'ELynx.Tree.Measurable.Length'. -newtype Support = Support {fromSupport :: Double} - deriving (Read, Show, Generic, NFData) - deriving (Enum, Eq, Floating, Fractional, Num, Ord, Real, RealFloat, RealFrac) via Double - deriving (Semigroup) via Min Double - -instance Splittable Support where - split = id - -instance ToJSON Support - -instance FromJSON Support - -instance Supported Support where - getSup = id - setSup = const - modSup f = f - --- | Nothing if support is negative. -toSupport :: Double -> Either String Support -toSupport x - | x < 0 = Left $ "toSupport: Support is negative: " ++ show x ++ "." - | otherwise = Right $ Support x - --- | Do not check if support value is negative. -toSupportUnsafe :: Double -> Support -toSupportUnsafe = Support - --- | A data type that supports extraction, setting and modifying of support --- values. -class Supported e where - getSup :: e -> Support - setSup :: Support -> e -> e - - -- For computational efficiency. - modSup :: (Support -> Support) -> e -> e - --- | 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 (modSup (/ 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) => Support -> 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 => Support -> Tree e a -> Bool -highP _ (Node _ _ []) = True -highP th (Node br _ _) = getSup br >= th - --- See 'collapse'. -collapse' :: Supported e => Support -> 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 diff --git a/src/ELynx/Tree/Zipper.hs b/src/ELynx/Tree/Zipper.hs --- a/src/ELynx/Tree/Zipper.hs +++ b/src/ELynx/Tree/Zipper.hs @@ -18,18 +18,24 @@ toTree, -- * Movement - goUp, + goParent, + goParentUnsafe, goRoot, goLeft, goRight, goChild, + goChildUnsafe, + + -- * Paths Path, goPath, + validPath, goPathUnsafe, getSubTreeUnsafe, -- * Modification insertTree, + modifyTree, insertBranch, insertLabel, ) @@ -63,8 +69,8 @@ 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 +goParent :: TreePos e a -> Maybe (TreePos e a) +goParent pos = case parents pos of (ls, br, lb, rs) : ps -> Just Pos @@ -75,9 +81,23 @@ } [] -> Nothing +-- | Go to parent. +-- +-- Call 'error' if no parent is found. +goParentUnsafe :: TreePos e a -> TreePos e a +goParentUnsafe pos = case parents pos of + (ls, br, lb, rs) : ps -> + Pos + { current = Node br lb $ getForest pos, + before = ls, + after = rs, + parents = ps + } + [] -> error "goUpUnsafe: No parent found." + -- | Go to root. goRoot :: TreePos e a -> TreePos e a -goRoot pos = maybe pos goRoot (goUp pos) +goRoot pos = maybe pos goRoot (goParent pos) -- | Go to left sibling in current forest. goLeft :: TreePos e a -> Maybe (TreePos e a) @@ -122,16 +142,6 @@ where (ls', rs') = splitAt n ts --- | Path from the root of a tree to the node of the tree. --- --- The position is specific to a tree topology. If the topology changes, the --- position becomes invalid. -type Path = [Int] - --- | Go to node with given path. -goPath :: Path -> 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. goChildUnsafe :: Int -> TreePos e a -> TreePos e a @@ -149,11 +159,29 @@ where (ls', rs') = splitAt n ts +-- | Path from the root of a tree to the node of the tree. +-- +-- The position is specific to a tree topology. If the topology changes, the +-- position becomes invalid. +type Path = [Int] + +-- | Go to node with given path. +goPath :: Path -> TreePos e a -> Maybe (TreePos e a) +goPath pos pth = foldlM (flip goChild) pth pos + +-- | Check if a path is valid in that it leads to a node on a tree. +validPath :: Tree e a -> Path -> Either String Path +validPath t p = case goPath p (fromTree t) of + Nothing -> Left "validPath: Path is invalid." + Just _ -> Right p + -- | Got to node with given path. -- -- Call 'error' if path is invalid. goPathUnsafe :: Path -> TreePos e a -> TreePos e a -goPathUnsafe pos pth = foldl (flip goChildUnsafe) pth pos +goPathUnsafe pos pth = + {-# SCC "goPathUnsafe" #-} + foldl (flip goChildUnsafe) pth pos -- | Get the sub tree at path. -- @@ -164,6 +192,12 @@ -- | 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} + +-- | Modify the tree at the current focus of the zipper. +modifyTree :: (Tree e a -> Tree e a) -> TreePos e a -> TreePos e a +modifyTree f pos = pos {current = f t} + where + t = current pos -- | Insert a new branch label into the current focus of the zipper. insertBranch :: e -> TreePos e a -> TreePos e a diff --git a/test/ELynx/Tree/PhylogenySpec.hs b/test/ELynx/Tree/PhylogenySpec.hs --- a/test/ELynx/Tree/PhylogenySpec.hs +++ b/test/ELynx/Tree/PhylogenySpec.hs @@ -25,9 +25,24 @@ import Test.Hspec.QuickCheck import Test.QuickCheck hiding (labels) -simpleTree :: Tree () String -simpleTree = Node () "i" [Node () "j" [Node () "x" [], Node () "y" []], Node () "z" []] +simpleTree1 :: Tree () String +simpleTree1 = Node () "i" [Node () "j" [Node () "x" [], Node () "y" []], Node () "z" []] +simpleTree2 :: Tree () String +simpleTree2 = Node () "i" [Node () "j" [Node () "y" [], Node () "x" []], Node () "z" []] + +simpleTree3 :: Tree () String +simpleTree3 = Node () "i" [Node () "j" [Node () "x" [], Node () "z" []], Node () "y" []] + +prop_commutative :: (Eq a, Ord a) => Tree () a -> Tree () a -> Bool +prop_commutative t1 t2 = case (s1, s2) of + (Left _, Left _) -> True + (Right x1, Right x2) -> x1 == x2 + _ -> False + where + s1 = t1 `equal` t2 + s2 = t2 `equal` t1 + simpleSol :: Forest () String simpleSol = [ Node () "i" [Node () "j" [Node () "x" [], Node () "y" []], Node () "z" []], @@ -91,23 +106,27 @@ -- -- TODO. -- prop_bifurcating_tree --- :: (Ord a, Measurable a, Named a, Supported a) => Tree a -> Bool +-- :: (Ord a, HasLength a, HasName a, HasSupport a) => Tree a -> Bool -- prop_bifurcating_tree t = partitions (resolve t) == empty prop_roots_total_length :: Tree Length a -> Bool prop_roots_total_length t@(Node _ _ [_, _]) = - all (\b -> abs (b - l) < 1e-8) $ - map totalBranchLength $ - either error id $ - roots t + all (\x -> abs (totalBranchLength x - l) < 1e-8) $ + either error id $ + roots t where l = totalBranchLength t prop_roots_total_length _ = True spec :: Spec spec = do - -- TODO: describe "Resolve" - + describe "equal" $ do + it "correctly handles some test cases" $ do + simpleTree1 `equal` simpleTree2 `shouldBe` Right True + simpleTree1 `equal` simpleTree3 `shouldBe` Right False + simpleTree2 `equal` simpleTree3 `shouldBe` Right False + it "is commutative" $ + property (prop_commutative :: Tree () Int -> Tree () Int -> Bool) describe "roots" $ do it "correctly handles leaves and cherries" $ do let tleaf = Node () 0 [] :: Tree () Int @@ -115,20 +134,19 @@ roots tleaf `shouldSatisfy` isLeft roots tcherry `shouldBe` Right [tcherry] it "correctly handles simple trees" $ - either error id (roots simpleTree) `shouldBe` simpleSol + either error id (roots simpleTree1) `shouldBe` simpleSol modifyMaxSize (* 100) $ it "returns the correct number of rooted trees for arbitrary trees" $ property (prop_roots :: (Tree () Int -> Bool)) - describe "rootAt" $ + describe "outgroup" $ modifyMaxSize (* 100) $ it "correctly handles simple trees" $ do - let p = either error id $ bipartition simpleTree - rootAt p simpleTree `shouldBe` Right simpleTree + let p = fst $ fromBipartition $ either error id $ bipartition simpleTree1 + outgroup p simpleTree1 `shouldBe` Right simpleTree1 let l = S.singleton "x" - r = S.fromList ["y", "z"] - p' = either error id $ bp l r - either error id (rootAt p' simpleTree) `shouldSatisfy` (`equal` (simpleSol !! 1)) + either error id (outgroup l simpleTree1) `equal` (simpleSol !! 1) + `shouldBe` Right True describe "rootsWithBranch" $ modifyMaxSize (* 100) $ it "does not change the tree height" $ diff --git a/test/ELynx/Tree/RootedSpec.hs b/test/ELynx/Tree/RootedSpec.hs --- a/test/ELynx/Tree/RootedSpec.hs +++ b/test/ELynx/Tree/RootedSpec.hs @@ -50,7 +50,7 @@ prune largeTree `shouldBe` largeTree it "correctly prunes a small example" $ prune smallSubTree `shouldBe` smallSubTreePruned - it "leaves height constant for Measurable trees" $ do + it "leaves height constant for trees with branch lengths" $ do let t' = either error id $ phyloToLengthTree subSampleLargeTree @@ -62,5 +62,3 @@ dropLeavesWith (const True) smallTree `shouldBe` Nothing it "returns the correct subtree for a small example" $ dropLeavesWith (== 2) smallTree `shouldBe` Just smallSubTree - --- TODO: intersect. diff --git a/test/ELynx/Tree/SupportSpec.hs b/test/ELynx/Tree/SupportSpec.hs new file mode 100644 --- /dev/null +++ b/test/ELynx/Tree/SupportSpec.hs @@ -0,0 +1,39 @@ +{-# LANGUAGE OverloadedStrings #-} + +-- | +-- Module : ELynx.Tree.SupportSpec +-- Description : Unit tests for ELynx.Tree.SupportSpec +-- Copyright : (c) Dominik Schrempf, 2020 +-- License : GPL-3.0-or-later +-- +-- Maintainer : dominik.schrempf@gmail.com +-- Stability : unstable +-- Portability : portable +-- +-- Creation date: Fri Aug 21 14:20:09 2020. +module ELynx.Tree.SupportSpec + ( spec, + ) +where + +import ELynx.Tools +import ELynx.Tree +import Test.Hspec + +collapseTree :: Tree Phylo Name +collapseTree = parseByteStringWith (oneNewick IqTree) "((a,b),(c,d));" + +collapseStarTree :: Tree Phylo Name +collapseStarTree = parseByteStringWith (oneNewick Standard) "(a[1.0],b[1.0],c[1.0],d[1.0])[1.0];" + +spec :: Spec +spec = do + describe "collapse" $ do + it "creates a star tree for 1.0" $ do + let t = phyloToSupportTreeUnsafe collapseTree + s = phyloToSupportTreeUnsafe collapseStarTree + collapse 0 t `shouldBe` t + collapse 0.01 t `shouldBe` t + collapse 0.99 t `shouldBe` t + collapse 1.0 t `shouldBe` t + collapse 1.1 t `shouldBe` s diff --git a/test/ELynx/Tree/SupportedSpec.hs b/test/ELynx/Tree/SupportedSpec.hs deleted file mode 100644 --- a/test/ELynx/Tree/SupportedSpec.hs +++ /dev/null @@ -1,39 +0,0 @@ -{-# LANGUAGE OverloadedStrings #-} - --- | --- Module : ELynx.Tree.SupportedSpec --- Description : Unit tests for ELynx.Tree.SupportedSpec --- Copyright : (c) Dominik Schrempf, 2020 --- License : GPL-3.0-or-later --- --- Maintainer : dominik.schrempf@gmail.com --- Stability : unstable --- Portability : portable --- --- Creation date: Fri Aug 21 14:20:09 2020. -module ELynx.Tree.SupportedSpec - ( spec, - ) -where - -import ELynx.Tools -import ELynx.Tree -import Test.Hspec - -collapseTree :: Tree Phylo Name -collapseTree = parseByteStringWith (oneNewick IqTree) "((a,b),(c,d));" - -collapseStarTree :: Tree Phylo Name -collapseStarTree = parseByteStringWith (oneNewick Standard) "(a[1.0],b[1.0],c[1.0],d[1.0])[1.0];" - -spec :: Spec -spec = do - describe "collapse" $ do - it "creates a star tree for 1.0" $ do - let t = phyloToSupportTreeUnsafe collapseTree - s = phyloToSupportTreeUnsafe collapseStarTree - collapse 0 t `shouldBe` t - collapse 0.01 t `shouldBe` t - collapse 0.99 t `shouldBe` t - collapse 1.0 t `shouldBe` t - collapse 1.1 t `shouldBe` s