packages feed

elynx-tree 0.2.2 → 0.3.0

raw patch · 53 files changed

+4003/−3118 lines, 53 filesdep +aesondep +comonaddep +deepseqdep −quickcheck-instancesdep ~QuickCheckdep ~basedep ~bytestring

Dependencies added: aeson, comonad, deepseq, elynx-nexus, vector

Dependencies removed: quickcheck-instances

Dependency ranges changed: QuickCheck, base, bytestring, containers, criterion, elynx-tools, hspec, hspec-megaparsec, math-functions, megaparsec, mwc-random, primitive, statistics

Files

ChangeLog.md view
@@ -2,6 +2,17 @@ # Changelog for ELynx  +## Version 0.3.1; in development+++## Version 0.3.0++-   **`elynx-nexus`:** library to import and export Nexus files.+-   **`elynx-tree`:** major refactor and big cleanup; use rose trees with branch+    labels.+-   **`elynx-tree`:** provide zippers.++ ## Version 0.2.2  -   Validation and repetition of previous analyses is finally possible with the
README.md view
@@ -2,9 +2,11 @@  # The ELynx Suite -Version: 0.2.2.+Version: 0.3.0. Reproducible evolution made easy. +<p align="center"><img src="https://travis-ci.org/dschrempf/elynx.svg?branch=master"/></p>+ A Haskell library and tool set for computational biology. The goal of ELynx is reproducible research. Evolutionary sequences and phylogenetic trees can be read, viewed, modified and simulated. The command line with all arguments is@@ -12,11 +14,12 @@ sums so that validation of past analyses is possible without the need to recompute the result. -The Elynx Suite consists of four library packages and three executables providing-a range of sub commands.+The Elynx Suite consists of library packages and executables providing a range+of sub commands.  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.@@ -43,7 +46,7 @@  2.  Clone the ELynx repository.     -        git clone clone https://github.com/dschrempf/elynx+        git clone https://github.com/dschrempf/elynx  3.  Navigate to the newly created `elynx` folder and build the binaries.     This will take a while.@@ -62,15 +65,34 @@     [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+    slynx --help | head -n -16 -    ELynx Suite version 0.2.1.+    ELynx Suite version 0.3.0.     Developed by Dominik Schrempf.-    Compiled on April 27, 2020, at 12:08 pm, UTC.+    Compiled on July 30, 2020, at 12:40 pm, UTC.          Usage: slynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME]                   [-f|--force] COMMAND@@ -95,6 +117,7 @@       sub-sample               Sub-sample columns from multi sequence alignments.       translate                Translate from DNA to Protein or DNAX to ProteinX.     +         Available sequence file formats:       - FASTA     @@ -105,23 +128,6 @@       - Protein (amino acids)       - ProteinX (amino acids; including gaps)       - ProteinS (amino acids; including gaps, and translation stops)-      - ProteinI (amino acids; including gaps, translation stops, and IUPAC codes)-    -    ELynx-    ------    A Haskell library and tool set for computational biology. The goal of ELynx is-    reproducible research. Evolutionary sequences and phylogenetic trees can be-    read, viewed, modified and simulated. The command line with all arguments is-    logged consistently, and automatically. Data integrity is verified using SHA256-    sums so that validation of past analyses is possible without the need to-    recompute the result.-    -    slynx     Analyze, modify, and simulate evolutionary sequences.-    tlynx     Analyze, modify, and simulate phylogenetic trees.-    elynx     Validate and redo past analyses.-    -    Get help for sub commands:-      slynx examine --help   ## Concatenate@@ -130,9 +136,9 @@      slynx concatenate --help -    ELynx Suite version 0.2.1.+    ELynx Suite version 0.3.0.     Developed by Dominik Schrempf.-    Compiled on April 27, 2020, at 12:08 pm, UTC.+    Compiled on July 30, 2020, at 12:40 pm, UTC.          Usage: slynx concatenate (-a|--alphabet NAME) INPUT-FILE       Concatenate sequences found in input files.@@ -143,23 +149,6 @@       -a,--alphabet NAME       Specify alphabet type NAME       INPUT-FILE               Read sequences from INPUT-FILE       -h,--help                Show this help text-    -    -    ELynx-    ------    A Haskell library and tool set for computational biology. The goal of ELynx is-    reproducible research. Evolutionary sequences and phylogenetic trees can be-    read, viewed, modified and simulated. The command line with all arguments is-    logged consistently, and automatically. Data integrity is verified using SHA256-    sums so that validation of past analyses is possible without the need to-    recompute the result.-    -    slynx     Analyze, modify, and simulate evolutionary sequences.-    tlynx     Analyze, modify, and simulate phylogenetic trees.-    elynx     Validate and redo past analyses.-    -    Get help for sub commands:-      slynx examine --help   ## Examine@@ -168,9 +157,9 @@      slynx examine --help -    ELynx Suite version 0.2.1.+    ELynx Suite version 0.3.0.     Developed by Dominik Schrempf.-    Compiled on April 27, 2020, at 12:08 pm, UTC.+    Compiled on July 30, 2020, at 12:40 pm, UTC.          Usage: slynx examine (-a|--alphabet NAME) INPUT-FILE [--per-site]       Examine sequences. If data is a multi sequence alignment, additionally analyze columns.@@ -182,23 +171,6 @@       INPUT-FILE               Read sequences from INPUT-FILE       --per-site               Report per site summary statistics       -h,--help                Show this help text-    -    -    ELynx-    ------    A Haskell library and tool set for computational biology. The goal of ELynx is-    reproducible research. Evolutionary sequences and phylogenetic trees can be-    read, viewed, modified and simulated. The command line with all arguments is-    logged consistently, and automatically. Data integrity is verified using SHA256-    sums so that validation of past analyses is possible without the need to-    recompute the result.-    -    slynx     Analyze, modify, and simulate evolutionary sequences.-    tlynx     Analyze, modify, and simulate phylogenetic trees.-    elynx     Validate and redo past analyses.-    -    Get help for sub commands:-      slynx examine --help   ## Filter@@ -207,9 +179,9 @@      slynx filter-rows --help -    ELynx Suite version 0.2.1.+    ELynx Suite version 0.3.0.     Developed by Dominik Schrempf.-    Compiled on April 27, 2020, at 12:08 pm, UTC.+    Compiled on July 30, 2020, at 12:40 pm, UTC.          Usage: slynx filter-rows (-a|--alphabet NAME) INPUT-FILE [--longer-than LENGTH]                               [--shorter-than LENGTH] [--standard-characters]@@ -225,31 +197,14 @@       --standard-characters    Only keep sequences containing at least one standard                                (i.e., non-IUPAC) character       -h,--help                Show this help text-    -    -    ELynx-    ------    A Haskell library and tool set for computational biology. The goal of ELynx is-    reproducible research. Evolutionary sequences and phylogenetic trees can be-    read, viewed, modified and simulated. The command line with all arguments is-    logged consistently, and automatically. Data integrity is verified using SHA256-    sums so that validation of past analyses is possible without the need to-    recompute the result.-    -    slynx     Analyze, modify, and simulate evolutionary sequences.-    tlynx     Analyze, modify, and simulate phylogenetic trees.-    elynx     Validate and redo past analyses.-    -    Get help for sub commands:-      slynx examine --help  Filter columns of multi sequence alignments with `filter-columns`.      slynx filter-columns --help -    ELynx Suite version 0.2.1.+    ELynx Suite version 0.3.0.     Developed by Dominik Schrempf.-    Compiled on April 27, 2020, at 12:08 pm, UTC.+    Compiled on July 30, 2020, at 12:40 pm, UTC.          Usage: slynx filter-columns (-a|--alphabet NAME) INPUT-FILE                                  [--standard-chars DOUBLE]@@ -263,23 +218,6 @@       --standard-chars DOUBLE  Keep columns with a proportion standard (non-IUPAC)                                characters larger than DOUBLE in [0,1]       -h,--help                Show this help text-    -    -    ELynx-    ------    A Haskell library and tool set for computational biology. The goal of ELynx is-    reproducible research. Evolutionary sequences and phylogenetic trees can be-    read, viewed, modified and simulated. The command line with all arguments is-    logged consistently, and automatically. Data integrity is verified using SHA256-    sums so that validation of past analyses is possible without the need to-    recompute the result.-    -    slynx     Analyze, modify, and simulate evolutionary sequences.-    tlynx     Analyze, modify, and simulate phylogenetic trees.-    elynx     Validate and redo past analyses.-    -    Get help for sub commands:-      slynx examine --help   ## Simulate@@ -288,9 +226,9 @@      slynx simulate --help -    ELynx Suite version 0.2.1.+    ELynx Suite version 0.3.0.     Developed by Dominik Schrempf.-    Compiled on April 27, 2020, at 12:08 pm, UTC.+    Compiled on July 30, 2020, at 12:40 pm, UTC.          Usage: slynx simulate (-t|--tree-file Name) [-s|--substitution-model MODEL]                            [-m|--mixture-model MODEL] [-e|--edm-file NAME] @@ -352,27 +290,13 @@        where XX is 10, 20, 30, 40, 50, or 60; CXX models, Quang et al., 2008.     -m "EDM(EXCHANGEABILITIES)"        Arbitrary empirical distribution mixture (EDM) models.-       Stationary distributions have to be provided with the -e option.+       Stationary distributions have to be provided with the -e or -p option.        For example,          LG exchangeabilities with stationary distributions given in FILE.          -m "EDM(LG-Custom)" -e FILE+         LG exchangeabilities with site profiles (Phylobayes) given in FILES.+         -m "EDM(LG-Custom)" -p FILES     For special mixture models, mixture weights are optional.-    -    ELynx-    ------    A Haskell library and tool set for computational biology. The goal of ELynx is-    reproducible research. Evolutionary sequences and phylogenetic trees can be-    read, viewed, modified and simulated. The command line with all arguments is-    logged consistently, and automatically. Data integrity is verified using SHA256-    sums so that validation of past analyses is possible without the need to-    recompute the result.-    -    slynx     Analyze, modify, and simulate evolutionary sequences.-    tlynx     Analyze, modify, and simulate phylogenetic trees.-    elynx     Validate and redo past analyses.-    -    Get help for sub commands:-      slynx examine --help   ## Sub-sample@@ -381,9 +305,9 @@      slynx sub-sample --help -    ELynx Suite version 0.2.1.+    ELynx Suite version 0.3.0.     Developed by Dominik Schrempf.-    Compiled on April 27, 2020, at 12:08 pm, UTC.+    Compiled on July 30, 2020, at 12:40 pm, UTC.          Usage: slynx sub-sample (-a|--alphabet NAME) INPUT-FILE                             (-n|--number-of-sites INT)@@ -403,22 +327,6 @@       -h,--help                Show this help text          Create a given number of multi sequence alignments, each of which contains a given number of random sites drawn from the original multi sequence alignment.-    -    ELynx-    ------    A Haskell library and tool set for computational biology. The goal of ELynx is-    reproducible research. Evolutionary sequences and phylogenetic trees can be-    read, viewed, modified and simulated. The command line with all arguments is-    logged consistently, and automatically. Data integrity is verified using SHA256-    sums so that validation of past analyses is possible without the need to-    recompute the result.-    -    slynx     Analyze, modify, and simulate evolutionary sequences.-    tlynx     Analyze, modify, and simulate phylogenetic trees.-    elynx     Validate and redo past analyses.-    -    Get help for sub commands:-      slynx examine --help   ## Translate@@ -427,9 +335,9 @@      slynx translate --help -    ELynx Suite version 0.2.1.+    ELynx Suite version 0.3.0.     Developed by Dominik Schrempf.-    Compiled on April 27, 2020, at 12:08 pm, UTC.+    Compiled on July 30, 2020, at 12:40 pm, UTC.          Usage: slynx translate (-a|--alphabet NAME) INPUT-FILE (-r|--reading-frame INT)                            (-u|--universal-code CODE)@@ -444,34 +352,17 @@       -u,--universal-code CODE universal code; one of: Standard,                                VertebrateMitochondrial.       -h,--help                Show this help text-    -    -    ELynx-    ------    A Haskell library and tool set for computational biology. The goal of ELynx is-    reproducible research. Evolutionary sequences and phylogenetic trees can be-    read, viewed, modified and simulated. The command line with all arguments is-    logged consistently, and automatically. Data integrity is verified using SHA256-    sums so that validation of past analyses is possible without the need to-    recompute the result.-    -    slynx     Analyze, modify, and simulate evolutionary sequences.-    tlynx     Analyze, modify, and simulate phylogenetic trees.-    elynx     Validate and redo past analyses.-    -    Get help for sub commands:-      slynx examine --help   # TLynx  Handle phylogenetic trees in Newick format. -    tlynx --help+    tlynx --help | head -n -16 -    ELynx Suite version 0.2.1.+    ELynx Suite version 0.3.0.     Developed by Dominik Schrempf.-    Compiled on April 27, 2020, at 12:08 pm, UTC.+    Compiled on July 30, 2020, at 12:40 pm, UTC.          Usage: tlynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME]                   [-f|--force] COMMAND@@ -496,26 +387,10 @@       shuffle                  Shuffle a phylogenetic tree (keep coalescent times, but shuffle topology and leaves).       simulate                 Simulate phylogenetic trees using birth and death processes (see also the 'coalesce' command for simulations using the coalescent process).     +         Available tree file formats:       - Newick Standard: Branch support values are stored in square brackets after branch lengths.       - Newick IqTree:   Branch support values are stored as node names after the closing bracket of forests.-      - Newick RevBayes  Key-value pairs is provided in square brackets after node names as well as branch lengths. XXX: Key value pairs are IGNORED at the moment.-    -    ELynx-    ------    A Haskell library and tool set for computational biology. The goal of ELynx is-    reproducible research. Evolutionary sequences and phylogenetic trees can be-    read, viewed, modified and simulated. The command line with all arguments is-    logged consistently, and automatically. Data integrity is verified using SHA256-    sums so that validation of past analyses is possible without the need to-    recompute the result.-    -    slynx     Analyze, modify, and simulate evolutionary sequences.-    tlynx     Analyze, modify, and simulate phylogenetic trees.-    elynx     Validate and redo past analyses.-    -    Get help for sub commands:-      slynx examine --help   ## Compare@@ -524,12 +399,12 @@      tlynx compare --help -    ELynx Suite version 0.2.1.+    ELynx Suite version 0.3.0.     Developed by Dominik Schrempf.-    Compiled on April 27, 2020, at 12:08 pm, UTC.+    Compiled on July 30, 2020, at 12:40 pm, UTC.          Usage: tlynx compare [-n|--normalize] [-b|--bipartitions] [-t|--intersect] -                         [-f|--newick-format FORMAT] NAME+                         [-f|--newick-format FORMAT] NAMES       Compare two phylogenetic trees (compute distances and branch-wise differences).          Available options:@@ -540,27 +415,11 @@       -t,--intersect           Compare intersections; i.e., before comparison, drop                                leaves that are not present in the other tree       -f,--newick-format FORMAT-                               Newick tree format; see 'tlynx-                               --help' (default: Standard)-      NAME                     Tree file+                               Newick tree format: Standard, IqTree, or RevBayes;+                               default: Standard; for detailed help, see 'tlynx+                               --help'+      NAMES                    Tree files       -h,--help                Show this help text-    -    -    ELynx-    ------    A Haskell library and tool set for computational biology. The goal of ELynx is-    reproducible research. Evolutionary sequences and phylogenetic trees can be-    read, viewed, modified and simulated. The command line with all arguments is-    logged consistently, and automatically. Data integrity is verified using SHA256-    sums so that validation of past analyses is possible without the need to-    recompute the result.-    -    slynx     Analyze, modify, and simulate evolutionary sequences.-    tlynx     Analyze, modify, and simulate phylogenetic trees.-    elynx     Validate and redo past analyses.-    -    Get help for sub commands:-      slynx examine --help   ## Examine@@ -569,9 +428,9 @@      tlynx examine --help -    ELynx Suite version 0.2.1.+    ELynx Suite version 0.3.0.     Developed by Dominik Schrempf.-    Compiled on April 27, 2020, at 12:08 pm, UTC.+    Compiled on July 30, 2020, at 12:40 pm, UTC.          Usage: tlynx examine INPUT-FILE [-f|--newick-format FORMAT]       Compute summary statistics of phylogenetic trees.@@ -581,26 +440,10 @@       -V,--version             Show version       INPUT-FILE               Read trees from INPUT-FILE       -f,--newick-format FORMAT-                               Newick tree format; see 'tlynx-                               --help' (default: Standard)+                               Newick tree format: Standard, IqTree, or RevBayes;+                               default: Standard; for detailed help, see 'tlynx+                               --help'       -h,--help                Show this help text-    -    -    ELynx-    ------    A Haskell library and tool set for computational biology. The goal of ELynx is-    reproducible research. Evolutionary sequences and phylogenetic trees can be-    read, viewed, modified and simulated. The command line with all arguments is-    logged consistently, and automatically. Data integrity is verified using SHA256-    sums so that validation of past analyses is possible without the need to-    recompute the result.-    -    slynx     Analyze, modify, and simulate evolutionary sequences.-    tlynx     Analyze, modify, and simulate phylogenetic trees.-    elynx     Validate and redo past analyses.-    -    Get help for sub commands:-      slynx examine --help   ## Simulate@@ -609,9 +452,9 @@      tlynx simulate --help -    ELynx Suite version 0.2.1.+    ELynx Suite version 0.3.0.     Developed by Dominik Schrempf.-    Compiled on April 27, 2020, at 12:08 pm, UTC.+    Compiled on July 30, 2020, at 12:40 pm, UTC.          Usage: tlynx simulate [-t|--nTrees INT] [-n|--nLeaves INT] [-H|--height DOUBLE]                            [-M|--condition-on-mrca] [-l|--lambda DOUBLE] @@ -640,33 +483,17 @@     Summary statistics only: only print (NumberOfExtantChildren BranchLength) pairs for each branch of each tree. The trees are separated by a newline character.     Sub-sampling: simulate one big tree with n'=round(n/rho), n'>=n, leaves, and randomly sample sub-trees with n leaves. Hence, with rho=1.0, the same tree is reported over and over again.     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-    -    ELynx-    ------    A Haskell library and tool set for computational biology. The goal of ELynx is-    reproducible research. Evolutionary sequences and phylogenetic trees can be-    read, viewed, modified and simulated. The command line with all arguments is-    logged consistently, and automatically. Data integrity is verified using SHA256-    sums so that validation of past analyses is possible without the need to-    recompute the result.-    -    slynx     Analyze, modify, and simulate evolutionary sequences.-    tlynx     Analyze, modify, and simulate phylogenetic trees.-    elynx     Validate and redo past analyses.-    -    Get help for sub commands:-      slynx examine --help   # ELynx  Validate and (optionally) redo past ELynx analyses. -    elynx --help+    elynx --help | head -n -16 -    ELynx Suite version 0.2.1.+    ELynx Suite version 0.3.0.     Developed by Dominik Schrempf.-    Compiled on April 27, 2020, at 12:08 pm, UTC.+    Compiled on July 30, 2020, at 12:40 pm, UTC.          Usage: elynx COMMAND       Validate and redo past ELynx analyses@@ -678,37 +505,4 @@     Available commands:       validate                 Validate an ELynx analysis       redo                     Redo an ELynx analysis-    -    -    ELynx-    ------    A Haskell library and tool set for computational biology. The goal of ELynx is-    reproducible research. Evolutionary sequences and phylogenetic trees can be-    read, viewed, modified and simulated. The command line with all arguments is-    logged consistently, and automatically. Data integrity is verified using SHA256-    sums so that validation of past analyses is possible without the need to-    recompute the result.-    -    slynx     Analyze, modify, and simulate evolutionary sequences.-    tlynx     Analyze, modify, and simulate phylogenetic trees.-    elynx     Validate and redo past analyses.-    -    Get help for sub commands:-      slynx examine --help---# Library documentation--Documentation of the libraries can be found on [Hackage](https://hackage.haskell.org/):---   [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)--Documentation of the executables is also available:---   [elynx](https://hackage.haskell.org/package/elynx)--   [slynx](https://hackage.haskell.org/package/slynx)--   [tlynx](https://hackage.haskell.org/package/tlynx)
bench/Bench.hs view
@@ -1,33 +1,30 @@-{- |-Module      :  Bench-Description :  Various benchmarks-Copyright   :  (c) Dominik Schrempf 2020-License     :  GPL-3.0-or-later--Maintainer  :  dominik.schrempf@gmail.com-Stability   :  unstable-Portability :  portable--Creation date: Mon Dec 16 13:33:27 2019.---}--import           Criterion.Main-import qualified Data.ByteString.Lazy.Char8    as L-import           Data.Tree+-- |+-- Module      :  Bench+-- Description :  Various benchmarks+-- Copyright   :  (c) Dominik Schrempf 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Mon Dec 16 13:33:27 2019.+module Main where -import           ELynx.Data.Tree-import           ELynx.Import.Tree.Newick-import           ELynx.Tools+import Criterion.Main+import Data.ByteString.Lazy.Char8 (ByteString)+import ELynx.Data.Tree+import ELynx.Import.Tree.Newick+import ELynx.Tools  treeFileMany :: FilePath treeFileMany = "data/Many.trees" -getManyTrees :: IO [Tree (PhyloLabel L.ByteString)]-getManyTrees = parseFileWith manyNewick treeFileMany+getManyTrees :: IO (Forest Phylo ByteString)+getManyTrees = parseFileWith (someNewick Standard) treeFileMany  main :: IO () main = do   ts <- getManyTrees   defaultMain-    [bgroup "bipartition" [bench "manyTrees" $ whnf (map bipartitions) ts]]+    [bgroup "bipartition" [bench "manyTrees" $ nf (map bipartitions) ts]]
elynx-tree.cabal view
@@ -1,115 +1,129 @@-cabal-version: 1.12-name: elynx-tree-version: 0.2.2-license: GPL-3-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 <https://github.com/dschrempf/elynx>.-category: Bioinformatics-build-type: Simple+cabal-version: 2.2++-- This file has been generated from package.yaml by hpack version 0.34.2.+--+-- see: https://github.com/sol/hpack+--+-- hash: df9f8f1cb0b1f1f6bac542cb0c44ceb8e3530eb04ea14e29c0b5bb4be0e831f3++name:           elynx-tree+version:        0.3.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+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 extra-source-files:     README.md     ChangeLog.md  source-repository head-    type: git-    location: https://github.com/dschrempf/elynx+  type: git+  location: https://github.com/dschrempf/elynx  library-    exposed-modules:-        ELynx.Data.Tree-        ELynx.Data.Tree.Bipartition-        ELynx.Data.Tree.BranchSupportTree-        ELynx.Data.Tree.Distance-        ELynx.Data.Tree.MeasurableTree-        ELynx.Data.Tree.Multipartition-        ELynx.Data.Tree.NamedTree-        ELynx.Data.Tree.PhyloTree-        ELynx.Data.Tree.SubSample-        ELynx.Data.Tree.SumStat-        ELynx.Data.Tree.Tree-        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.Import.Tree.Newick-        ELynx.Simulate.Coalescent-        ELynx.Simulate.PointProcess-    hs-source-dirs: src-    other-modules:-        Paths_elynx_tree-    default-language: Haskell2010-    ghc-options: -Wall -fllvm-    build-depends:-        base >=4.13.0.0 && <4.14,-        bytestring >=0.10.10.0 && <0.11,-        containers >=0.6.2.1 && <0.7,-        elynx-tools >=0.2.2 && <0.3,-        math-functions >=0.3.3.0 && <0.4,-        megaparsec >=8.0.0 && <8.1,-        mwc-random >=0.14.0.0 && <0.15,-        primitive >=0.7.0.1 && <0.8,-        statistics >=0.15.2.0 && <0.16+  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+  other-modules:+      Paths_elynx_tree+  autogen-modules:+      Paths_elynx_tree+  hs-source-dirs:+      src+  ghc-options: -Wall -fllvm+  build-depends:+      aeson+    , base >=4.7 && <5+    , bytestring+    , comonad+    , containers+    , deepseq+    , elynx-nexus+    , elynx-tools+    , math-functions+    , megaparsec+    , mwc-random+    , primitive+    , statistics+    , vector+  default-language: Haskell2010  test-suite tree-test-    type: exitcode-stdio-1.0-    main-is: Spec.hs-    hs-source-dirs: test-    other-modules:-        ELynx.Data.Tree.BipartitionSpec-        ELynx.Data.Tree.DistanceSpec-        ELynx.Data.Tree.MultipartitionSpec-        ELynx.Data.Tree.PhyloTreeArbitraryInstance-        ELynx.Data.Tree.TreeSpec-        ELynx.Export.Tree.NewickSpec-        ELynx.Import.Tree.NewickSpec-        Paths_elynx_tree-    default-language: Haskell2010-    ghc-options: -Wall -fllvm -threaded -rtsopts -with-rtsopts=-N-    build-depends:-        QuickCheck >=2.13.2 && <2.14,-        base >=4.13.0.0 && <4.14,-        bytestring >=0.10.10.0 && <0.11,-        containers >=0.6.2.1 && <0.7,-        elynx-tools >=0.2.2 && <0.3,-        elynx-tree -any,-        hspec >=2.7.1 && <2.8,-        hspec-megaparsec >=2.1.0 && <2.2,-        math-functions >=0.3.3.0 && <0.4,-        megaparsec >=8.0.0 && <8.1,-        mwc-random >=0.14.0.0 && <0.15,-        primitive >=0.7.0.1 && <0.8,-        quickcheck-instances >=0.3.22 && <0.4,-        statistics >=0.15.2.0 && <0.16+  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.Export.Tree.NewickSpec+      ELynx.Export.Tree.NexusSpec+      ELynx.Import.Tree.NewickSpec+      ELynx.Import.Tree.NexusSpec+      Paths_elynx_tree+  hs-source-dirs:+      test+  ghc-options: -Wall -fllvm+  build-depends:+      QuickCheck+    , base >=4.7 && <5+    , bytestring+    , containers+    , elynx-tools+    , elynx-tree+    , hspec+    , hspec-megaparsec+    , megaparsec+  default-language: Haskell2010  benchmark tree-bench-    type: exitcode-stdio-1.0-    main-is: Bench.hs-    hs-source-dirs: bench-    other-modules:-        Paths_elynx_tree-    default-language: Haskell2010-    ghc-options: -Wall -fllvm -threaded -rtsopts -with-rtsopts=-N-    build-depends:-        base >=4.13.0.0 && <4.14,-        bytestring >=0.10.10.0 && <0.11,-        containers >=0.6.2.1 && <0.7,-        criterion >=1.5.6.2 && <1.6,-        elynx-tools >=0.2.2 && <0.3,-        elynx-tree -any,-        math-functions >=0.3.3.0 && <0.4,-        megaparsec >=8.0.0 && <8.1,-        mwc-random >=0.14.0.0 && <0.15,-        primitive >=0.7.0.1 && <0.8,-        statistics >=0.15.2.0 && <0.16+  type: exitcode-stdio-1.0+  main-is: Bench.hs+  other-modules:+      Paths_elynx_tree+  hs-source-dirs:+      bench+  ghc-options: -Wall -fllvm+  build-depends:+      base >=4.7 && <5+    , bytestring+    , containers+    , criterion+    , elynx-tools+    , elynx-tree+  default-language: Haskell2010
+ src/ELynx/Data/Topology/Phylogeny.hs view
@@ -0,0 +1,115 @@+-- |+-- 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
+ src/ELynx/Data/Topology/Rooted.hs view
@@ -0,0 +1,207 @@+{-# 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
src/ELynx/Data/Tree.hs view
@@ -1,41 +1,50 @@-{- |-Module      :  ELynx.Data.Tree-Description :  Phylogenetic trees-Copyright   :  (c) Dominik Schrempf 2020-License     :  GPL-3.0-or-later+-- TODO: Topology data type.+-- data Topology a = Node (NonEmptySet (Topology a)) | Leaf a -Maintainer  :  dominik.schrempf@gmail.com-Stability   :  unstable-Portability :  portable+-- |+-- 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, -Creation date: Sat Mar 21 16:27:20 2020.+    -- * 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, -module ELynx.Data.Tree-  ( module ELynx.Data.Tree.Tree-  , module ELynx.Data.Tree.PhyloTree-  , module ELynx.Data.Tree.Bipartition-  , module ELynx.Data.Tree.Multipartition-  , module ELynx.Data.Tree.Distance-  , module ELynx.Data.Tree.BranchSupportTree-  , module ELynx.Data.Tree.MeasurableTree-  , module ELynx.Data.Tree.NamedTree-  , module ELynx.Data.Tree.SubSample-  , module ELynx.Data.Tree.SumStat+    -- * 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.BranchSupportTree-import           ELynx.Data.Tree.Distance-import           ELynx.Data.Tree.MeasurableTree-import           ELynx.Data.Tree.Multipartition-import           ELynx.Data.Tree.NamedTree-import           ELynx.Data.Tree.PhyloTree-import           ELynx.Data.Tree.SubSample-import           ELynx.Data.Tree.SumStat-import           ELynx.Data.Tree.Tree+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
src/ELynx/Data/Tree/Bipartition.hs view
@@ -1,198 +1,198 @@-{- |-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@.+{-# LANGUAGE GeneralizedNewtypeDeriving #-} -That's why we have to make sure that for-> Bipartition x y-we always have @x >= y@.+-- |+-- 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, -module ELynx.Data.Tree.Bipartition-  ( -- * The 'Bipartition' data type.-    Bipartition()-  , bps-  , bp-  , bpmap-  , bphuman-  , sshow-    -- * Working with 'Bipartition's.-  , bipartition-  , bipartitions-  , bipartitionToBranchLength-  , bpcompatible+    -- * Work with 'Bipartition's+    bipartition,+    bipartitions,+    getComplementaryLeaves,+    bipartitionToBranch,   ) where -import           Data.List                      ( intercalate )-import qualified Data.Map                      as M-import           Data.Maybe-import qualified Data.Set                      as S-import           Data.Tree+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 -import           ELynx.Data.Tree.Tree+-- | 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. Also the order of the two partitions of the 'Bipartition'--- is not important (see the 'Eq' instance).-newtype Bipartition a = Bipartition {bps :: (S.Set a, S.Set a) -- ^ Tuple of partitions-                                    }-  deriving (Show, Read)+-- 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) --- I decided to NOT provide a human readable show instance because in this case,--- I need the following identity to hold:+-- | Create a bipartition from two sets. ----- > read . show = id+-- Ensure that the smaller set comes first. ----- 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 form. Use a provided function to--- extract the valuable information.-bphuman :: (a -> String) -> Bipartition a -> String-bphuman f (Bipartition (x, y)) = "(" ++ sshow f x ++ "|" ++ sshow f y ++ ")"---- | Show the elements of a subset in a human readable way.-sshow :: (a -> String) -> S.Set a -> String-sshow f = intercalate "," . map f . S.toList+-- 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 subsets.-bp :: Ord a => S.Set a -> S.Set a -> Bipartition a-bp xs ys = if xs >= ys then Bipartition (xs, ys) else Bipartition (ys, xs)+-- | 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) --- | Map a function over all elements in the 'Bipartition's.-bpmap :: (Ord a, Ord b) => (a -> b) -> Bipartition a -> Bipartition b-bpmap f (Bipartition (x, y)) = bp (S.map f x) (S.map f y)+-- | Conversion to a set containing both partitions.+toSet :: Ord a => Bipartition a -> Set a+toSet (Bipartition (x, y)) = S.union x y --- | Create a bipartition from two 'S.Set's.-bpwith :: (Ord a, Ord b) => (a -> b) -> S.Set a -> S.Set a -> Bipartition b-bpwith f x y = bpmap f $ bp 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. -instance (Eq a) => Eq (Bipartition a) where-  Bipartition x == Bipartition y = x == y+-- | 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 ++ ")" -instance (Ord a) => Ord (Bipartition a) where-  Bipartition x `compare` Bipartition y = x `compare` y+-- Show the elements of a set in a human readable format.+setShow :: Show a => Set a -> String+setShow = intercalate "," . map show . S.toList --- Check if a bipartition is valid. For now, only checks if one set is empty.-valid :: Bipartition a -> Bool-valid (Bipartition (xs, ys)) = not $ S.null xs || S.null ys+-- -- | 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.-bipartition :: Ord a => Tree a -> Bipartition a-bipartition (Node _ [x, y]) =-  bp (S.fromList $ leaves x) (S.fromList $ leaves y)-bipartition _ = error "Root node is not bifurcating."+--+-- 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.-bipartitions :: Ord a => Tree a -> S.Set (Bipartition a)-bipartitions t = if S.size (S.fromList ls) == length ls-  then S.filter valid $ bipartitionsUnsafe S.empty pTree-  else error "bipartitions: The tree contains duplicate leaves."- where-  ls    = leaves t-  pTree = partitionTree t+--+-- 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 --- See 'bipartitions', but do not check if leaves are unique.-bipartitionsUnsafe-  :: Ord a => S.Set a -> Tree (S.Set a) -> S.Set (Bipartition a)-bipartitionsUnsafe p (Node l [] ) = S.singleton $ bp p l--- Degree two nodes do not induce additional bipartitions.-bipartitionsUnsafe p (Node _ [x]) = bipartitionsUnsafe p x--- Go through the list of children and combine each of them with the rest.-bipartitionsUnsafe p t@(Node ls xs) =-  S.unions-    $ S.singleton (bp p ls)-    : [ bipartitionsUnsafe lvs x | (lvs, x) <- zip lsOthers xs ]-  where lsOthers = subForestGetSubsets p 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 --- | For a given rose 'Tree', remove all degree two nodes and reconnect the--- resulting disconnected pairs of branches and sum their branch lengths. Since--- the induced bipartitions of the daughter branches of a bifurcating root node--- are equal, the branches are also combined in this case. See--- http://evolution.genetics.washington.edu/phylip/doc/treedist.html and how--- unrooted trees should be handled.------ For this operation, a combining binary function and a unity element is--- required. These requirements are encoded in the 'Monoid' type class--- constraint. Now, each branch on the tree defines a unique 'Bipartition' of--- leaves. Convert a tree into a 'M.Map' from each 'Bipartition' to the length--- of the branch inducing the respective 'Bipartition'. The relevant information--- about the leaves is extracted from the (leaf) nodes with a given function.--- Also check if leaves are unique.-bipartitionToBranchLength-  :: (Ord a, Ord b, Monoid c)-  => (a -> b)      -- ^ Convert node labels to leaves (usually-                                     -- leaf names)-  -> (a -> c)      -- ^ Get length of branch attached to node-  -> Tree a        -- ^ Tree to dissect-  -> M.Map (Bipartition b) c-bipartitionToBranchLength f g t = if S.size (S.fromList ls) == length ls-  then M.filterWithKey (const . valid)-    $ bipartitionToBranchLengthUnsafe (mempty, S.empty) f lAndPTree-  else error "bipartitionToBranchLength: The tree contains duplicate leaves."- where-  ls        = leaves t-  bTree     = fmap g t-  pTree     = partitionTree t-  lAndPTree = fromJust $ merge bTree pTree+-- 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 --- | See 'bipartitionToBranchLength'. When calculating the map, branches--- separated by various degree two nodes have to be combined. Hence, not only--- the complementary partition towards the stem, but also the node label itself--- have to be passed along.-type Info c a = (c, S.Set a)+-- TODO: Unrooted? See module comment of Distance.hs. --- | See 'bipartitionToBranchLength', but does not check if leaves are unique.--- We need information about the nodes, and also about the leaves of the induced--- sub trees. Hence, we need a somewhat complicated node label type+-- | Convert a tree into a 'Map' from each 'Bipartition' to the branch inducing+-- the respective 'Bipartition'. ----- > (a, S.Set a)-bipartitionToBranchLengthUnsafe-  :: (Ord a, Ord b, Monoid c)-  => Info c a-  -> (a -> b)        -- ^ Value to compare on-  -> Tree (Info c a) -- ^ Tree to dissect-  -> M.Map (Bipartition b) c-bipartitionToBranchLengthUnsafe (l, p) f (Node (l', p') []) =-  M.singleton (bpwith f p p') (l <> l')--- The branch length has to be added for degree two nodes.-bipartitionToBranchLengthUnsafe (l, p) f (Node (l', _) [x]) =-  bipartitionToBranchLengthUnsafe (l <> l', p) f x--- Go through the list of children and combine each of them with the rest.-bipartitionToBranchLengthUnsafe (l, p) f t@(Node (l', p') xs) =-  M.unionsWith (<>)-    $ M.singleton (bpwith f p p') (l <> l')-    : [ bipartitionToBranchLengthUnsafe (mempty, lvs) f x-      | (lvs, x) <- zip lvsOthers xs-      ]-  where lvsOthers = subForestGetSubsets p (fmap snd t)+-- 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 --- | Determine compatibility between an bipartition and a subset. If both--- subsets of the bipartition share elements with the given subset, the--- bipartition is incompatible with this subset. If all elements of the subset--- are either not in the bipartition or mapping to one of the two subsets of the--- bipartition, the bipartition and the subset are compatible. See also--- 'ELynx.Data.Tree.Multipartition.compatible'.-bpcompatible :: (Show a, Ord a) => Bipartition a -> S.Set a -> Bool--- compatible (Bipartition (l, r)) ss = sintersection l ss `sdisjoint` sintersection r ss-bpcompatible (Bipartition (l, r)) ss = S.null lOverlap || S.null rOverlap- where-  lOverlap = S.intersection l ss-  rOverlap = S.intersection r ss+-- 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
− src/ELynx/Data/Tree/BranchSupportTree.hs 
@@ -1,77 +0,0 @@-{- |-Module      :  ELynx.Data.Tree.BranchSupportTree-Description :  Node label with branch support-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.Data.Tree.BranchSupportTree-  ( BranchSupport-  , BranchSupported(..)-  , normalizeBranchSupport-  , collapse-  )-where--import           Data.List-import           Data.Tree---- XXX: Implementation of a proper data type. When parsing Newick files, we--- need to check if an Integer can be parsed, otherwise parse a Double. Also, if--- one branch support value is a Double, all others have to be converted to--- Double, so that we do not have a mixed data type. In this case, it may be--- better to use a type family? Or just use Maybe a, and don't restrict the a--- type?------ type BranchSupport = Maybe BSValue------ data BSValue = BSInt Int---              | BSDouble Double---- | At the moment, just use 'Double'. It would be preferable to use a wrapper--- data type that can handle 'Int' or 'Double'.-type BranchSupport = Maybe Double---- | A label that supports extraction and setting of branch support values.-class BranchSupported a where-  -- | For now, branch support is a Double, but one could also think about-  -- bootstrap values, which are integers.-  getBranchSupport :: a -> BranchSupport-  setBranchSupport :: BranchSupport -> a -> a--apply :: BranchSupported a => (Double -> Double) -> a -> a-apply f l = setBranchSupport (f <$> s) l where s = getBranchSupport l---- | Normalize branch support values. The maximum branch support value will be--- set to 1.0.-normalizeBranchSupport :: BranchSupported a => Tree a -> Tree a-normalizeBranchSupport t = case mm of-  Nothing -> t-  Just m  -> fmap (apply (/ m)) t-  where mm = maximum $ fmap getBranchSupport t--accept :: Double -> Maybe Double -> Bool-accept _      Nothing  = True-accept thresh (Just s) = s >= thresh---- | Collapse branches with support lower than given value. Note, branch length--- of collapsed branches is ignored at the moment. Continue collapsing until a--- fix point is reached.-collapse :: (Show a, Eq a, BranchSupported a) => Double -> Tree a -> Tree a-collapse th tr = if tr == tr' then tr else collapse th tr'-  where tr' = collapse' th tr---- | See 'collapse'.-collapse' :: BranchSupported a => Double -> Tree a -> Tree a-collapse' _  t@(Node _ []) = t-collapse' th (  Node l xs) = Node l $ map (collapse' th) (highS ++ lowSubForest)- where-  (highS, lowS) = partition (accept th . getBranchSupport . rootLabel) xs-  lowSubForest  = concatMap subForest lowS
src/ELynx/Data/Tree/Distance.hs view
@@ -1,197 +1,139 @@-{- |-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 binary trees in general).--TODO: All trees are assumed to be UNROOTED. See comments of 'symmetricWith' and-'branchScoreWith', as well as 'bipartitionToBranchLength'.-http://evolution.genetics.washington.edu/phylip/doc/treedist.html. However, this-disagrees with the statement in 'ELynx.Data.Tree.Tree', and should be changed.-I definitely need separate data types for rooted and unrooted trees.--TODO: Use subset module.---}-+-- |+-- 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-  , symmetricWith-  , incompatibleSplits-  , incompatibleSplitsWith-  , branchScore-  , branchScoreWith-  , pairwise-  , adjacent+  ( symmetric,+    incompatibleSplits,+    branchScore,   ) where -import           Data.List-import qualified Data.Map                      as M-import           Data.Monoid-import qualified Data.Set                      as S-import           Data.Tree+-- adjacent, -import           ELynx.Data.Tree.Bipartition    ( Bipartition-                                                , bipartitionToBranchLength-                                                , bipartitions-                                                )-import           ELynx.Data.Tree.MeasurableTree ( Measurable-                                                , getLen-                                                )-import           ELynx.Data.Tree.Multipartition ( Multipartition-                                                , mpcompatible-                                                , fromBipartition-                                                , multipartitions-                                                )-import           ELynx.Data.Tree.NamedTree-import           ELynx.Data.Tree.Tree           ( leaves )+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 => S.Set a -> S.Set a -> S.Set a+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. Before comparing--- the leaf labels, apply a given function. This is useful, for example, to--- compare the labels of 'ELynx.Data.Tree.NamedTree.Named' trees on their names--- only. The tree is assumed to be UNROOTED! See--- http://evolution.genetics.washington.edu/phylip/doc/treedist.html.+-- | Symmetric (Robinson-Foulds) distance between two trees. ----- XXX: Comparing a list of trees with this function recomputes bipartitions.-symmetricWith :: (Show a, Ord b) => (a -> b) -> Tree a -> Tree a -> Int-symmetricWith f t1 t2-  | S.fromList (leaves t1') /= S.fromList (leaves t2') = error-    "symmetricWith: trees do not have equal leaf sets."-  | otherwise = length-  $ symmetricDifference (bipartitions t1') (bipartitions t2')- where-  t1' = fmap f t1-  t2' = fmap f t2---- | See 'symmetricWith', but with 'getName' for comparisons.-symmetric :: (Show a, Ord a, Named a) => Tree a -> Tree a -> Int-symmetric = symmetricWith getName+-- 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-  :: (Ord a, Show a) => S.Set (Bipartition a) -> S.Set (Multipartition a) -> Int-countIncompatibilities bs ms = foldl'-  (\i b -> if any (mpcompatible (fromBipartition b)) ms then i else i + 1)-  0-  bs+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 'symmetricWith' but all--- bipartitions induced by multifurcations are considered. For a detailed--- description of how the distance is calculated, see--- 'ELynx.Data.Tree.Bipartition.bpcompatible'.+-- | 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 tree+-- to get a set of compatible bifurcations of the tree. For example, the star tree ----- > (A,(B,C,D))+-- > (A,B,C,D); -- -- induces the following bipartitions: -- -- > A|BCD -- > B|ACD--- -- > C|ABD -- > D|ABC ----- Those are also reported by the function 'bipartitions'. However, the tree is--- additionally compatible with the following hidden bipartitions:+-- However, the tree is additionally compatible with the following hidden+-- bipartitions: -- -- > AB|CD -- > AC|BD -- > AD|BC ----- For an explanation of how compatibility of a bipartition with a--- multipartition is checked, see 'ELynx.Data.Tree.Multipartition.mpcompatible'.--- Before using 'mpcompatible', bipartitions are simply converted to--- multipartitions with two subsets.+-- For an explanation of how compatibility of partitions is checked, see+-- 'compatible'. Before using 'compatible', bipartitions are simply converted to+-- partitions with two subsets. ----- Only if a bipartition is not compatible with all induced multifurcations of--- the other tree, it is incompatible.+-- 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.-incompatibleSplitsWith :: (Ord b, Show b) => (a -> b) -> Tree a -> Tree a -> Int-incompatibleSplitsWith f t1 t2-  | S.fromList (leaves t1') /= S.fromList (leaves t2')-  = error "incompatibleSplitsWith: trees do not have equal leaf sets."-  | otherwise-  = countIncompatibilities putIncBs1 ms2 + countIncompatibilities putIncBs2 ms1- where-  t1'       = fmap f t1-  t2'       = fmap f t2-  bs1       = bipartitions t1'-  bs2       = bipartitions t2'-  -- Putative incompatible bipartitions of trees one and two, respectively.-  putIncBs1 = bs1 S.\\ bs2-  putIncBs2 = bs2 S.\\ bs1-  -- Multipartitions.-  ms t = multipartitions $ fmap f t-  ms1 = ms t1-  ms2 = ms t2---- | See 'incompatibleSplitsWith', use 'getName' for comparisons.-incompatibleSplits :: (Ord a, Named a) => Tree a -> Tree a -> Int-incompatibleSplits = incompatibleSplitsWith getName+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. Before comparing the leaf--- labels, apply a function. This is useful, for example, to compare the labels--- of 'ELynx.Data.Tree.NamedTree.Named' trees on their names only. The branch--- information which is compared to compute the distance is extracted from the--- nodes with a given function. Assumes that the trees are UNROOTED. See--- http://evolution.genetics.washington.edu/phylip/doc/treedist.html.+-- | 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.-branchScoreWith-  :: (Ord a, Ord b, Floating c)-  => (a -> b) -- ^ Label to compare on-  -> (a -> c) -- ^ Branch information (e.g., length)-                                    -- associated with a node-  -> Tree a-  -> Tree a-  -> c-branchScoreWith f g t1 t2-  | S.fromList (leaves . fmap f $ t1) /= S.fromList (leaves . fmap f $ t2)-  = error "branchScoreWith: trees do not have equal leaf sets."-  | otherwise-  = sqrt dsSquared- where-  bs        = bipartitionToBranchLength f (Sum . g)-  dBs       = M.map getSum $ M.unionWith (-) (bs t1) (bs t2)-  dsSquared = foldl' (\acc e -> acc + e * e) 0 dBs---- | See 'branchScoreWith', use 'getName' and 'getLen' for comparisons.-branchScore :: (Ord a, Named a, Measurable a) => Tree a -> Tree a -> Double-branchScore = branchScoreWith getName getLen---- | Compute pairwise distances of a list of input trees. Use given distance--- measure. Returns a triple, the first two elements are the indices of the--- compared trees, the third is the distance.-pairwise-  :: (a -> a -> b)   -- ^ Distance function-  -> [a]             -- ^ Input trees-  -> [(Int, Int, b)] -- ^ (index i, index j, distance i j)-pairwise dist trs =-  [ (i, j, dist x y)-  | (i : is, x : xs) <- zip (tails [0 ..]) (tails trs)-  , (j     , y     ) <- zip is xs-  ]---- | Compute distances between adjacent pairs of a list of input trees. Use--- given distance measure.-adjacent-  :: (Tree a -> Tree a -> b) -- ^ Distance function-  -> [Tree a]                -- ^ Input trees-  -> [b]-adjacent dist trs = [ dist x y | (x, y) <- zip trs (tail trs) ]+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
+ src/ELynx/Data/Tree/Measurable.hs view
@@ -0,0 +1,148 @@+{-# LANGUAGE OverloadedStrings #-}++-- |+-- 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.+module ELynx.Data.Tree.Measurable+  ( BranchLength,+    Measurable (..),+    applyStem,+    getStem,+    setStem,+    height,+    rootHeight,+    distancesOriginLeaves,+    totalBranchLength,+    summarizeBranchLengths,+    normalizeBranchLengths,+    normalizeHeight,+    ultrametric,+    makeUltrametric,+  )+where++import Data.Bifoldable+import Data.Bifunctor+import Data.ByteString.Lazy.Char8 (ByteString)+import qualified Data.ByteString.Lazy.Char8 as L+import ELynx.Data.Tree.Rooted+import ELynx.Tools+import Text.Printf++-- | Branch length.+type BranchLength = Double++pretty :: BranchLength -> String+pretty = printf "%.5f"++-- | 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)++-- -- | The maximum distance between root and leaves.+-- --+-- -- Does not include the length of the stem, see 'height'.+-- rootHeight :: Measurable e => Tree e a -> BranchLength+-- rootHeight t = subtract (getStem t) $ height 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++prettyRow :: String -> String -> ByteString+prettyRow name val = alignLeft 33 n <> alignRight 8 v+  where+    n = L.pack name+    v = L.pack val++-- | Examine branches of a tree.+summarizeBranchLengths :: Measurable e => Tree e a -> ByteString+summarizeBranchLengths t =+  L.intercalate+    "\n"+    [ prettyRow "Origin height: " $ pretty h,+      prettyRow "Average distance origin to leaves: " $ pretty h',+      prettyRow "Total branch length: " $ pretty b+    ]+  where+    n = length $ leaves t+    h = height t+    h' = sum (distancesOriginLeaves t) / fromIntegral n+    b = totalBranchLength t++-- | 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++-- | 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
− src/ELynx/Data/Tree/MeasurableTree.hs 
@@ -1,158 +0,0 @@-{-# LANGUAGE OverloadedStrings #-}--{- |-Module      :  ELynx.Data.Tree.MeasurableTree-Description :  Functions on trees with branch lengths-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.---}---module ELynx.Data.Tree.MeasurableTree-  ( Measurable(..)-  , extendBranchLength-  , distancesRootLeaves-  , distancesOriginLeaves-  , averageDistanceOriginLeaves-  , height-  , rootHeight-  , lengthenStem-  , shortenStem-  , summarize-  , totalBranchLength-  , normalizeBranchLength-  , prune-  , removeMultifurcations-  , ultrametric-  )-where--import qualified Data.ByteString.Lazy.Char8    as L-import           Data.Foldable-import           Data.Tree-import           Text.Printf--import           ELynx.Data.Tree.Tree-import           ELynx.Tools---- | A 'Node' label with measurable and modifiable branch length to the parent.-class Measurable a where-  -- | Length of attached branch.-  getLen :: a -> Double-  -- | Set attached branch length.-  setLen :: Double -> a -> a--lengthen :: Measurable a => Double -> a -> a-lengthen dl l = setLen (dl + getLen l) l---- | @extendBranchLength daughter parent@ takes the daughter node and extends--- the branch by the length obtained from the parent node.-extendBranchLength :: Measurable a => a -> a -> a-extendBranchLength da pa = lengthen (getLen pa) da---- | Elongate branch length.--- -- | Shorten branch length.--- shorten :: Double -> a -> a--- shorten dl = lengthen (-dl)---- | Distances from the root node of a tree to its leaves (this are not the--- distances from the origin to the leaves, see 'distancesOriginLeaves').-distancesRootLeaves :: (Measurable a) => Tree a -> [Double]-distancesRootLeaves (Node _ []) = [0]-distancesRootLeaves (Node _ f) =-  concat [ map (+ getLen (rootLabel d)) (distancesRootLeaves d) | d <- f ]---- -- | Distances from the origin of a tree to its leaves (this is not the distance--- -- from the root node to the leaves, which would be @distanceOriginLeaves t ---- -- (getLen $ rootLabel t)@.).--- distancesOriginLeaves :: (Measurable a) => Tree a -> [Double]--- distancesOriginLeaves (Node l []) = [getLen l]--- distancesOriginLeaves (Node l f ) = concatMap (map (+ getLen l) . distancesOriginLeaves) f---- | Distances from the origin of a tree to its leaves (this is not the distance--- from the root node to the leaves, which would be @distanceOriginLeaves t ---- (getLen $ rootLabel t)@.). See also 'distancesRootLeaves'.-distancesOriginLeaves :: (Measurable a) => Tree a -> [Double]-distancesOriginLeaves t@(Node l _) = map (+ getLen l) (distancesRootLeaves t)---- | Average distance from the origin of a tree to its leaves, see--- 'distancesOriginLeaves'.-averageDistanceOriginLeaves :: (Measurable a) => Tree a -> Double-averageDistanceOriginLeaves tr = sum ds / fromIntegral n- where-  ds = distancesOriginLeaves tr-  n  = length ds---- | Height (max distance between origin and leaves) of a tree. Return 0 if the--- tree is empty.-height :: (Measurable a) => Tree a -> Double-height = maximum . distancesOriginLeaves---- | Height of root node. Return 0 if the tree is empty.-rootHeight :: (Measurable a) => Tree a -> Double-rootHeight = maximum . distancesRootLeaves---- | Lengthen the distance between root and origin.-lengthenStem :: (Measurable a) => Double -> Tree a -> Tree a-lengthenStem dl (Node lbl chs) = Node (lengthen dl lbl) chs---- | Lengthen the distance between root and origin.-shortenStem :: (Measurable a) => Double -> Tree a -> Tree a-shortenStem dl = lengthenStem (-dl)--pRow :: String -> String -> L.ByteString-pRow name val = alignLeft 33 n <> alignRight 8 v- where-  n = L.pack name-  v = L.pack val---- | Summarize a tree with measureable branch lengths.-summarize :: (Measurable a) => Tree a -> L.ByteString-summarize t = L.intercalate-  "\n"-  [ pRow "Leaves: " $ show n-  , pRow "Height: " $ printf "%.5f" h-  , pRow "Average distance root to leaves: " $ printf "%.5f" h'-  , pRow "Total branch length: " $ printf "%.5f" b-  ]- where-  n  = length . leaves $ t-  h  = height t-  b  = totalBranchLength t-  h' = sum (distancesOriginLeaves t) / fromIntegral n---- | Total branch length of a tree.-totalBranchLength :: (Measurable a) => Tree a -> Double-totalBranchLength = foldl' (\acc n -> acc + getLen n) 0---- | Normalize tree so that sum of branch lengths is 1.0.-normalizeBranchLength :: (Measurable a) => Tree a -> Tree a-normalizeBranchLength t = fmap (\n -> setLen (getLen n / s) n) t-  where s = totalBranchLength t---- | Prune degree 2 nodes. Add branch lengths but forget pruned node label. See--- 'pruneWith'.-prune :: (Measurable a) => Tree a -> Tree a-prune = pruneWith f where f da pa = lengthen (getLen pa) da---- | Remove multifurcations by copying multifurcating nodes and introducing--- branches of length 0.-removeMultifurcations :: Measurable a => Tree a -> Tree a-removeMultifurcations t@(Node _ [] ) = t-removeMultifurcations (  Node l [x]) = Node l [removeMultifurcations x]-removeMultifurcations (Node l [x, y]) =-  Node l $ map removeMultifurcations [x, y]-removeMultifurcations (Node l (x : xs)) = Node l-  $ map removeMultifurcations [x, Node l' xs]-  where l' = setLen 1.0 l---- | Check if a tree is ultrametric.-ultrametric :: Measurable a => Tree a -> Bool-ultrametric = allNearlyEqual . distancesOriginLeaves
− src/ELynx/Data/Tree/Multipartition.hs 
@@ -1,166 +0,0 @@-{- |-Module      :  ELynx.Data.Tree.Multipartition-Description :  Multipartitions 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 'Multipartition', similar to branches inducing-'ELynx.Data.Tree.Bipartition's.---}--module ELynx.Data.Tree.Multipartition-  (-          -- * The 'Multipartition' data type.-    Multipartition()-  , mps-  , mp-  , mpmap-  , mphuman-  , fromBipartition-        -- * Working with 'Multipartition's.-  , multipartitions-  , findSubset-  , mpcompatible-  )-where--import           Data.List                      ( find-                                                , foldl'-                                                , intercalate-                                                )-import           Data.Maybe                     ( fromMaybe )-import qualified Data.Set                      as S-import           Data.Tree--import           ELynx.Data.Tree.Bipartition    ( Bipartition-                                                , bps-                                                , sshow-                                                )-import           ELynx.Data.Tree.Tree---- | Each branch of a bifurcating tree partitions the leaves of the tree into--- three subsets, see 'ELynx.Data.Tree.Bipartition'. In a similar way, each--- internal node induces a tripartition. Tripartitions are not yet implemented--- (December 2019) because it is usually sufficient to work with bipartitions.--- If, however, 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,--- and a 'Multipartition'. Multipartitions are interesting in that we can use--- them for calculating incompatible splits, see 'ELynx.Data.Tree.Distance'. The--- order of the partitions within a multipartition is unimportant, .-newtype Multipartition a = Multipartition { mps :: S.Set (S.Set a) -- ^ Set of partitions-                                          }-  deriving (Show, Read)---- | Show a multipartition in a human readable form. Use a provided function to--- extract the valuable information.-mphuman :: (a -> String) -> Multipartition a -> String-mphuman f (Multipartition xs) =-  "(" ++ intercalate "|" (map (sshow f) (S.toList xs)) ++ ")"---- | Create a multipartition.-mp :: Ord a => [S.Set a] -> Multipartition a-mp = mp' . filter (not . S.null)--mp' :: Ord a => [S.Set a] -> Multipartition a-mp' xs = Multipartition (S.fromList xs)---- | Map a function over all elements in the multipartitions.-mpmap :: (Ord a, Ord b) => (a -> b) -> Multipartition a -> Multipartition b-mpmap f (Multipartition xs) = Multipartition $ S.map (S.map f) xs--instance (Eq a) => Eq (Multipartition a) where-  Multipartition xs == Multipartition ys = xs == ys--instance (Ord a) => Ord (Multipartition a) where-  Multipartition xs `compare` Multipartition ys = xs `compare` ys---- | Convert bipartition to multipartition.-fromBipartition :: Ord a => Bipartition a -> Multipartition a-fromBipartition bp = mp [l, r] where (l, r) = bps bp---- | Get all multipartitions of a tree.-multipartitions :: Ord a => Tree a -> S.Set (Multipartition a)-multipartitions t = if S.size (S.fromList lvs) == length lvs-  then multipartitionsUnsafe S.empty (partitionTree t)-  else error "multipartitions: The tree contains duplicate leaves."-  where lvs = leaves t---- | See 'multipartitions', but do not check if leaves are unique.-multipartitionsUnsafe-  :: Ord a => S.Set a -> Tree (S.Set a) -> S.Set (Multipartition a)-multipartitionsUnsafe _  (Node _ []    ) = S.empty-multipartitionsUnsafe xs (Node _ [x]   ) = multipartitionsUnsafe xs x-multipartitionsUnsafe xs (Node _ [x, y]) = S.union l r- where-  l = multipartitionsUnsafe (S.union xs (rootLabel x)) y-  r = multipartitionsUnsafe (S.union xs (rootLabel y)) x-multipartitionsUnsafe xs t@(Node _ ys) =-  S.unions-    $ S.singleton (mp (xs : map rootLabel ys))-    : zipWith multipartitionsUnsafe lvsOthers ys-  where lvsOthers = subForestGetSubsets xs t---- | Find the subset of a multipartition containing a given element.-findSubset :: Ord a => a -> Multipartition a -> S.Set a-findSubset l m = -- Return the empty subset if nothing is found. This corresponds-                 -- to having no information about the leaf in question.-                 fromMaybe S.empty (find (S.member l) ss)-  where ss = mps m---- Add the subset of a bipartition which contains a given element.-addSubset-  :: Ord a => Multipartition a -> S.Set (S.Set a) -> a -> S.Set (S.Set a)-addSubset m ss l = if not $ S.null s then s `S.insert` ss else ss-  where s = findSubset l m---- Each subset overlaps with a number of subsets of a bipartition which are--- returned by this function.-overlap :: Ord a => Multipartition a -> S.Set a -> S.Set (S.Set a)-overlap m = foldl' (addSubset m) S.empty---- | Multipartitions are compatible if they do not contain conflicting--- information. This function checks if two multipartitions are compatible with--- each other. Thereby, following algorithm is used:------ 1. Take each subset of the first multipartition.------ 2a. Determine the overlap: For each leaf of the chosen subset, add the subset--- of the second multipartition containing the leaf. The result is a set of--- subsets, which is the union of the added subsets.------ The data type "set of subsets" is actually the same data type as a--- multipartition. However, it is not a partition, because it may and will not--- span the whole set of leaves, and so, I use @S.Set (S.Set a)@. One could--- define a multiset data type to improve comprehensibility.------ 2b. Collect the set of subsets from point 1.------ 3. Each set of subsets needs to be either equal or disjoint with any other--- set of subsets in the collection. If so, the first multipartition is--- compatible with the second.------ 4. Exchange the first with the second multipartition and go through steps 1--- to 3.------ See also 'ELynx.Data.Tree.Bipartition.compatible'.-mpcompatible :: (Ord a, Show a) => Multipartition a -> Multipartition a -> Bool-mpcompatible l r =-  and-    $  [ x `S.disjoint` y | x <- lOverlaps, y <- lOverlaps, x /= y ]-    ++ [ x `S.disjoint` y | x <- rOverlaps, y <- rOverlaps, x /= y ]- where-  ls        = S.toList $ mps l-  rs        = S.toList $ mps r-  -- The subsets on the left multipartition overlap the subsets of the-  -- right multipartition.-  lOverlaps = map (overlap r) ls-  -- The subsets on the left multipartition overlap the subsets of the-  -- right multipartition.-  rOverlaps = map (overlap l) rs
+ src/ELynx/Data/Tree/Named.hs view
@@ -0,0 +1,34 @@+-- |+-- 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 Data.ByteString.Lazy.Builder (char8, doubleDec, intDec, toLazyByteString)+import Data.ByteString.Lazy.Char8 (ByteString)++-- | Data types with names.+class Named a where+  getName :: a -> ByteString++instance Named Int where+  getName = toLazyByteString . intDec++instance Named Double where+  getName = toLazyByteString . doubleDec++instance Named Char where+  getName = toLazyByteString . char8++instance Named ByteString where+  getName = id
− src/ELynx/Data/Tree/NamedTree.hs 
@@ -1,38 +0,0 @@-{- |-Module      :  ELynx.Data.Tree.NamedTree-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.NamedTree-  ( Named(..)-  )-where--import qualified Data.ByteString.Lazy.Builder  as L-import qualified Data.ByteString.Lazy.Char8    as L---- | Data types with names.-class Named a where-  getName :: a -> L.ByteString--instance Named Int where-  getName = L.toLazyByteString . L.intDec--instance Named Double where-  getName = L.toLazyByteString . L.doubleDec--instance Named Char where-  getName = L.toLazyByteString . L.char8--instance Named L.ByteString where-  getName = id-
+ src/ELynx/Data/Tree/Partition.hs view
@@ -0,0 +1,137 @@+-- |+-- 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`)
− src/ELynx/Data/Tree/PhyloTree.hs 
@@ -1,67 +0,0 @@-{-# LANGUAGE FlexibleInstances #-}--{- |-Module      :  ELynx.Data.Tree.PhyloTree-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.--Phylogenetic nodes labels, aka 'PhyloLabel's, have a branch length and an-arbitrary label type, e.g., of type 'Int'.---}---module ELynx.Data.Tree.PhyloTree-  ( PhyloLabel(..)-  , removeBrInfo-  )-where--import           Data.Function-import           Data.Maybe                     ( fromMaybe )-import           Data.Tree--- import           Test.QuickCheck         hiding ( label )--import           ELynx.Data.Tree.BranchSupportTree-import           ELynx.Data.Tree.MeasurableTree-import           ELynx.Data.Tree.NamedTree---- | A primitive label type for phylogenetic trees with a name, possibly a--- branch support value, and possibly a branch length.-data PhyloLabel a = PhyloLabel { label :: a-                               , brSup :: Maybe Double-                               , brLen :: Maybe Double }-                  deriving (Read, Show, Eq)--instance Ord a => Ord (PhyloLabel a) where-  compare = compare `on` label---- | If no branch length is available, 0 is returned. This is probably not the--- best (and final) behavior.-instance Measurable (PhyloLabel a) where-  getLen = fromMaybe 0 . brLen-  setLen l x-    | l >= 0    = x { brLen = Just l }-    | otherwise = error $ "Branch lengths cannot be negative: " <> show l--instance BranchSupported (PhyloLabel a) where-  getBranchSupport = brSup-  setBranchSupport Nothing l = l { brSup = Nothing }-  setBranchSupport (Just s) l-    | s > 0     = l { brSup = Just s }-    | otherwise = error "Branch support cannot be negative."--instance Named a => Named (PhyloLabel a) where-  getName = getName . label---- | Remove all branch relevant information from all nodes of the tree; only--- retain the labels.-removeBrInfo :: Tree (PhyloLabel a) -> Tree a-removeBrInfo = fmap label
+ src/ELynx/Data/Tree/Phylogeny.hs view
@@ -0,0 +1,476 @@+{-# 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 enforced 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 'roots', or 'rootAt'.+--+-- Trees with multifurcating root nodes can be properly rooted using 'outgroup'.+module ELynx.Data.Tree.Phylogeny+  ( -- * Functions+    equal,+    intersect,+    bifurcating,+    outgroup,+    midpoint,+    roots,+    rootAt,++    -- * Branch labels+    Phylo (..),+    Length (..),+    phyloToLengthTree,+    lengthToPhyloTree,+    Support (..),+    phyloToSupportTree,+    phyloToSupportTreeUnsafe,+    PhyloStrict (..),+    toStrictTree,+    fromStrictTree,+  )+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++-- TODO. Adapt for trees with branches.+-- -- | 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 e a -> Tree e 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.+--+-- 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)]++-- XXX: 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++-- | 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++-- | Set all branch support values to 'Nothing'.+--+-- Useful, for example, to export a tree with branch lengths in Newick format.+lengthToPhyloTree :: Tree Length a -> Tree Phylo a+lengthToPhyloTree = first fromLengthLabel++fromLengthLabel :: Length -> Phylo+fromLengthLabel (Length b) = Phylo (Just b) Nothing++-- | 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++-- TODO: Change name. Strict is reserved for "not lazy".++-- | Strict branch label for phylogenetic trees.+data PhyloStrict = PhyloStrict+  { sBrLen :: BranchLength,+    sBrSup :: BranchSupport+  }+  deriving (Read, Show, Eq, Ord, Generic)++instance Semigroup PhyloStrict where+  PhyloStrict bL sL <> PhyloStrict bR sR = PhyloStrict (bL + bR) (min sL sR)++instance Measurable PhyloStrict where+  getLen = sBrLen+  setLen b l = l {sBrLen = b}++instance Splittable PhyloStrict where+  split l = l {sBrLen = b'}+    where+      b' = sBrLen l / 2.0++instance Supported PhyloStrict where+  getSup = sBrSup+  setSup s l = l {sBrSup = s}++instance ToJSON PhyloStrict++instance FromJSON PhyloStrict++-- | Conversion to a 'PhyloStrict' tree.+--+-- See 'phyloToLengthTree' and 'phyloToSupportTree'.+toStrictTree :: Tree Phylo a -> Either String (Tree PhyloStrict a)+toStrictTree t = do+  lt <- first fromLength <$> phyloToLengthTree t+  st <- first fromSupport <$> phyloToSupportTree t+  case zipTreesWith PhyloStrict const lt st of+    Nothing -> error "toStrictTree: This is a bug. Can not zip two trees with the same topology."+    Just zt -> return zt++-- | Set all branch length and support values to 'Just' the value.+--+-- Useful, for example, to export a tree with branch lengths in Newick format.+fromStrictTree :: Tree PhyloStrict a -> Tree Phylo a+fromStrictTree = first fromStrictLabel++fromStrictLabel :: PhyloStrict -> Phylo+fromStrictLabel (PhyloStrict b s) = Phylo (Just b) (Just s)
+ src/ELynx/Data/Tree/Rooted.hs view
@@ -0,0 +1,355 @@+{-# 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
+ src/ELynx/Data/Tree/Splittable.hs view
@@ -0,0 +1,29 @@+-- |+-- 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)
− src/ELynx/Data/Tree/SubSample.hs 
@@ -1,64 +0,0 @@-{- |-Module      :  ELynx.Data.Tree.SubSample-Description :  Sub sample trees-Copyright   :  (c) Dominik Schrempf 2020-License     :  GPL-3.0-or-later--Maintainer  :  dominik.schrempf@gmail.com-Stability   :  unstable-Portability :  portable--Creation date: Tue Feb 11 15:56:01 2020.--TODO: Move this to the sub-sample command file? Somehow this does not belong in-the library.---}--module ELynx.Data.Tree.SubSample-  ( subSample-  , nSubSamples-  )-where--import           Control.Monad                  ( replicateM )-import           Control.Monad.Primitive-import qualified Data.Sequence                 as Seq-import qualified Data.Set                      as Set-import           Data.Tree-import           System.Random.MWC--import           ELynx.Data.Tree.Tree           ( subTree )-import           ELynx.Tools---- | Extract a random subtree with @N@ leaves of a tree with @M@ leaves, where--- @M>N@ (otherwise error). The complete list of leaves (names are assumed to be--- unique) has to be provided as a 'Seq.Seq', and a 'Seq.Set', so that fast--- sub-sampling as well as lookup are fast and so that these data structures do--- not have to be recomputed when many sub-samples are requested.-subSample-  :: (PrimMonad m, Ord a)-  => Seq.Seq a-  -> Int-  -> Tree a-  -> Gen (PrimState m)-  -> m (Maybe (Tree a))-subSample lvs n tree g-  | Seq.length lvs < n = error-    "Given list of leaves is shorter than requested number of leaves."-  | otherwise = do-    sampledLs <- sample lvs n g-    let ls = Set.fromList sampledLs-    return $ subTree (`Set.member` ls) tree---- | See 'subSample', but @n@ times.-nSubSamples-  :: (PrimMonad m, Ord a)-  => Int-  -> Seq.Seq a-  -> Int-  -> Tree a-  -> Gen (PrimState m)-  -> m [Maybe (Tree a)]-nSubSamples nS lvs nL tree g = replicateM nS $ subSample lvs nL tree g-
− src/ELynx/Data/Tree/SumStat.hs 
@@ -1,60 +0,0 @@-{- |-Module      :  ELynx.Data.Tree.SumStat-Description :  Summary statistics for phylogenetic trees-Copyright   :  (c) Dominik Schrempf 2018-License     :  GPL-3.0-or-later--Maintainer  :  dominik.schrempf@gmail.com-Stability   :  unstable-Portability :  portable--Creation date: Thu May 17 14:05:45 2018.--TODO: Move this to the sub-sample command file? Somehow this does not belong in-the library.---}--module ELynx.Data.Tree.SumStat-  ( BrLnNChildren-  , NChildSumStat-  , toNChildSumStat-  , formatNChildSumStat-  )-where--import qualified Data.ByteString.Builder       as L-import qualified Data.ByteString.Lazy.Char8    as L-import           Data.Tree-import           ELynx.Data.Tree.MeasurableTree---- This may be too specific, but I only change it if necessary. E.g., use types--- a (for node labels) and b (for branch lengths).---- | Pair of branch length with number of extant children.-type BrLnNChildren = (Double, Int)---- | Possible summary statistic of phylogenetic trees. A list of tuples--- (BranchLength, NumberOfExtantChildrenBelowThisBranch).-type NChildSumStat = [BrLnNChildren]---- | Format the summary statistics in the following form:--- @---    nLeaves1 branchLength1---    nLeaves2 branchLength2---    ....-formatNChildSumStat :: NChildSumStat -> L.ByteString-formatNChildSumStat s =-  L.toLazyByteString . mconcat $ map formatNChildSumStatLine s--formatNChildSumStatLine :: BrLnNChildren -> L.Builder-formatNChildSumStatLine (l, n) =-  L.intDec n <> L.char8 ' ' <> L.doubleDec l <> L.char8 '\n'---- | Compute NChilSumStat for a phylogenetic tree.-toNChildSumStat :: Measurable a => Tree a -> NChildSumStat-toNChildSumStat (Node lbl []) = [(getLen lbl, 1)]-toNChildSumStat (Node lbl ts) = (getLen lbl, sumCh) : concat nChSS- where-  nChSS = map toNChildSumStat ts-  sumCh = sum $ map (snd . head) nChSS
+ src/ELynx/Data/Tree/Supported.hs view
@@ -0,0 +1,61 @@+-- |+-- 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.+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++-- TODO: Something was wrong here. @collapse 1.0 t@ should be a star tree but it+-- was a leaf. Is this still so?++-- | 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'++-- See 'collapse'.+collapse' :: Supported e => BranchSupport -> Tree e a -> Tree e a+collapse' _ t@(Node _ _ []) = t+collapse' th (Node br lb ts) = Node br lb $ map (collapse' th) (highSupport ++ lowSupportForest)+  where+    (highSupport, lowSupport) = partition ((>= th) . getSup . branch) ts+    lowSupportForest = concatMap forest lowSupport
− src/ELynx/Data/Tree/Tree.hs 
@@ -1,268 +0,0 @@-{- |-Module      :  ELynx.Data.Tree.Tree-Description :  Functions related to 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 09:57:29 2019.--Functions to work with rooted, rose 'Tree's with labeled, unique leaves. The-order of children is not important. However, for the moment, equality checks and-other comparisons are sensitive to the order of children. See TODO note below.--Comment about nomenclature:--- In "Data.Tree", a 'Tree' is defined as--@-data Tree a = Node {-        rootLabel :: a,         -- ^ label value-        subForest :: Forest a   -- ^ zero or more child trees-    }-@--This means, that the word 'Node' is reserved for the constructor of a tree, and-that a 'Node' has a label and a children. The terms 'Node' and /label/ are not-to be confused.--- Branches have /lengths/. For example, a branch length can be a distances or a-  time.--NOTE: Trees in this library are all rooted. Unrooted trees can be treated with a-rooted data structure equally well. However, in these cases, some functions have-no meaning. For example, functions measuring the distance from the root to the-leaves (the height of a rooted tree).--TODO: The 'Tree' data type is a rose tree with an ordered sub-forest. However,-the order of the sub-forest does not matter for phylogenetic trees. Equality-checks will throw false negatives the compared trees only differ in their orders-of sub-trees.--NOTE: Try fgl or alga. Use functional graph library for unrooted trees see also-the book /Haskell high performance programming from Thomasson/, p. 344.---}---module ELynx.Data.Tree.Tree-  ( singleton-  , degree-  , leaves-  , pruneWith-  , dropLeafWith-  , intersectWith-  , merge-  , tZipWith-  , partitionTree-  , subForestGetSubsets-  , subTree-  , bifurcating-  , roots-  , connect-  , clades-  )-where--import           Control.Monad                  ( zipWithM )-import           Data.List                      ( foldl'-                                                , foldl1'-                                                )-import           Data.Maybe                     ( mapMaybe )-import qualified Data.Set                      as S-import           Data.Traversable-import           Data.Tree----- | The simplest tree. Usually an extant leaf.-singleton :: a -> Tree a-singleton l = Node l []---- | The degree of the root node of a tree.-degree :: Tree a -> Int-degree = (+ 1) . length . subForest---- | Get leaves of tree.-leaves :: Tree a -> [a]-leaves (Node l []) = [l]-leaves (Node _ f ) = concatMap leaves f---- | Prune or remove degree two inner nodes. The information stored in a pruned--- node can be used to change the daughter node. To discard this information,--- use, @pruneWith const@, otherwise @pruneWith (\daughter parent -> combined)@.-pruneWith :: (a -> a -> a) -> Tree a -> Tree a-pruneWith _ n@(Node _ []) = n-pruneWith f (Node paLbl [ch]) =-  let lbl = f (rootLabel ch) paLbl in pruneWith f $ Node lbl (subForest ch)-pruneWith f (Node paLbl chs) = Node paLbl (map (pruneWith f) chs)---- | Drop a leaf from a tree with unique leaf names. The possibly resulting--- degree two node is pruned with 'pruneWith'. Two functions are given for node--- name extraction, and for the combination of possibly resulting degree two--- nodes.-dropLeafWith-  :: (Show b, Ord b) => (a -> b) -> (a -> a -> a) -> b -> Tree a -> Tree a-dropLeafWith f g l t-  | l `notElem` lvs = error "dropLeafWith: leaf not found on tree."-  | S.size (S.fromList lvs) < length lvs = error-    "dropLeafWith: tree does not have unique leaves."-  | otherwise = dropLeafWithUnsafe f g l t-  where lvs = leaves $ fmap f t---- See 'dropLeafWith'.-dropLeafWithUnsafe :: Eq b => (a -> b) -> (a -> a -> a) -> b -> Tree a -> Tree a-dropLeafWithUnsafe f g lf (Node x xs)-  | length xs' == 1 = let Node z zs = head xs' in Node (g z x) zs-  | otherwise       = Node x xs'- where-  isThisLeaf y = null (subForest y) && f (rootLabel y) == lf-  xs' = map (dropLeafWithUnsafe f g lf) (filter (not . isThisLeaf) xs)---- | Compute the intersection of trees. The intersections are the largest--- subtrees sharing the same leaf set. Leaf names used for comparison are--- extracted by a given function. Leaves are dropped with 'dropLeafWith', and--- degree two nodes are pruned with 'pruneWith'.-intersectWith-  :: (Show b, Ord b) => (a -> b) -> (a -> a -> a) -> [Tree a] -> [Tree a]-intersectWith f g ts = if null ls-  then error "intersect: intersection of leaves is empty."-  else map (retainLeavesWith f g ls) ts- where -- Leaf sets.-  lss = map (S.fromList . leaves . fmap f) ts-  -- Common leaf set.-  ls  = foldl1' S.intersection lss---- Retain all leaves in a provided set; or conversely, drop all leaves not in a--- provided set.-retainLeavesWith-  :: (Show b, Ord b) => (a -> b) -> (a -> a -> a) -> S.Set b -> Tree a -> Tree a-retainLeavesWith f g ls t = foldl' (flip (dropLeafWith f g)) t leavesToDrop-  where leavesToDrop = filter (`S.notMember` ls) $ leaves $ fmap f t---- | Merge two trees with the same topology. Returns 'Nothing' if the topologies--- are different.-merge :: Tree a -> Tree b -> Maybe (Tree (a, b))-merge (Node l xs) (Node r ys) = if length xs == length ys-                                -- I am proud of that :)).-  then zipWithM merge xs ys >>= Just . Node (l, r)-  else Nothing---- | Apply a function with different effect on each node to a 'Traversable'.--- Based on https://stackoverflow.com/a/41523456.-tZipWith :: Traversable t => (a -> b -> c) -> [a] -> t b -> Maybe (t c)-tZipWith f xs = sequenceA . snd . mapAccumL pair xs- where-  pair []       _ = ([], Nothing)-  pair (y : ys) z = (ys, Just (f y z))---- | Each node of a tree is root of a subtree. Get the leaves of the subtree of--- each node.-partitionTree :: (Ord a) => Tree a -> Tree (S.Set a)-partitionTree (Node l []) = Node (S.singleton l) []-partitionTree (Node _ xs) = Node (S.unions $ map rootLabel xs') xs'-  where xs' = map partitionTree xs---- | Get subtree of 'Tree' with nodes satisfying predicate. Return 'Nothing', if--- no leaf satisfies predicate. At the moment: recursively, for each child, take--- the child if any leaf in the child satisfies the predicate.-subTree :: (a -> Bool) -> Tree a -> Maybe (Tree a)-subTree p leaf@(Node lbl []) | p lbl     = Just leaf-                             | otherwise = Nothing-subTree p (Node lbl chs) = if null subTrees-  then Nothing-  else Just $ Node lbl subTrees-  where subTrees = mapMaybe (subTree p) chs---- | Loop through each tree in a forest to report the complementary leaf sets.-subForestGetSubsets-  :: (Ord a)-  => S.Set a          -- ^ Complementary partition at the stem-  -> Tree (S.Set a)   -- ^ Tree with partition nodes-  -> [S.Set a]-subForestGetSubsets lvs t = lvsOthers- where-  xs          = subForest t-  nChildren   = length xs-  lvsChildren = map rootLabel xs-  lvsOtherChildren =-    [ S.unions $ lvs : take i lvsChildren ++ drop (i + 1) lvsChildren-    | i <- [0 .. (nChildren - 1)]-    ]-  lvsOthers = map (S.union lvs) lvsOtherChildren---- | Check if a tree is bifurcating. A Bifurcating tree only contains degree one--- and degree three nodes. I know, one should use a proper data structure to--- encode bifurcating trees.-bifurcating :: Tree a -> Bool-bifurcating (Node _ []    ) = True-bifurcating (Node _ [_]   ) = False-bifurcating (Node _ [x, y]) = bifurcating x && bifurcating y-bifurcating (Node _ _     ) = False---- TODO: This bifurcating stuff irritates me. There are two solutions:------ 1. Use a bifurcating data structure.------ 2. Do not move down multifurcations.------ Solution 1 is pretty, but doesn't allow for what we actually want to do,--- which is solution 2. We want to encode clades that for sure do not contain--- the root as multifurcations.---- | For a rooted, bifurcating tree, get all possible rooted trees. For a tree--- with @n>2@ leaves, there are @(2n-3)@ rooted trees. Beware, a bifurcating--- tree without degree two nodes is assumed (see 'bifurcating'). The root node--- is moved.-roots :: Tree a -> [Tree a]--- Leaves, and cherries have to be handled separately, because they cannot be--- rotated.-roots t@(Node _ []                    ) = [t]-roots t@(Node _ [Node _ [], Node _ []]) = [t]-roots t | bifurcating t = t : left t ++ right t-        | otherwise     = error "roots: Tree is not bifurcating."---- Move the root to the left.-left :: Tree a -> [Tree a]-left (Node i [Node j [x], z]) =-  let t' = Node i [x, Node j [z]] in t' : left t'-left (Node i [Node j [x, y], z]) =-  let tll = Node i [x, Node j [y, z]]-      tlr = Node i [Node j [x, z], y]-  in  tll : tlr : left tll ++ right tlr-left (Node _ [Node _ [], _]) = []-left (Node _ []            ) = error "left: Encountered a leaf."-left _                       = error "left: Tree is not bifurcating."---- Move the root to the right.-right :: Tree a -> [Tree a]-right (Node i [x, Node j [z]]) =-  let t' = Node i [Node j [x], z] in t' : right t'-right (Node i [x, Node j [y, z]]) =-  let trl = Node i [y, Node j [x, z]]-      trr = Node i [Node j [x, y], z]-  in  trl : trr : left trl ++ right trr-right (Node _ [_, Node _ []]) = []-right (Node _ []            ) = error "right: Encountered a leaf."-right (Node _ [_]) =-  error "right: TODO; this case has to be handled separately."-right _ = error "left: Tree is not bifurcating."---- | Connect two trees with a branch in all possible ways.------ Basically, introduce a branch between two trees. If the trees have n, and m--- branches, respectively, there are n*m ways to connect them.------ A base node has to be given which will be used wherever the new node is--- introduced.-connect :: a -> Tree a -> Tree a -> [Tree a]-connect n l r = [ Node n [x, y] | x <- roots l, y <- roots r ]---- | Get clades induced by multifurcations.-clades :: Ord a => Tree a -> [S.Set a]-clades (Node _ []    ) = []-clades (Node _ [x]   ) = clades x-clades (Node _ [x, y]) = clades x ++ clades y-clades t               = S.fromList (leaves t) : concatMap clades (subForest t)
+ src/ELynx/Data/Tree/Zipper.hs view
@@ -0,0 +1,160 @@+-- |+-- 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}
src/ELynx/Distribution/BirthDeath.hs view
@@ -1,92 +1,95 @@ {-# 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.---}+{-# 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+  ( BirthDeathDistribution (..),+    cumulative,+    density,+    quantile,   ) where -import           Data.Data                      ( Data-                                                , Typeable-                                                )-import           GHC.Generics                   ( Generic )-import qualified Statistics.Distribution       as D--import           ELynx.Distribution.Types+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-  { bddTOr :: Time         -- ^ Time to origin of the tree.-  , bddLa  :: Rate         -- ^ Birth rate.-  , bddMu  :: Rate         -- ^ Death rate.-  } deriving (Eq, Typeable, Data, Generic)+  { -- | 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))+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+  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))+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 [0,1] range. 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))+  | 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
src/ELynx/Distribution/BirthDeathCritical.hs view
@@ -1,81 +1,84 @@ {-# 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.---}+{-# 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+  ( BirthDeathCriticalDistribution (..),+    cumulative,+    density,+    quantile,   ) where -import           Data.Data                      ( Data-                                                , Typeable-                                                )-import           GHC.Generics                   ( Generic )-import qualified Statistics.Distribution       as D--import           ELynx.Distribution.Types+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-  { bdcdTOr :: Time         -- ^ Time to origin of the tree.-  , bdcdLa  :: Rate    -- ^ Birth and death rate.-  } deriving (Eq, Typeable, Data, Generic)+  { -- | 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+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+  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)+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)+  | 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
src/ELynx/Distribution/BirthDeathCriticalNoTime.hs view
@@ -1,75 +1,76 @@ {-# 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.---}+{-# 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+  ( BirthDeathCriticalNoTimeDistribution (..),+    cumulative,+    density,+    quantile,   ) where -import           Data.Data                      ( Data-                                                , Typeable-                                                )-import           GHC.Generics                   ( Generic )-import qualified Statistics.Distribution       as D--import           ELynx.Distribution.Types+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-  { bdcntdLa  :: Rate    -- ^ Birth and death rate.-  } deriving (Eq, Typeable, Data, Generic)+  { -- | 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)+cumulative (BDCNTD l) x+  | x <= 0 = 0+  | otherwise = x * l / (1.0 + x * l)  instance D.ContDistr BirthDeathCriticalNoTimeDistribution where-  density  = density+  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)+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 ++ "."+  | 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
src/ELynx/Distribution/BirthDeathNearlyCritical.hs view
@@ -1,106 +1,109 @@ {-# 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).---}+{-# 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+  ( BirthDeathNearlyCriticalDistribution (..),+    cumulative,+    density,+    quantile,   ) where -import           Data.Data                      ( Data-                                                , Typeable-                                                )-import           GHC.Generics                   ( Generic )-import qualified Statistics.Distribution       as D--import           ELynx.Distribution.Types+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-  { bdncdTOr :: Time    -- ^ Time to origin of the tree.-  , bdncdLa  :: Rate    -- ^ Birth and death rate.-  , bdncdMu  :: Rate    -- ^ Birth and death rate.-  } deriving (Eq, Typeable, Data, Generic)+  { -- | 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)+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+  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)+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+  | 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-        * m-        * m-        )-  t2  = t2Nom / (den ** 2)-  res = 0.5 * (t1 + sqrt t2)+        * ( 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
src/ELynx/Distribution/CoalescentContinuous.hs view
@@ -1,28 +1,26 @@-{- |-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+-- 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+  ( coalescentDistributionCont,   ) where -import           Numeric.SpecFunctions          ( choose )-import           Statistics.Distribution.Exponential+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-  :: Int -- ^ Sample size.-  -> ExponentialDistribution+coalescentDistributionCont ::+  -- | Sample size.+  Int ->+  ExponentialDistribution coalescentDistributionCont n = exponential (choose n 2)
src/ELynx/Distribution/TimeOfOrigin.hs view
@@ -1,93 +1,95 @@ {-# 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.---}+{-# 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+  ( TimeOfOriginDistribution (..),+    cumulative,+    density,+    quantile,   ) where -import           Data.Data                      ( Data-                                                , Typeable-                                                )-import           GHC.Generics                   ( Generic )-import qualified Statistics.Distribution       as D--import           ELynx.Distribution.Types+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-  { todTN :: Int           -- ^ Number of leaves of the tree.-  , todLa :: Rate          -- ^ Birth rate.-  , todMu :: Rate          -- ^ Death rate.-  } deriving (Eq, Typeable, Data, Generic)+  { -- | 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))+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+  density = density   quantile = quantile  -- | The density function Eq. (5). density :: TimeOfOriginDistribution -> Time -> Double density (TOD nn l m) x-  | x < 0     = 0+  | 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+  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+  | 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
src/ELynx/Distribution/TimeOfOriginNearCritical.hs view
@@ -1,93 +1,96 @@ {-# 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.---}+{-# 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+  ( TimeOfOriginNearCriticalDistribution (..),+    cumulative,+    density,+    quantile,   ) where -import           Data.Data                      ( Data-                                                , Typeable-                                                )-import           GHC.Generics                   ( Generic )-import qualified Statistics.Distribution       as D--import           ELynx.Distribution.Types+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-  { todTN :: Int           -- ^ Number of leaves of the tree.-  , todLa :: Rate          -- ^ Birth rate.-  , todMu :: Rate          -- ^ Death rate.-  } deriving (Eq, Typeable, Data, Generic)+  { -- | 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))+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+  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+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+  | 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
src/ELynx/Distribution/Types.hs view
@@ -1,21 +1,17 @@-{- |-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+-- 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+  ( Time,+    Rate,   ) where @@ -24,5 +20,3 @@  -- | Birth or death rates. type Rate = Double--
src/ELynx/Export/Tree/Newick.hs view
@@ -1,54 +1,53 @@-{- |-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+-- 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+  ( toNewick,   ) where -import qualified Data.ByteString.Lazy.Builder  as L-import qualified Data.ByteString.Lazy.Char8    as L-import           Data.List                      ( intersperse )-import           Data.Tree+import Data.ByteString.Lazy.Builder (Builder)+import qualified Data.ByteString.Lazy.Builder as L+import Data.ByteString.Lazy.Char8 (ByteString)+import Data.List (intersperse)+import ELynx.Data.Tree.Named+import ELynx.Data.Tree.Phylogeny+import ELynx.Data.Tree.Rooted+import ELynx.Tools -import           ELynx.Data.Tree-import           ELynx.Tools+toNewickBuilder :: Named a => Tree Phylo a -> Builder+toNewickBuilder t = go t <> L.word8 (c2w ';')+  where+    go (Node b l []) = lbl b l+    go (Node b l ts) =+      L.word8 (c2w '(')+        <> mconcat (intersperse (L.word8 $ c2w ',') $ map go ts)+        <> L.word8 (c2w ')')+        <> lbl b l+    mBrSupBuilder x = maybe mempty (\bs -> L.word8 (c2w '[') <> L.doubleDec bs <> L.word8 (c2w ']')) (brSup x)+    mBrLenBuilder x = maybe mempty (\bl -> L.word8 (c2w ':') <> L.doubleDec bl) (brLen x)+    lbl x y =+      L.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 'L.Bytestring'. Use provided -- functions to extract node labels and branch lengths builder objects. See also -- Biobase.Newick.Export.-toNewick :: Named a => Tree (PhyloLabel a) -> L.ByteString-toNewick t = L.toLazyByteString $ go t <> L.word8 (c2w ';')- where-  go (Node l []) = lbl l-  go (Node l ts) =-    L.word8 (c2w '(')-      <> mconcat (intersperse (L.word8 $ c2w ',') $ map go ts)-      <> L.word8 (c2w ')')-      <> lbl l-  brSupStr bs = L.word8 (c2w '[') <> L.doubleDec bs <> L.word8 (c2w ']')-  mBrSup l = maybe mempty brSupStr (brSup l)-  brLenStr bl = L.word8 (c2w ':') <> L.doubleDec bl-  mBrLen l = maybe mempty brLenStr (brLen l)-  lbl l =-    L.lazyByteString (getName l)-      <> mBrLen l-             -- After reading several discussion, I go for the "more semantical-             -- form" with branch support values in square brackets.-      <> mBrSup l+toNewick :: Named a => Tree Phylo a -> ByteString+toNewick = L.toLazyByteString . toNewickBuilder
+ src/ELynx/Export/Tree/Nexus.hs view
@@ -0,0 +1,31 @@+{-# 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 Data.ByteString.Lazy (ByteString)+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 => [(ByteString, Tree Phylo a)] -> ByteString+toNexusTrees ts = toNexus "TREES" (map tree ts)++tree :: Named a => (ByteString, Tree Phylo a) -> ByteString+tree (n, t) = "  TREE " <> n <> " = " <> toNewick t
src/ELynx/Import/Tree/Newick.hs view
@@ -1,71 +1,66 @@ {-# 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.--TODO: Use 'between' for forests.---}-+-- TODO: Use 'between' for forests. +-- |+-- 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-  ( Parser-  -- * Newick tree format-  , NewickFormat(..)-  , description-  , newick-  , oneNewick-  , manyNewick-  , forest-  , leaf-  , node-  , name-  , branchLength+  ( Parser,++    -- * Newick tree format+    NewickFormat (..),+    description,+    newick,+    oneNewick,+    someNewick,   ) where -import qualified Data.ByteString.Lazy          as L-import           Data.Tree-import           Data.Void-import           Data.Word-import           Text.Megaparsec-import           Text.Megaparsec.Byte-import           Text.Megaparsec.Byte.Lexer     ( decimal-                                                , float-                                                )--import           ELynx.Data.Tree-import           ELynx.Tools+import Data.ByteString.Lazy (ByteString)+import qualified Data.ByteString.Lazy as L+import Data.Void+import Data.Word+import ELynx.Data.Tree.Measurable+import ELynx.Data.Tree.Phylogeny+import ELynx.Data.Tree.Rooted hiding (label)+import ELynx.Data.Tree.Supported+import ELynx.Tools+import Text.Megaparsec+import Text.Megaparsec.Byte+import Text.Megaparsec.Byte.Lexer+  ( decimal,+    float,+  )  -- | Shortcut.-type Parser = Parsec Void L.ByteString+type Parser = Parsec Void ByteString  -- | 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. Key value pairs are IGNORED at the moment.-data NewickFormat =-  Standard-  | IqTree-  | 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.+data NewickFormat = Standard | IqTree | RevBayes   deriving (Eq, Show, Read, Bounded, Enum, Generic)  instance FromJSON NewickFormat@@ -76,147 +71,138 @@ 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."+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 (PhyloLabel L.ByteString))+newick :: NewickFormat -> Parser (Tree Phylo ByteString) newick Standard = newickStandard-newick IqTree   = newickIqTree+newick IqTree = newickIqTree newick RevBayes = newickRevBayes  -- | Parse a single Newick tree. Fails when end of file is not reached.-oneNewick :: NewickFormat -> Parser (Tree (PhyloLabel L.ByteString))+oneNewick :: NewickFormat -> Parser (Tree Phylo ByteString) oneNewick Standard = oneNewickStandard-oneNewick IqTree   = oneNewickIqTree+oneNewick IqTree = oneNewickIqTree oneNewick RevBayes = oneNewickRevBayes --- | Parse many Newick trees until end of file.-manyNewick :: NewickFormat -> Parser [Tree (PhyloLabel L.ByteString)]-manyNewick Standard = manyNewickStandard-manyNewick IqTree   = manyNewickIqTree-manyNewick RevBayes = manyNewickRevBayes+-- | Parse one or more Newick trees until end of file.+someNewick :: NewickFormat -> Parser (Forest Phylo ByteString)+someNewick Standard = someNewickStandard+someNewick IqTree = someNewickIqTree+someNewick RevBayes = someNewickRevBayes  w :: Char -> Parser Word8 w = char . c2w --- | Parse a single Newick tree. Also succeeds when more trees follow.-newickStandard :: Parser (Tree (PhyloLabel L.ByteString))-newickStandard = space *> tree <* w ';' <* space <?> "newick"+-- Parse a single Newick tree. Also succeeds when more trees follow.+newickStandard :: Parser (Tree Phylo ByteString)+newickStandard = space *> tree <* w ';' <* space <?> "newickStandard" --- | Parse a single Newick tree. Fails when end of file is not reached.-oneNewickStandard :: Parser (Tree (PhyloLabel L.ByteString))-oneNewickStandard = newickStandard <* eof <?> "oneNewick"+-- Parse a single Newick tree. Fails when end of file is not reached.+oneNewickStandard :: Parser (Tree Phylo ByteString)+oneNewickStandard = newickStandard <* eof <?> "oneNewickStandard" --- | Parse many Newick trees until end of file.-manyNewickStandard :: Parser [Tree (PhyloLabel L.ByteString)]-manyNewickStandard = some newickStandard <* eof <?> "manyNewick"+-- Parse one ore more Newick trees until end of file.+someNewickStandard :: Parser (Forest Phylo ByteString)+someNewickStandard = some newickStandard <* eof <?> "someNewickStandard" -tree :: Parser (Tree (PhyloLabel L.ByteString))+tree :: Parser (Tree Phylo ByteString) tree = branched <|> leaf <?> "tree" -branched :: Parser (Tree (PhyloLabel L.ByteString))-branched = do-  f <- forest-  n <- node <?> "branched"-  return $ Node n f---- | A 'forest' is a set of trees separated by @,@ and enclosed by parentheses.-forest :: Parser [Tree (PhyloLabel L.ByteString)]-forest = between (w '(') (w ')') (tree `sepBy1` w ',') <?> "forest"+branched :: Parser (Tree Phylo ByteString)+branched = label "branched" $ do+  f <- forestP+  n <- name+  p <- phylo+  return $ Node p n f --- TODO: Why try?-branchSupport :: Parser (Maybe Double)-branchSupport = optional $ do-  _ <- try $ w '['-  s <- try float <|> try decimalAsDouble-  _ <- try $ w ']'-  return s+-- A 'forest' is a set of trees separated by @,@ and enclosed by parentheses.+forestP :: Parser (Forest Phylo ByteString)+forestP = between (w '(') (w ')') (tree `sepBy1` w ',') <?> "forestP" --- | A 'leaf' is a 'node' without children.-leaf :: Parser (Tree (PhyloLabel L.ByteString))-leaf = do-  n <- node <?> "leaf"-  return $ Node n []+branchSupport :: Parser (Maybe BranchSupport)+branchSupport = label "branchSupport" $+  optional $+    try $+      do+        _ <- w '['+        s <- float <|> decimalAsDouble+        _ <- w ']'+        return s --- | A 'node' has a name and a 'branchLength'.-node :: Parser (PhyloLabel L.ByteString)-node = do+-- A 'leaf' has a 'name' and a 'phylo' branch.+leaf :: Parser (Tree Phylo ByteString)+leaf = label "leaf" $ do   n <- name-  b <- branchLength-  s <- branchSupport <?> "node"-  return $ PhyloLabel n s b+  p <- phylo+  return $ Node p n []  checkNameCharacter :: Word8 -> Bool checkNameCharacter c = c `notElem` map c2w " :;()[]," --- | A name can be any string of printable characters except blanks, colons,+-- A name can be any string of printable characters except blanks, colons, -- semicolons, parentheses, and square brackets (and commas).-name :: Parser L.ByteString+name :: Parser ByteString name = L.pack <$> many (satisfy checkNameCharacter) <?> "name" --- | Branch length.-branchLength :: Parser (Maybe Double)+phylo :: Parser Phylo+phylo = Phylo <$> branchLength <*> branchSupport <?> "phylo"++-- Branch length.+branchLength :: Parser (Maybe BranchLength) branchLength = optional (w ':' *> branchLengthGiven) <?> "branchLength"  branchLengthGiven :: Parser Double-branchLengthGiven = try float <|> decimalAsDouble+branchLengthGiven = try float <|> decimalAsDouble <?> "branchLengthGiven"  decimalAsDouble :: Parser Double-decimalAsDouble = fromIntegral <$> (decimal :: Parser Int)+decimalAsDouble = fromIntegral <$> (decimal :: Parser Int) <?> "decimalAsDouble"  -------------------------------------------------------------------------------- -- 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 (PhyloLabel L.ByteString))+newickIqTree :: Parser (Tree Phylo ByteString) newickIqTree = space *> treeIqTree <* w ';' <* space <?> "newickIqTree"  -- See 'newickIqTree'. Parse a single Newick tree. Fails when end of file is not -- reached.-oneNewickIqTree :: Parser (Tree (PhyloLabel L.ByteString))+oneNewickIqTree :: Parser (Tree Phylo ByteString) oneNewickIqTree = newickIqTree <* eof <?> "oneNewickIqTree" --- See 'newickIqTree'. Parse many Newick trees until end of file.-manyNewickIqTree :: Parser [Tree (PhyloLabel L.ByteString)]-manyNewickIqTree = some newickIqTree <* eof <?> "manyNewickIqTree"+-- See 'newickIqTree'. Parse one ore more Newick trees until end of file.+someNewickIqTree :: Parser (Forest Phylo ByteString)+someNewickIqTree = some newickIqTree <* eof <?> "someNewickIqTree"  -- IQ-TREE stores the branch support as node names after the closing bracket of a forest.-treeIqTree :: Parser (Tree (PhyloLabel L.ByteString))+treeIqTree :: Parser (Tree Phylo ByteString) treeIqTree = branchedIqTree <|> leaf <?> "treeIqTree"  -- IQ-TREE stores the branch support as node names after the closing bracket of a forest.-branchedIqTree :: Parser (Tree (PhyloLabel L.ByteString))-branchedIqTree = do+branchedIqTree :: Parser (Tree Phylo ByteString)+branchedIqTree = label "branchedIqTree" $ do   f <- forestIqTree   s <- branchSupportIqTree-  n <- nodeIqTree <?> "branchedIqTree"-  let n' = n { brSup = s }-  return $ Node n' f+  n <- name+  b <- 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 [Tree (PhyloLabel L.ByteString)]-forestIqTree = do+forestIqTree :: Parser (Forest Phylo ByteString)+forestIqTree = label "forestIqTree" $ do   _ <- w '('   f <- treeIqTree `sepBy1` w ','-  _ <- w ')' <?> "forestIqTree"+  _ <- w ')'   return f  -- TODO: Same here, why try? -- IQ-TREE stores the branch support as node names after the closing bracket of a forest.-branchSupportIqTree :: Parser (Maybe Double)-branchSupportIqTree = optional $ try float <|> try decimalAsDouble---- IQ-TREE stores the branch support as node names after the closing bracket of a forest.-nodeIqTree :: Parser (PhyloLabel L.ByteString)-nodeIqTree = do-  n <- name-  b <- branchLength <?> "nodeIqTree"-  return $ PhyloLabel n Nothing b-+branchSupportIqTree :: Parser (Maybe BranchSupport)+branchSupportIqTree = label "branchSupportIqTree" $ optional $ try float <|> try decimalAsDouble  -------------------------------------------------------------------------------- -- RevBayes.@@ -225,53 +211,55 @@ -- 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 (PhyloLabel L.ByteString))+-- XXX: Key value pairs are ignored at the moment.+newickRevBayes :: Parser (Tree Phylo ByteString) newickRevBayes =-  space *> brackets *> treeRevBayes <* w ';' <* space <?> "newickRevBayes"+  space *> optional brackets *> treeRevBayes <* w ';' <* space <?> "newickRevBayes"  -- See 'newickRevBayes'. Parse a single Newick tree. Fails when end of file is -- not reached.-oneNewickRevBayes :: Parser (Tree (PhyloLabel L.ByteString))+oneNewickRevBayes :: Parser (Tree Phylo ByteString) oneNewickRevBayes = newickRevBayes <* eof <?> "oneNewickRevBayes" --- See 'newickRevBayes'. Parse many Newick trees until end of file.-manyNewickRevBayes :: Parser [Tree (PhyloLabel L.ByteString)]-manyNewickRevBayes = some newickRevBayes <* eof <?> "manyNewickRevBayes"+-- See 'newickRevBayes'. Parse one ore more Newick trees until end of file.+someNewickRevBayes :: Parser (Forest Phylo ByteString)+someNewickRevBayes = some newickRevBayes <* eof <?> "someNewickRevBayes" -treeRevBayes :: Parser (Tree (PhyloLabel L.ByteString))+treeRevBayes :: Parser (Tree Phylo ByteString) treeRevBayes = branchedRevBayes <|> leafRevBayes <?> "treeRevBayes" -branchedRevBayes :: Parser (Tree (PhyloLabel L.ByteString))-branchedRevBayes = do+branchedRevBayes :: Parser (Tree Phylo ByteString)+branchedRevBayes = label "branchedRevgBayes" $ do   f <- forestRevBayes-  n <- nodeRevBayes <?> "branchedRevBayes"-  return $ Node n f+  n <- nameRevBayes+  b <- branchLengthRevBayes+  return $ Node (Phylo b Nothing) n f -forestRevBayes :: Parser [Tree (PhyloLabel L.ByteString)]-forestRevBayes = do+forestRevBayes :: Parser (Forest Phylo ByteString)+forestRevBayes = label "forestRevBayes" $ do   _ <- w '('   f <- treeRevBayes `sepBy1` w ','-  _ <- w ')' <?> "forestRevBayes"+  _ <- w ')'   return f -nodeRevBayes :: Parser (PhyloLabel L.ByteString)-nodeRevBayes = do-  n <- name-  _ <- optional brackets-  b <- branchLength-  _ <- optional brackets <?> "nodeRevBayes"-  return $ PhyloLabel n Nothing b+nameRevBayes :: Parser ByteString+nameRevBayes = name <* optional brackets <?> "nameRevBayes" -leafRevBayes :: Parser (Tree (PhyloLabel L.ByteString))-leafRevBayes = do-  n <- nodeRevBayes <?> "leafRevBayes"-  return $ Node n []+branchLengthRevBayes :: Parser (Maybe BranchLength)+branchLengthRevBayes = branchLength <* optional brackets <?> "branchLengthRevBayes" +leafRevBayes :: Parser (Tree Phylo ByteString)+leafRevBayes = label "leafRevBayes" $ do+  n <- nameRevBayes+  b <- branchLengthRevBayes+  return $ Node (Phylo b Nothing) n []+ -- Drop anything between brackets. brackets :: Parser ()-brackets = do-  _ <- between (w '[')-               (w ']')-               (takeWhileP (Just "allCharsButBracketEnd") (/= c2w ']'))+brackets = label "brackets" $ do+  _ <-+    between+      (w '[')+      (w ']')+      (takeWhileP (Just "allCharsButBracketEnd") (/= c2w ']'))   return ()
+ src/ELynx/Import/Tree/Nexus.hs view
@@ -0,0 +1,45 @@+{-# 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 Data.ByteString.Internal (c2w)+import Data.ByteString.Lazy (ByteString, pack)+import ELynx.Data.Tree.Phylogeny+import ELynx.Data.Tree.Rooted+import ELynx.Import.Nexus hiding (Parser)+import ELynx.Import.Tree.Newick+import Text.Megaparsec+import Text.Megaparsec.Byte++-- | Parse a Nexus files with a TREES block.+nexusTrees :: NewickFormat -> Parser [(ByteString, Tree Phylo ByteString)]+nexusTrees = nexus . trees++trees :: NewickFormat -> Block [(ByteString, Tree Phylo ByteString)]+trees f = Block "TREES" (some $ namedNewick f)++namedNewick :: NewickFormat -> Parser (ByteString, Tree Phylo ByteString)+namedNewick f = do+  _ <- space+  _ <- string "TREE"+  _ <- space+  n <- some alphaNumChar+  _ <- space+  _ <- char $ c2w '='+  _ <- space+  t <- newick f+  return (pack n, t)
src/ELynx/Simulate/Coalescent.hs view
@@ -1,62 +1,60 @@-{- |-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+-- 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+  ( simulate,   ) where -import           Control.Monad.Primitive-import           Data.Tree-import           Statistics.Distribution-import           System.Random.MWC--import           ELynx.Data.Tree-import           ELynx.Distribution.CoalescentContinuous+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)-  => Int -- ^ Number of leaves.-  -> Gen (PrimState m)-  -> m (Tree (PhyloLabel Int))+simulate ::+  (PrimMonad m) =>+  -- | Number of leaves.+  Int ->+  Gen (PrimState m) ->+  m (Tree Length Int) simulate n = simulate' n 0 trs-  where trs = [ singleton (PhyloLabel i Nothing Nothing) | i <- [0 .. n - 1] ]+  where+    trs = [Node (Length 0) i [] | i <- [0 .. n - 1]] -simulate'-  :: (PrimMonad m)-  => Int-  -> Int-  -> [Tree (PhyloLabel Int)]-  -> Gen (PrimState m)-  -> m (Tree (PhyloLabel Int))+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.+    -- 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 (lengthenStem t) trs -- Move time 't' up on the tree.-        tl    = trs' !! (i - 1)-        tr    = trs' !! i+    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 (PhyloLabel a Nothing Nothing) [tl, tr]+        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
src/ELynx/Simulate/PointProcess.hs view
@@ -1,52 +1,52 @@ {-# 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+--   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+  ( PointProcess (..),+    TimeSpec,+    simulate,+    toReconstructedTree,+    simulateReconstructedTree,+    simulateNReconstructedTrees,   ) where -import           Control.Monad-import           Control.Monad.Primitive-import           Data.List                      ( mapAccumL )-import           Data.Tree-import qualified Statistics.Distribution       as D-                                                ( genContVar )-import           System.Random.MWC--import           ELynx.Data.Tree-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           ELynx.Tools+import Control.Monad+import Control.Monad.Primitive+import Data.List (mapAccumL)+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 ELynx.Tools+import qualified Statistics.Distribution as D+  ( genContVar,+  )+import System.Random.MWC  epsNearCriticalPointProcess :: Double epsNearCriticalPointProcess = 1e-5@@ -61,9 +61,11 @@ -- 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)+  { 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@@ -72,57 +74,57 @@  -- | Sample a point process using the 'BirthDeathDistribution'. The names of the -- points will be integers.-simulate-  :: (PrimMonad m)-  => Int        -- ^ Number of points (samples)-  -> TimeSpec   -- ^ Time of origin or MRCA-  -> Rate       -- ^ Birth rate-  -> Rate       -- ^ Death rate-  -> Gen (PrimState m)   -- ^ Generator (see 'System.Random.MWC')-  -> m (PointProcess Int Double)+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+  | -- 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 :).+  | -- 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@@ -150,45 +152,58 @@ -- that they will have while creating the tree). sort :: (Ord b) => PointProcess a b -> ([b], [Int]) sort (PointProcess _ vs _) = (vsSorted, isSorted)- where-  vsIsSorted = sortListWithIndices vs-  vsSorted   = map fst vsIsSorted-  isSorted   = flattenIndices $ map snd vsIsSorted+  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)-  => Int        -- ^ Number of trees-  -> Int        -- ^ Number of points (samples)-  -> TimeSpec   -- ^ Time of origin or MRCA-  -> Rate       -- ^ Birth rate-  -> Rate       -- ^ Death rate-  -> Gen (PrimState m)   -- ^ Generator (see 'System.Random.MWC')-  -> m [Tree (PhyloLabel Int)]+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 []+  | 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)-  => Int        -- ^ Number of points (samples)-  -> TimeSpec   -- ^ Time of origin or MRCA-  -> Rate       -- ^ Birth rate-  -> Rate       -- ^ Death rate-  -> Gen (PrimState m)   -- ^ Generator (see 'System.Random.MWC')-  -> m (Tree (PhyloLabel Int))+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 @@ -202,50 +217,51 @@ -- 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 (PhyloLabel a)+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."-  |-  -- -- XXX: Test is deactivated.-  -- -- | otherwise = if isReconstructed treeOrigin then treeOrigin else error "Error in algorithm."-    otherwise                  = treeOrigin- where-  (vsSorted, isSorted) = sort pp-  !lvs                 = [ singleton (PhyloLabel p Nothing Nothing) | p <- ps ]-  !heights             = replicate (length ps) 0-  !treeRoot            = toReconstructedTree' isSorted vsSorted l lvs heights-  !h                   = last vsSorted-  !treeOrigin          = lengthenStem (o - h) treeRoot+  | 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) = sort 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.-  -> [Tree (PhyloLabel a)] -- Leaves with accumulated root branch lengths.-  -> [Double]              -- Accumulated heights of the leaves.-  -> Tree (PhyloLabel a)-toReconstructedTree' [] [] _ trs _  = head trs+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'-  -- For the algorithm, see 'ELynx.Coalescent.simulate', but index starts-  -- at zero.- where-  !i     = head is-  !is'   = tail is-  !v     = head vs-  !vs'   = tail vs-  -- Left: l, right: r.-  !hl    = hs !! i-  !hr    = hs !! (i + 1)-  !dvl   = v - hl-  !dvr   = v - hr-  !tl    = lengthenStem dvl $ trs !! i-  !tr    = lengthenStem dvr $ trs !! (i + 1)-  !h'    = hl + dvl       -- Should be the same as 'hr + dvr'.-  !tm    = Node (PhyloLabel l Nothing Nothing) [tl, tr]-  !trs'' = take i trs ++ [tm] ++ drop (i + 2) trs-  !hs'   = take i hs ++ [h'] ++ drop (i + 2) 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
+ test/ELynx/Data/Topology/RootedSpec.hs view
@@ -0,0 +1,23 @@+-- |+-- Module      :  ELynx.Data.Topology.RootedSpec+-- Description :  Unit tests for ELynx.Data.Topology.Rooted+-- 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:58:16 2020.+module ELynx.Data.Topology.RootedSpec+  ( spec,+  )+where++import Test.Hspec++spec :: Spec+spec =+  describe "TODO" $+    it "returns the original number when given a positive input" $+      (1 :: Int) `shouldBe` 1
+ test/ELynx/Data/Tree/Arbitrary.hs view
@@ -0,0 +1,66 @@+{-# OPTIONS_GHC -Wno-orphans #-}++-- |+-- Module      :  ELynx.Data.Tree.Arbitrary+-- Description :  Arbitrary instance, needed for QuickCheck+-- Copyright   :  (c) Dominik Schrempf 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Tue Apr 21 17:13:23 2020.+module ELynx.Data.Tree.Arbitrary+  (+  )+where++import Data.Traversable+import ELynx.Data.Tree+import Test.QuickCheck++-- Of course, the boundaries for branch support and length have been chosen+-- pretty arbitrarily :).+--+-- XXX: This instance does not produce values without branch lengths nor branch+-- supports.+instance Arbitrary Phylo where+  arbitrary = Phylo <$> (Just <$> choose (1, 100)) <*> (Just <$> choose (0.001, 10))++instance Arbitrary2 Tree where+  liftArbitrary2 arbB arbN = go+    where+      go = sized $ \n -> do+        -- Sized is the size of the trees.+        br <- arbB+        val <- arbN+        pars <- frequency [(1, pure [1, 1]), (3, arbPartition (n - 1))] -- can go negative!+        frst <- for pars $ \i -> resize i go+        return $ Node br val frst+      arbPartition :: Int -> Gen [Int]+      arbPartition k = case compare k 1 of+        LT -> pure []+        EQ -> pure [1]+        GT -> do+          first <- elements [1 .. k]+          rest <- arbPartition $ k - first+          return $ first : rest++  liftShrink2 _ shrN = go+    where+      go (Node br val frst) =+        frst+          ++ [ Node br e fs+               | (e, fs) <- liftShrink2 shrN (liftShrink go) (val, frst)+             ]++instance (Arbitrary e, Arbitrary a) => Arbitrary (Tree e a) where+  arbitrary = arbitrary2++instance (CoArbitrary e, CoArbitrary a) => CoArbitrary (Tree e a) where+  coarbitrary (Node br val frst) =+    coarbitrary br . coarbitrary val . coarbitrary frst++instance Arbitrary Length where+  arbitrary = Length . getPositive <$> arbitrary
test/ELynx/Data/Tree/BipartitionSpec.hs view
@@ -1,89 +1,88 @@ {-# LANGUAGE OverloadedStrings #-} -{- |-Module      :  ELynx.Data.Tree.BipartitionSpec-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 09:38:50 2019.---}-+-- |+-- Module      :  ELynx.Data.Tree.BipartitionSpec+-- 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 09:38:50 2019. module ELynx.Data.Tree.BipartitionSpec-  ( spec+  ( spec,   ) where -import qualified Data.ByteString.Lazy          as L-import qualified Data.Map                      as M-import           Data.Monoid-import qualified Data.Set                      as S-import           Data.Tree-import           Test.Hspec--import           ELynx.Data.Tree-import           ELynx.Import.Tree.Newick-import           ELynx.Tools+import Data.ByteString.Lazy (ByteString)+import Data.Map (Map)+import qualified Data.Map as M+import Data.Set (Set)+import qualified Data.Set as S+import ELynx.Data.Tree+import ELynx.Import.Tree.Newick+import ELynx.Tools+import Test.Hspec -sfrom :: [L.ByteString] -> S.Set L.ByteString+sfrom :: [ByteString] -> S.Set ByteString sfrom = S.fromList  treeFileSimple :: FilePath treeFileSimple = "data/TreeDist.trees" -getSimpleTrees :: IO [Tree (PhyloLabel L.ByteString)]-getSimpleTrees = parseFileWith (manyNewick Standard) treeFileSimple+getSimpleTrees :: IO (Forest Phylo ByteString)+getSimpleTrees = parseFileWith (someNewick Standard) treeFileSimple -bipartitionToBranchAnswer :: M.Map (Bipartition L.ByteString) (Sum Double)-bipartitionToBranchAnswer = M.fromList-  [ (bp (sfrom ["B"]) (sfrom ["A", "C", "D", "E"]), Sum { getSum = 0.3 })-  , (bp (sfrom ["B", "C", "D", "E"]) (sfrom ["A"]), Sum { getSum = 0.1 })-  , (bp (sfrom ["B", "C", "E"]) (sfrom ["A", "D"]), Sum { getSum = 5.0e-2 })-  , (bp (sfrom ["B", "E"]) (sfrom ["A", "C", "D"]), Sum { getSum = 0.4 })-  , (bp (sfrom ["C"]) (sfrom ["A", "B", "D", "E"]), Sum { getSum = 1.0e-2 })-  , (bp (sfrom ["D"]) (sfrom ["A", "B", "C", "E"]), Sum { getSum = 0.25 })-  , (bp (sfrom ["E"]) (sfrom ["A", "B", "C", "D"]), Sum { getSum = 0.8 })-  ]+bipartitionToBranchAnswer :: Map (Bipartition ByteString) Length+bipartitionToBranchAnswer =+  M.fromList+    [ (bpUnsafe (sfrom ["B"]) (sfrom ["A", "C", "D", "E"]), 0.3),+      (bpUnsafe (sfrom ["B", "C", "D", "E"]) (sfrom ["A"]), 0.1),+      (bpUnsafe (sfrom ["B", "C", "E"]) (sfrom ["A", "D"]), 5.0e-2),+      (bpUnsafe (sfrom ["B", "E"]) (sfrom ["A", "C", "D"]), 0.4),+      (bpUnsafe (sfrom ["C"]) (sfrom ["A", "B", "D", "E"]), 1.0e-2),+      (bpUnsafe (sfrom ["D"]) (sfrom ["A", "B", "C", "E"]), 0.25),+      (bpUnsafe (sfrom ["E"]) (sfrom ["A", "B", "C", "D"]), 0.8)+    ] -bipartitionsFirstTree :: S.Set (Bipartition L.ByteString)-bipartitionsFirstTree = S.fromList-  [ bp (sfrom ["B"])                (sfrom ["A", "C", "D", "E"])-  , bp (sfrom ["B", "C", "D", "E"]) (sfrom ["A"])-  , bp (sfrom ["B", "D", "E"])      (sfrom ["A", "C"])-  , bp (sfrom ["B", "E"])           (sfrom ["A", "C", "D"])-  , bp (sfrom ["C"])                (sfrom ["A", "B", "D", "E"])-  , bp (sfrom ["D"])                (sfrom ["A", "B", "C", "E"])-  , bp (sfrom ["E"])                (sfrom ["A", "B", "C", "D"])-  ]+bipartitionsFirstTree :: Set (Bipartition ByteString)+bipartitionsFirstTree =+  S.fromList+    [ bpUnsafe (sfrom ["B"]) (sfrom ["A", "C", "D", "E"]),+      bpUnsafe (sfrom ["B", "C", "D", "E"]) (sfrom ["A"]),+      bpUnsafe (sfrom ["B", "D", "E"]) (sfrom ["A", "C"]),+      bpUnsafe (sfrom ["B", "E"]) (sfrom ["A", "C", "D"]),+      bpUnsafe (sfrom ["C"]) (sfrom ["A", "B", "D", "E"]),+      bpUnsafe (sfrom ["D"]) (sfrom ["A", "B", "C", "E"]),+      bpUnsafe (sfrom ["E"]) (sfrom ["A", "B", "C", "D"])+    ] -bipartitionsSecondTree :: S.Set (Bipartition L.ByteString)-bipartitionsSecondTree = S.fromList-  [ bp (sfrom ["B"])                (sfrom ["A", "C", "D", "E"])-  , bp (sfrom ["B", "C", "D", "E"]) (sfrom ["A"])-  , bp (sfrom ["B", "C", "E"])      (sfrom ["A", "D"])-  , bp (sfrom ["B", "E"])           (sfrom ["A", "C", "D"])-  , bp (sfrom ["C"])                (sfrom ["A", "B", "D", "E"])-  , bp (sfrom ["D"])                (sfrom ["A", "B", "C", "E"])-  , bp (sfrom ["E"])                (sfrom ["A", "B", "C", "D"])-  ]+bipartitionsSecondTree :: Set (Bipartition ByteString)+bipartitionsSecondTree =+  S.fromList+    [ bpUnsafe (sfrom ["B"]) (sfrom ["A", "C", "D", "E"]),+      bpUnsafe (sfrom ["B", "C", "D", "E"]) (sfrom ["A"]),+      bpUnsafe (sfrom ["B", "C", "E"]) (sfrom ["A", "D"]),+      bpUnsafe (sfrom ["B", "E"]) (sfrom ["A", "C", "D"]),+      bpUnsafe (sfrom ["C"]) (sfrom ["A", "B", "D", "E"]),+      bpUnsafe (sfrom ["D"]) (sfrom ["A", "B", "C", "E"]),+      bpUnsafe (sfrom ["E"]) (sfrom ["A", "B", "C", "D"])+    ]  spec :: Spec spec = do-  describe "bipartitions"-    $ it "calculates correct bipartitions for sample trees"-    $ do-        simpleTrees <- map removeBrInfo <$> getSimpleTrees+  describe "bipartitions" $+    it "calculates correct bipartitions for sample trees" $+      do+        simpleTrees <- getSimpleTrees         let t1 = head simpleTrees             t2 = simpleTrees !! 1-        bipartitions t1 `shouldBe` bipartitionsFirstTree-        bipartitions t2 `shouldBe` bipartitionsSecondTree-  describe "bipartitionToBranch"-    $ it "creates a map from bipartitions to branch lengths"-    $ do+        bipartitions t1 `shouldBe` Right bipartitionsFirstTree+        bipartitions t2 `shouldBe` Right bipartitionsSecondTree+  describe "bipartitionToBranch" $+    it "creates a map from bipartitions to branch lengths" $+      do         simpleTrees <- getSimpleTrees-        bipartitionToBranchLength label (Sum . getLen) (simpleTrees !! 2)-          `shouldBe` bipartitionToBranchAnswer+        (phyloToLengthTree (simpleTrees !! 2) >>= bipartitionToBranch)+          `shouldBe` Right bipartitionToBranchAnswer
test/ELynx/Data/Tree/DistanceSpec.hs view
@@ -1,257 +1,262 @@ {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE OverloadedStrings #-} -{- |-Module      :  ELynx.Data.Tree.DistanceSpec-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 09:38:50 2019.---}-+-- |+-- Module      :  ELynx.Data.Tree.DistanceSpec+-- 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 09:38:50 2019. module ELynx.Data.Tree.DistanceSpec-  ( spec+  ( spec,   ) where -import qualified Data.ByteString.Lazy.Char8    as L-import           Data.Tree-import           Test.Hspec-import           Test.QuickCheck-import           Test.QuickCheck.Instances.Containers-                                                ( )-import           ELynx.Data.Tree-import           ELynx.Data.Tree.PhyloTreeArbitraryInstance-                                                ( )-import           ELynx.Import.Tree.Newick-import           ELynx.Tools+-- import Data.Bifunctor+import Data.ByteString.Lazy.Char8 (ByteString)+import ELynx.Data.Tree+-- import qualified Data.Set as S+import ELynx.Data.Tree.Arbitrary ()+import ELynx.Import.Tree.Newick+-- import ELynx.Export.Tree.Newick+import ELynx.Tools+import Test.Hspec+import Test.QuickCheck  treeFileSimple :: FilePath treeFileSimple = "data/TreeDist.trees" -getSimpleTrees :: IO [Tree (PhyloLabel L.ByteString)]-getSimpleTrees = parseFileWith (manyNewick Standard) treeFileSimple+getSimpleTrees :: IO (Forest Phylo ByteString)+getSimpleTrees = parseFileWith (someNewick Standard) treeFileSimple  treeFileMany :: FilePath treeFileMany = "data/Many.trees" -getManyTrees :: IO [Tree (PhyloLabel L.ByteString)]-getManyTrees = parseFileWith (manyNewick Standard) treeFileMany+getManyTrees :: IO (Forest Phylo ByteString)+getManyTrees = parseFileWith (someNewick Standard) treeFileMany  -- I used treedist from Phylip to get the correct results. -- See http://evolution.genetics.washington.edu/phylip/doc/treedist.html. symmetricAnswers :: [Int] symmetricAnswers =-  [ 6-  , 8-  , 0-  , 0-  , 12-  , 20-  , 18-  , 20-  , 10-  , 2-  , 10-  , 4-  , 4-  , 4-  , 4-  , 4-  , 10-  , 16-  , 8-  , 2-  , 4-  , 0-  , 0-  , 0-  , 10-  , 4-  , 0-  , 0-  , 2-  , 2-  , 0-  , 0-  , 4-  , 0-  , 2-  , 0-  , 8-  , 6-  , 2-  , 6-  , 4-  , 4-  , 8-  , 0-  , 0-  , 4-  , 2-  , 0-  , 10-  , 0-  , 0-  , 10+  [ 6,+    8,+    0,+    0,+    12,+    20,+    18,+    20,+    10,+    2,+    10,+    4,+    4,+    4,+    4,+    4,+    10,+    16,+    8,+    2,+    4,+    0,+    0,+    0,+    10,+    4,+    0,+    0,+    2,+    2,+    0,+    0,+    4,+    0,+    2,+    0,+    8,+    6,+    2,+    6,+    4,+    4,+    8,+    0,+    0,+    4,+    2,+    0,+    10,+    0,+    0,+    10   ]  branchScoreAnswers :: [Double] branchScoreAnswers =-  [ 8.567916e-02-  , 9.570577e-02-  , 1.704571e-02-  , 7.603990e-03-  , 6.149761e-01-  , 3.557070e-01-  , 2.329811e-01-  , 3.820208e-01-  , 1.895421e-02-  , 6.302364e-03-  , 2.083286e-02-  , 1.023777e-03-  , 2.138244e-02-  , 1.444380e-02-  , 1.958628e-02-  , 6.089461e-03-  , 2.551873e-02-  , 8.041220e-02-  , 4.123102e-02-  , 8.241811e-03-  , 2.623805e-02-  , 2.109278e-02-  , 1.953769e-02-  , 4.459926e-03-  , 6.594537e-02-  , 7.040703e-02-  , 8.603133e-03-  , 3.878009e-03-  , 2.969969e-02-  , 2.505262e-02-  , 2.095988e-02-  , 8.461041e-03-  , 5.228005e-02-  , 6.001320e-02-  , 8.276652e-03-  , 6.966115e-03-  , 7.701581e-02-  , 4.946339e-02-  , 2.548024e-02-  , 5.800598e-03-  , 3.875927e-02-  , 2.836737e-02-  , 9.059706e-02-  , 1.333325e-02-  , 5.071356e-02-  , 7.433056e-02-  , 3.854717e-02-  , 3.255993e-02-  , 1.581909e-01-  , 6.813096e-02-  , 8.210513e-02-  , 7.664642e-02+  [ 8.567916e-02,+    9.570577e-02,+    1.704571e-02,+    7.603990e-03,+    6.149761e-01,+    3.557070e-01,+    2.329811e-01,+    3.820208e-01,+    1.895421e-02,+    6.302364e-03,+    2.083286e-02,+    1.023777e-03,+    2.138244e-02,+    1.444380e-02,+    1.958628e-02,+    6.089461e-03,+    2.551873e-02,+    8.041220e-02,+    4.123102e-02,+    8.241811e-03,+    2.623805e-02,+    2.109278e-02,+    1.953769e-02,+    4.459926e-03,+    6.594537e-02,+    7.040703e-02,+    8.603133e-03,+    3.878009e-03,+    2.969969e-02,+    2.505262e-02,+    2.095988e-02,+    8.461041e-03,+    5.228005e-02,+    6.001320e-02,+    8.276652e-03,+    6.966115e-03,+    7.701581e-02,+    4.946339e-02,+    2.548024e-02,+    5.800598e-03,+    3.875927e-02,+    2.836737e-02,+    9.059706e-02,+    1.333325e-02,+    5.071356e-02,+    7.433056e-02,+    3.854717e-02,+    3.255993e-02,+    1.581909e-01,+    6.813096e-02,+    8.210513e-02,+    7.664642e-02   ] -prop_dist_same_tree-  :: (Num b, Eq b) => (Tree a -> Tree a -> b) -> Tree a -> Bool-prop_dist_same_tree distanceMeasure t = distanceMeasure t t == 0+prop_dist_same_tree :: (Num b, Eq b) => (Tree e a -> Tree e a -> Either String b) -> Tree e a -> Bool+prop_dist_same_tree distanceMeasure t = distanceMeasure t t == Right 0  each :: Int -> [a] -> [a] each n = map head . takeWhile (not . null) . iterate (drop n) -multifurcating :: Tree Char-multifurcating = Node-  ' '-  [Node 'A' [], Node 'B' [], Node ' ' [Node 'C' [], Node 'D' [], Node 'E' []]]+multifurcating :: Tree () Char+multifurcating =+  Node+    ()+    ' '+    [Node () 'A' [], Node () 'B' [], Node () ' ' [Node () 'C' [], Node () 'D' [], Node () 'E' []]] -bifurcatingComp :: Tree Char-bifurcatingComp = Node-  ' '-  [ Node ' ' [Node 'A' [], Node 'B' []]-  , Node ' ' [Node 'C' [], Node ' ' [Node 'D' [], Node 'E' []]]-  ]+bifurcatingComp :: Tree () Char+bifurcatingComp =+  Node+    ()+    ' '+    [ Node () ' ' [Node () 'A' [], Node () 'B' []],+      Node () ' ' [Node () 'C' [], Node () ' ' [Node () 'D' [], Node () 'E' []]]+    ] -bifurcatingIncomp :: Tree Char-bifurcatingIncomp = Node-  ' '-  [ Node ' ' [Node 'A' [], Node 'C' []]-  , Node ' ' [Node 'B' [], Node ' ' [Node 'D' [], Node 'E' []]]-  ]+bifurcatingIncomp :: Tree () Char+bifurcatingIncomp =+  Node+    ()+    ' '+    [ Node () ' ' [Node () 'A' [], Node () 'C' []],+      Node () ' ' [Node () 'B' [], Node () ' ' [Node () 'D' [], Node () 'E' []]]+    ] -incSplitTree1a :: Tree (PhyloLabel L.ByteString)+incSplitTree1a :: Tree Phylo ByteString incSplitTree1a = parseByteStringWith "" (oneNewick IqTree) "((a,b)0.7,(c,d));" -incSplitTree1b :: Tree (PhyloLabel L.ByteString)+incSplitTree1b :: Tree Phylo ByteString incSplitTree1b = parseByteStringWith "" (oneNewick IqTree) "((a,b)0.7,c,d);" -incSplitTree2 :: Tree (PhyloLabel L.ByteString)+incSplitTree2 :: Tree Phylo ByteString incSplitTree2 = parseByteStringWith "" (oneNewick IqTree) "((a,c),(b,d));" -incSplitTree3 :: Tree (PhyloLabel L.ByteString)-incSplitTree3 =-  parseByteStringWith "" (oneNewick IqTree) "(((a,b)0.7,c),(d,e));"+incSplitTree3 :: Tree Phylo ByteString+incSplitTree3 = parseByteStringWith "" (oneNewick IqTree) "(((a,b)0.7,c),(d,e));" -incSplitTree4 :: Tree (PhyloLabel L.ByteString)+incSplitTree4 :: Tree Phylo ByteString incSplitTree4 = parseByteStringWith "" (oneNewick IqTree) "(((a,c),b),(d,e));" +-- Compute distances between adjacent pairs of a list of input trees. Use given+-- distance measure.+adjacent :: (a -> a -> b) -> [a] -> [b]+adjacent dist trs = [dist x y | (x, y) <- zip trs (tail trs)]++-- noPL :: Phylo+-- noPL = Phylo Nothing Nothing+ spec :: Spec spec = do   describe "symmetric" $ do     it "calculates correct distances for sample trees" $ do       simpleTrees <- getSimpleTrees-      symmetric (head simpleTrees) (simpleTrees !! 1) `shouldBe` 2+      symmetric (head simpleTrees) (simpleTrees !! 1) `shouldBe` Right 2       manyTrees <- getManyTrees       -- Since treedist computes the distance between adjacent pairs, in the       -- following manner: [tr0, tr1, tr2, tr3] -> [dist tr0 tr1, dist tr2 tr3],       -- we have to skip some distances.-      each 2 (adjacent (symmetricWith getName) manyTrees)-        `shouldBe` symmetricAnswers-    it "is zero for a collection of random trees"-      $ property-      $ prop_dist_same_tree-          (symmetric :: Tree (PhyloLabel Double)-            -> Tree (PhyloLabel Double)-            -> Int-          )-+      each 2 (adjacent symmetric manyTrees)+        `shouldBe` map Right symmetricAnswers+    it "is zero for a collection of random trees" $+      property $+        prop_dist_same_tree+          (symmetric :: Tree Phylo Double -> Tree Phylo Double -> Either String Int)   describe "incompatibleSplit" $ do     it "calculates correct distances for sample trees" $ do-      incompatibleSplits multifurcating bifurcatingComp `shouldBe` 0-      incompatibleSplits bifurcatingComp multifurcating `shouldBe` 0-      incompatibleSplits bifurcatingIncomp multifurcating `shouldBe` 2-      incompatibleSplits multifurcating bifurcatingIncomp `shouldBe` 2+      incompatibleSplits multifurcating bifurcatingComp `shouldBe` Right 0+      incompatibleSplits bifurcatingComp multifurcating `shouldBe` Right 0+      -- print $ S.map bpHuman <$> bipartitions bifurcatingIncomp+      -- print $ S.map bpHuman <$> bipartitions multifurcating+      -- print $ S.map mpHuman <$> partitions bifurcatingIncomp+      -- print $ S.map mpHuman <$> partitions multifurcating+      -- print $ toNewick $ first (const noPL) bifurcatingIncomp+      -- print $ toNewick $ first (const noPL) multifurcating+      incompatibleSplits bifurcatingIncomp multifurcating `shouldBe` Right 2+      incompatibleSplits multifurcating bifurcatingIncomp `shouldBe` Right 2     it "calculates correct distances for sample trees with branch support" $ do-      incompatibleSplits incSplitTree1a incSplitTree2 `shouldBe` 2-      incompatibleSplits incSplitTree1b incSplitTree2 `shouldBe` 2-      incompatibleSplits (collapse 0.71 incSplitTree1a) incSplitTree2-        `shouldBe` 2-      incompatibleSplits (collapse 0.71 incSplitTree1b) incSplitTree2-        `shouldBe` 0-      incompatibleSplits (collapse 0.71 incSplitTree3) incSplitTree4-        `shouldBe` 0-    it "is zero for a collection of random trees"-      $ property-      $ prop_dist_same_tree-          (incompatibleSplits :: Tree (PhyloLabel Double)-            -> Tree (PhyloLabel Double)-            -> Int-          )-+      incompatibleSplits incSplitTree1a incSplitTree2 `shouldBe` Right 2+      incompatibleSplits incSplitTree1b incSplitTree2 `shouldBe` Right 2+      let t1a = phyloToSupportTreeUnsafe incSplitTree1a+          t1b = phyloToSupportTreeUnsafe incSplitTree1b+          tr2 = phyloToSupportTreeUnsafe incSplitTree2+          tr3 = phyloToSupportTreeUnsafe incSplitTree3+          tr4 = phyloToSupportTreeUnsafe incSplitTree4+      incompatibleSplits (collapse 0.7 t1a) tr2 `shouldBe` Right 2+      incompatibleSplits (collapse 0.71 t1b) tr2 `shouldBe` Right 0+      incompatibleSplits (collapse 0.71 tr3) tr4 `shouldBe` Right 0+    it "is zero for a collection of random trees" $+      property $+        prop_dist_same_tree+          (incompatibleSplits :: Tree Phylo Double -> Tree Phylo Double -> Either String Int)   describe "branchScore" $ do     it "calculates correct distances for sample trees" $ do       manyTrees <- getManyTrees-      -- print branchScoreAnswers-      each 2 (adjacent branchScore manyTrees)-        `shouldSatisfy` nearlyEqListWith 1e-5 branchScoreAnswers-    it "is zero for a collection of random trees"-      $ property-      $ prop_dist_same_tree-          (branchScore :: Tree (PhyloLabel Double)-            -> Tree (PhyloLabel Double)-            -> Double-          )+      let ts = map (either error id . phyloToLengthTree) manyTrees+      let ds = map (either error id) $ each 2 $ adjacent branchScore ts+      ds `shouldSatisfy` nearlyEqListWith 1e-5 branchScoreAnswers+    it "is zero for a collection of random trees" $+      property $+        prop_dist_same_tree+          (branchScore :: Tree Length Double -> Tree Length Double -> Either String Double)
− test/ELynx/Data/Tree/MultipartitionSpec.hs 
@@ -1,63 +0,0 @@-{-# LANGUAGE OverloadedStrings #-}--{- |-Module      :  ELynx.Data.Tree.MultipartitionSpec-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 09:38:50 2019.---}--module ELynx.Data.Tree.MultipartitionSpec-  ( spec-  )-where--import           Data.Set                       ( Set-                                                , empty-                                                , fromList-                                                )-import           Data.Tree-import           Test.Hspec-import           Test.QuickCheck-import           Test.QuickCheck.Instances.Containers-                                                ( )--import           ELynx.Data.Tree-import           ELynx.Data.Tree.PhyloTreeArbitraryInstance-                                                ( )--ex1 :: Tree Int-ex1 = Node 0 [Node 1 [], Node 2 [Node 4 [], Node 5 [], Node 6 []], Node 3 []]--sol1 :: Set (Multipartition Int)-sol1 = fromList-  [ mp [fromList [1], fromList [3], fromList [4, 5, 6]]-  , mp [fromList [1, 3], fromList [4], fromList [5], fromList [6]]-  ]--ex2 :: Tree Int-ex2 = Node 0 [Node 1 [], Node 2 [], Node 0 [Node 3 [], Node 4 []], Node 5 []]--sol2 :: Set (Multipartition Int)-sol2 =-  fromList [mp [fromList [1], fromList [2], fromList [3, 4], fromList [5]]]--prop_bifurcating_tree-  :: (Ord a, Measurable a, Named a, BranchSupported a) => Tree a -> Bool-prop_bifurcating_tree t = multipartitions (removeMultifurcations t) == empty--spec :: Spec-spec = describe "multipartitions" $ do-  it "calculates correct multipartitions for a sample tree" $ do-    -- t <- removeBrInfo <$> getMultifurcatingTree-    -- print $ multipartitions ex1-    multipartitions ex1 `shouldBe` sol1-    multipartitions ex2 `shouldBe` sol2-  it "is empty for a collection of random bifurcating trees"-    $ property (prop_bifurcating_tree :: Tree (PhyloLabel Double) -> Bool)
+ test/ELynx/Data/Tree/PartitionSpec.hs view
@@ -0,0 +1,53 @@+{-# LANGUAGE OverloadedStrings #-}++-- |+-- Module      :  ELynx.Data.Tree.PartitionSpec+-- 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 09:38:50 2019.+module ELynx.Data.Tree.PartitionSpec+  ( spec,+  )+where++import Data.Set (Set, fromList)+import ELynx.Data.Tree+import ELynx.Data.Tree.Arbitrary ()+import Test.Hspec++ex1 :: Tree () Int+ex1 = Node () 0 [Node () 1 [], Node () 2 [Node () 4 [], Node () 5 [], Node () 6 []], Node () 3 []]++sol1 :: Set (Partition Int)+sol1 =+  fromList+    [ mpUnsafe [fromList [1], fromList [3], fromList [4, 5, 6]],+      mpUnsafe [fromList [1, 3], fromList [4], fromList [5], fromList [6]]+    ]++ex2 :: Tree () Int+ex2 = Node () 0 [Node () 1 [], Node () 2 [], Node () 0 [Node () 3 [], Node () 4 []], Node () 5 []]++sol2 :: Set (Partition Int)+sol2 =+  fromList+    [ mpUnsafe [fromList [1], fromList [2], fromList [3, 4], fromList [5]],+      mpUnsafe [fromList [1, 2, 5], fromList [3], fromList [4]]+    ]++spec :: Spec+spec = describe "partitions" $+  it "calculates correct partitions for a sample tree" $+    do+      -- t <- removeBrInfo <$> getMultifurcatingTree+      -- print $ partitions ex1+      partitions ex1 `shouldBe` Right sol1+      partitions ex2 `shouldBe` Right sol2++-- it "is empty for a collection of random bifurcating trees"+--   $ property (prop_bifurcating_tree :: Tree (PhyloLabel Double) -> Bool)
− test/ELynx/Data/Tree/PhyloTreeArbitraryInstance.hs 
@@ -1,33 +0,0 @@-{- |-Module      :  ELynx.Data.Tree.PhyloTreeArbitraryInstance-Description :  Arbitrary instance, needed for QuickCheck-Copyright   :  (c) Dominik Schrempf 2020-License     :  GPL-3.0-or-later--Maintainer  :  dominik.schrempf@gmail.com-Stability   :  unstable-Portability :  portable--Creation date: Tue Apr 21 17:13:23 2020.---}--module ELynx.Data.Tree.PhyloTreeArbitraryInstance-  ()-where--import           Test.QuickCheck--import           ELynx.Data.Tree.PhyloTree---- Of course, the boundaries for branch support and length are chosen pretty--- arbitrarily :).------ XXX: This instance does not produce values without branch lengths nor branch--- supports.-instance Arbitrary a => Arbitrary (PhyloLabel a) where-  arbitrary =-    PhyloLabel-      <$> arbitrary-      <*> (Just <$> choose (0, 100))-      <*> (Just <$> choose (0, 10))
+ test/ELynx/Data/Tree/PhylogenySpec.hs view
@@ -0,0 +1,159 @@+{-# LANGUAGE OverloadedStrings #-}+{-# OPTIONS_GHC -fno-warn-orphans #-}++-- |+-- Module      :  ELynx.Data.Tree.PhylogenySpec+-- Description :  Unit tests for ELynx.Data.Tree.Phylogeny+-- Copyright   :  (c) Dominik Schrempf, 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Wed Jul 15 11:05:32 2020.+module ELynx.Data.Tree.PhylogenySpec+  ( spec,+  )+where++-- import Data.Bifunctor+-- import Data.ByteString.Lazy.Char8 (ByteString)+import Data.Either+import qualified Data.Set as S+import ELynx.Data.Tree+import ELynx.Data.Tree.Arbitrary ()+-- import ELynx.Import.Tree.Newick+-- import ELynx.Tools+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck hiding (labels)++simpleTree :: Tree () String+simpleTree = Node () "i" [Node () "j" [Node () "x" [], Node () "y" []], Node () "z" []]++simpleSol :: Forest () String+simpleSol =+  [ Node () "i" [Node () "j" [Node () "x" [], Node () "y" []], Node () "z" []],+    Node () "i" [Node () "j" [Node () "z" [], Node () "y" []], Node () "x" []],+    Node () "i" [Node () "j" [Node () "z" [], Node () "x" []], Node () "y" []]+  ]++instance Splittable () where+  split = id++-- Skip leaves and trees with multifurcating root nodes.+prop_roots :: Tree () a -> Bool+prop_roots t@(Node _ _ [_, _])+  | length (leaves t) == 2 = (length <$> roots t) == Right 1+  | otherwise = (length <$> roots t) == (Right $ length (labels t) - 2)+prop_roots _ = True++-- -- Skip leaves and trees with multifurcating root nodes.+-- prop_connect :: a -> Tree () a -> Tree () a -> Bool+-- prop_connect n l@(Node _ _ [_, _]) r@(Node _ _ [_, _])+--   | length (leaves l) < 3 = (length <$> connect n l r) == Right (length (flatten r) - 2)+--   | length (leaves r) < 3 = (length <$> connect n l r) == Right (length (flatten l) - 2)+--   | otherwise =+--     (length <$> connect n l r)+--       == (Right $ (length (flatten l) - 2) * (length (flatten r) - 2))+-- prop_connect _ _ _ = True++-- -- | Determine compatibility between a bipartition and a set.+-- --+-- -- If both subsets of the bipartition share elements with the given set, the+-- -- bipartition is incompatible with this subset. If all elements of the subset+-- -- are either not in the bipartition or mapping to one of the two subsets of the+-- -- bipartition, the bipartition and the subset are compatible.+-- --+-- -- See also 'ELynx.Data.Tree.Partition.compatible'.+-- bipartitionCompatible :: (Show a, Ord a) => Either String (Bipartition a) -> Set a -> Bool+-- -- compatible (Bipartition (l, r)) ss = sintersection l ss `sdisjoint` sintersection r ss+-- bipartitionCompatible (Left _) _ = False+-- bipartitionCompatible (Right p) s = S.null lOverlap || S.null rOverlap+--   where+--     (l, r) = fromBipartition p+--     lOverlap = S.intersection l s+--     rOverlap = S.intersection r s++-- compatibleAll :: (Show a, Ord a) => Tree e a -> [Set a] -> Bool+-- compatibleAll (Node _ _ [l, r]) cs =+--   all (bipartitionCompatible (bipartition l)) cs && all (bipartitionCompatible (bipartition r)) cs+-- compatibleAll _ _ = error "Tree is not bifurcating."++-- compatibleWith ::+--   (Show b, Ord b) => (a -> b) -> [Set a] -> Tree e a -> Bool+-- compatibleWith f cs t = compatibleAll (fmap f t) (map (S.map f) cs)++-- -- Get groups induced by multifurcations. Collect the leaves of all trees+-- -- induced by multifurcations.+-- multifurcatingGroups :: Tree e a -> [[a]]+-- multifurcatingGroups (Node _ _ []) = []+-- multifurcatingGroups (Node _ _ [x]) = multifurcatingGroups x+-- multifurcatingGroups (Node _ _ [x, y]) = multifurcatingGroups x ++ multifurcatingGroups y+-- multifurcatingGroups t = leaves t : concatMap multifurcatingGroups (forest t)++-- -- TODO.+-- prop_bifurcating_tree+--   :: (Ord a, Measurable a, Named a, BranchSupported 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-10) $+    map totalBranchLength $+      either error id $+        roots t+  where+    l = totalBranchLength t+prop_roots_total_length _ = True++spec :: Spec+spec = do+  -- TODO: describe "Resolve"++  describe "roots" $ do+    it "correctly handles leaves and cherries" $ do+      let tleaf = Node () 0 [] :: Tree () Int+          tcherry = Node () 0 [Node () 1 [], Node () 2 []] :: Tree () Int+      roots tleaf `shouldSatisfy` isLeft+      roots tcherry `shouldBe` Right [tcherry]+    it "correctly handles simple trees" $+      either error id (roots simpleTree) `shouldBe` simpleSol+    modifyMaxSize (* 100) $+      it "returns the correct number of rooted trees for arbitrary trees" $+        property (prop_roots :: (Tree () Int -> Bool))+  describe "rootAt" $+    modifyMaxSize (* 100) $+      it "correctly handles simple trees" $+        do+          let p = either error id $ bipartition simpleTree+          rootAt p simpleTree `shouldBe` Right simpleTree+          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))+  describe "rootsWithBranch" $+    modifyMaxSize (* 100) $+      it "does not change the tree height" $+        property (prop_roots_total_length :: Tree Length Int -> Bool)++-- -- TODO: Move this test to the executable.+-- describe "connect" $+--   modifyMaxSize (* 100) $ do+--     it "returns the correct number of rooted trees for arbitrary trees" $+--       property (prop_connect :: Int -> Tree () Int -> Tree () Int -> Bool)+--     it "correctly connects sample trees without and with constraints" $ do+--       a <- parseFileWith (oneNewick Standard) "data/ConnectA.tree"+--       b <- parseFileWith (oneNewick Standard) "data/ConnectB.tree"+--       c <- parseFileWith (someNewick Standard) "data/ConnectConstraints.tree"+--       let ts =+--             either error id $+--               connect "ROOT" (first (const ()) a) (first (const ()) b)+--           cs =+--             map S.fromList $+--               concatMap (multifurcatingGroups . first (const ())) c ::+--               [Set ByteString]+--           ts' = filter (compatibleWith getName cs) ts+--       length ts `shouldBe` 63+--       length ts' `shouldBe` 15
+ test/ELynx/Data/Tree/RootedSpec.hs view
@@ -0,0 +1,70 @@+{-# LANGUAGE OverloadedStrings #-}++-- |+-- Module      :  ELynx.Data.Tree.RootedSpec+-- Description :  Unit tests for ELynx.Data.Tree.Rooted+-- Copyright   :  (c) Dominik Schrempf 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Mon May  6 14:04:05 2019.+module ELynx.Data.Tree.RootedSpec+  ( spec,+  )+where++import Data.ByteString.Lazy.Char8 (ByteString)+import qualified Data.ByteString.Lazy.Char8 as L+import Data.Maybe+import ELynx.Data.Tree+import ELynx.Import.Tree.Newick+import ELynx.Tools+import Test.Hspec++node :: Int -> Tree () Int+node n = Node () n []++smallTree :: Tree () Int+smallTree = Node () 0 [node 1, node 2]++smallSubTree :: Tree () Int+smallSubTree = Node () 0 [node 1]++smallSubTreePruned :: Tree () Int+smallSubTreePruned = node 1++sampleTreeBS :: ByteString+sampleTreeBS =+  L.pack+    "(Aeropyrum0:0.5478645225,(((((((((Arabidopsi:0.0701001024,Oryza_sati:0.0765988261):0.0309636193,Gymnosperm:0.0520325624):0.0338982245,Physcomitr:0.0768008916):0.0895714685,(Chlamydomo:0.1136227755,Dunaliella:0.1406347323):0.1117340620):0.0818876186,Rhodophyta:0.3405656487):0.0363527066,((((((Babesia_bo:0.1646969208,Theileria0:0.1519889486):0.1908081096,Plasmodium:0.3250696762):0.0637865908,(Toxoplasma:0.1153570425,Eimeria000:0.1671916078):0.0980136930):0.0518956330,Cryptospor:0.3175062809):0.1607708388,Ciliophora:0.5687502950):0.0624078848,(Phytophtho:0.2016424948,((Thalassios:0.1202730781,Phaeodacty:0.1290341329):0.1772775509,Phaeophyce:0.1989260715):0.0312359673):0.1154768302):0.0311952864):0.0149160316,(((((((((Candida_al:0.1027755272,Saccharomy:0.1190206560):0.1333487870,Neurospora:0.1977309079):0.0522926266,Schizosacc:0.2019603227):0.0567441011,(Cryptococc:0.1948614959,Ustilago_m:0.1564451295):0.0775729694):0.0323959951,Glomus_int:0.1573670796):0.0194701292,Chytridiom:0.2228415254):0.0384370601,Encephalit:1.4622174644):0.0416231688,(((Drosophila:0.2160627753,(Mammalians:0.1080484094,Tunicates0:0.1739253014):0.0289624371):0.0346633757,Hydrozoa00:0.2058137032):0.0480963050,Monosiga_b:0.3020637584):0.0654894239):0.0380915725,(Dictyostel:0.3453588998,Mastigamoe:0.3844779231):0.0478795653):0.0129578395):1.7592083381,((Archaeoglo:0.5402784445,Methanococ:0.4088567459):0.0993669265,Pyrococcus:0.4058713829):0.1734405968):0.2193511807,Pyrobaculu:0.7507718047):0.1646616482,Sulfolobus:0.5404967897);"++largeTree :: Tree Phylo ByteString+largeTree = parseByteStringWith "Sample newick byte string" (newick Standard) sampleTreeBS++subSampleLargeTree :: Tree Phylo ByteString+subSampleLargeTree = fromJust $ dropLeavesWith ((/= 'P') . L.head) largeTree++spec :: Spec+spec = do+  describe "prune" $ do+    it "leaves a normal tree untouched" $+      prune largeTree `shouldBe` largeTree+    it "correctly prunes a small example" $+      prune smallSubTree `shouldBe` smallSubTreePruned+    it "leaves height constant for Measurable trees" $ do+      let t' =+            either error id $+              phyloToLengthTree subSampleLargeTree+      height (prune t') `shouldBe` height t'+  describe "dropLeavesWith" $ do+    it "returns the same tree if no leaves satisfy predicate" $+      dropLeavesWith (const False) smallTree `shouldBe` Just smallTree+    it "returns nothing if all leaves satisfy predicate" $+      dropLeavesWith (const True) smallTree `shouldBe` Nothing+    it "returns the correct subtree for a small example" $+      dropLeavesWith (== 2) smallTree `shouldBe` Just smallSubTree++-- TODO: intersect.
− test/ELynx/Data/Tree/TreeSpec.hs 
@@ -1,142 +0,0 @@-{-# LANGUAGE OverloadedStrings #-}--{- |-Module      :  ELynx.Data.Tree.TreeSpec-Copyright   :  (c) Dominik Schrempf 2020-License     :  GPL-3.0-or-later--Maintainer  :  dominik.schrempf@gmail.com-Stability   :  unstable-Portability :  portable--Creation date: Mon May  6 14:04:05 2019.---}--module ELynx.Data.Tree.TreeSpec-  ( spec-  )-where--import qualified Data.ByteString.Lazy.Char8    as L-import           Data.Maybe-import qualified Data.Set                      as S-import           Data.Tree-import           Test.Hspec                     ( Spec-                                                , describe-                                                , it-                                                , shouldBe-                                                )-import           Test.Hspec.QuickCheck          ( modifyMaxSize )-import           Test.QuickCheck         hiding ( label )-import           Test.QuickCheck.Instances.Containers-                                                ( )--import           ELynx.Data.Tree-import           ELynx.Import.Tree.Newick-                                         hiding ( node )-import           ELynx.Tools--node :: Int -> Tree Int-node n = Node n []--smallTree :: Tree Int-smallTree = Node 0 [node 1, node 2]--smallSubTree :: Tree Int-smallSubTree = Node 0 [node 1]--smallSubTreePruned :: Tree Int-smallSubTreePruned = node 1--sampleTreeBS :: L.ByteString-sampleTreeBS =-  L.pack-    "(Aeropyrum0:0.5478645225,(((((((((Arabidopsi:0.0701001024,Oryza_sati:0.0765988261):0.0309636193,Gymnosperm:0.0520325624):0.0338982245,Physcomitr:0.0768008916):0.0895714685,(Chlamydomo:0.1136227755,Dunaliella:0.1406347323):0.1117340620):0.0818876186,Rhodophyta:0.3405656487):0.0363527066,((((((Babesia_bo:0.1646969208,Theileria0:0.1519889486):0.1908081096,Plasmodium:0.3250696762):0.0637865908,(Toxoplasma:0.1153570425,Eimeria000:0.1671916078):0.0980136930):0.0518956330,Cryptospor:0.3175062809):0.1607708388,Ciliophora:0.5687502950):0.0624078848,(Phytophtho:0.2016424948,((Thalassios:0.1202730781,Phaeodacty:0.1290341329):0.1772775509,Phaeophyce:0.1989260715):0.0312359673):0.1154768302):0.0311952864):0.0149160316,(((((((((Candida_al:0.1027755272,Saccharomy:0.1190206560):0.1333487870,Neurospora:0.1977309079):0.0522926266,Schizosacc:0.2019603227):0.0567441011,(Cryptococc:0.1948614959,Ustilago_m:0.1564451295):0.0775729694):0.0323959951,Glomus_int:0.1573670796):0.0194701292,Chytridiom:0.2228415254):0.0384370601,Encephalit:1.4622174644):0.0416231688,(((Drosophila:0.2160627753,(Mammalians:0.1080484094,Tunicates0:0.1739253014):0.0289624371):0.0346633757,Hydrozoa00:0.2058137032):0.0480963050,Monosiga_b:0.3020637584):0.0654894239):0.0380915725,(Dictyostel:0.3453588998,Mastigamoe:0.3844779231):0.0478795653):0.0129578395):1.7592083381,((Archaeoglo:0.5402784445,Methanococ:0.4088567459):0.0993669265,Pyrococcus:0.4058713829):0.1734405968):0.2193511807,Pyrobaculu:0.7507718047):0.1646616482,Sulfolobus:0.5404967897);"--largeTree :: Tree (PhyloLabel L.ByteString)-largeTree =-  parseByteStringWith "Sample newick byte string" (newick Standard) sampleTreeBS--subSampleLargeTree :: Tree (PhyloLabel L.ByteString)-subSampleLargeTree = fromJust $ subTree ((== 'P') . L.head . label) largeTree---- XXX: Skip not bifurcating trees. This is ugly, I know.-prop_roots :: Tree a -> Bool-prop_roots t | not $ bifurcating t = True-             | length (leaves t) < 3 = length (roots t) == 1-             | otherwise = length (roots t) == 2 * length (leaves t) - 3---- XXX: Skip not bifurcating trees. This is ugly, I know.-prop_connect :: a -> Tree a -> Tree a -> Bool-prop_connect n l r-  | not (bifurcating l) || not (bifurcating r) = True-  | length (leaves l) < 3 || length (leaves r) < 3 = length (connect n l r) == 1-  | otherwise = length (connect n l r) == length (leaves l) * length (leaves r)--type Constraint a = S.Set a--compatibleAll :: (Show a, Ord a) => Tree a -> [Constraint a] -> Bool-compatibleAll (Node _ [l, r]) cs =-  all (bpcompatible (bipartition l)) cs && all (bpcompatible (bipartition r)) cs-compatibleAll _ _ = error "Tree is not bifurcating."--compatibleWith-  :: (Show b, Ord b) => (a -> b) -> [Constraint a] -> Tree a -> Bool-compatibleWith f cs t = compatibleAll (fmap f t) (map (S.map f) cs)--spec :: Spec-spec = do-  describe "subTree" $ do-    it "returns nothing if no leaf satisfies prediacte"-      $          subTree (== 3) smallTree-      `shouldBe` Nothing-    it "returns the correct subtree for a small example"-      $          subTree (== 1) smallTree-      `shouldBe` Just smallSubTree--  describe "pruneWith" $ do-    it "leaves a normal tree untouched"-      $          pruneWith const largeTree-      `shouldBe` largeTree-    it "correctly prunes a small example"-      $          pruneWith const smallSubTree-      `shouldBe` smallSubTreePruned-    it "leaves height constant for Measurable trees"-      $          height (prune subSampleLargeTree)-      `shouldBe` height subSampleLargeTree--  describe "roots" $ do-    it "correctly handles leaves and cherries" $ do-      let tleaf   = Node 0 [] :: Tree Int-          tcherry = Node 0 [Node 1 [], Node 2 []] :: Tree Int-      roots tleaf `shouldBe` [tleaf]-      roots tcherry `shouldBe` [tcherry]-    it "correctly handles simple trees" $ do-      let simpleTre =-            Node "i" [Node "j" [Node "x" [], Node "y" []], Node "z" []] :: Tree-                String-          simpleSol =-            [ Node "i" [Node "j" [Node "x" [], Node "y" []], Node "z" []]-            , Node "i" [Node "x" [], Node "j" [Node "y" [], Node "z" []]]-            , Node "i" [Node "j" [Node "x" [], Node "z" []], Node "y" []]-            ]-      roots simpleTre `shouldBe` simpleSol-    modifyMaxSize (* 100)-      $ it "returns the correct number of rooted trees for arbitrary trees"-      $ property (prop_roots :: (Tree Int -> Bool))--  -- TODO: dropLeafWith, intersect.--  describe "connect" $ modifyMaxSize (* 100) $ do-    it "returns the correct number of rooted trees for arbitrary trees"-      $ property (prop_connect :: Int -> Tree Int -> Tree Int -> Bool)-    it "correctly connects sample trees without and with constraints" $ do-      a <- parseFileWith (oneNewick Standard) "data/ConnectA.tree"-      b <- parseFileWith (oneNewick Standard) "data/ConnectB.tree"-      c <- parseFileWith (manyNewick Standard) "data/ConnectConstraints.tree"-      let ts  = connect (PhyloLabel "" Nothing (Just 1.0)) a b-          cs  = concatMap clades c :: [Constraint (PhyloLabel L.ByteString)]-          ts' = filter (compatibleWith getName cs) ts-      length ts `shouldBe` 63-      length ts' `shouldBe` 15
test/ELynx/Export/Tree/NewickSpec.hs view
@@ -1,61 +1,53 @@-{- |-Module      :  ELynx.Export.Tree.NewickSpec-Copyright   :  (c) Dominik Schrempf 2020-License     :  GPL-3.0-or-later--Maintainer  :  dominik.schrempf@gmail.com-Stability   :  unstable-Portability :  portable--Creation date: Wed Feb 13 11:01:53 2019.---}-+-- |+-- Module      :  ELynx.Export.Tree.NewickSpec+-- Copyright   :  (c) Dominik Schrempf 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Wed Feb 13 11:01:53 2019. module ELynx.Export.Tree.NewickSpec-  ( spec+  ( spec,   ) where -import qualified Data.ByteString.Lazy.Char8    as L-import           Data.Tree-import           Test.Hspec-import           Test.QuickCheck         hiding ( label )-import           Test.QuickCheck.Instances.Containers-                                                ( )--import           ELynx.Data.Tree-import           ELynx.Data.Tree.PhyloTreeArbitraryInstance-                                                ( )-import           ELynx.Export.Tree.Newick-import           ELynx.Import.Tree.Newick-import           ELynx.Tools+import Data.Bifunctor+import Data.ByteString.Lazy.Char8 (ByteString)+import qualified Data.ByteString.Lazy.Char8 as L+import ELynx.Data.Tree+import ELynx.Data.Tree.Arbitrary ()+import ELynx.Export.Tree.Newick+import ELynx.Import.Tree.Newick+import ELynx.Tools+import Test.Hspec+import Test.QuickCheck hiding (label) -samplePhyloByteStringTree :: Tree (PhyloLabel L.ByteString)-samplePhyloByteStringTree = parseByteStringWith "Sample newick byte string"-                                                (newick Standard)-                                                sampleNewickByteString1+samplePhyloByteStringTree :: Tree Phylo ByteString+samplePhyloByteStringTree = parseByteStringWith "Sample newick byte string" (newick Standard) sampleNewickByteString1 -sampleNewickByteString1 :: L.ByteString+sampleNewickByteString1 :: ByteString sampleNewickByteString1 =   L.pack     "(Aeropyrum0:0.5478645225,(((((((((Arabidopsi:0.0701001024,Oryza_sati:0.0765988261):0.0309636193,Gymnosperm:0.0520325624):0.0338982245,Physcomitr:0.0768008916):0.0895714685,(Chlamydomo:0.1136227755,Dunaliella:0.1406347323):0.1117340620):0.0818876186,Rhodophyta:0.3405656487):0.0363527066,((((((Babesia_bo:0.1646969208,Theileria0:0.1519889486):0.1908081096,Plasmodium:0.3250696762):0.0637865908,(Toxoplasma:0.1153570425,Eimeria000:0.1671916078):0.0980136930):0.0518956330,Cryptospor:0.3175062809):0.1607708388,Ciliophora:0.5687502950):0.0624078848,(Phytophtho:0.2016424948,((Thalassios:0.1202730781,Phaeodacty:0.1290341329):0.1772775509,Phaeophyce:0.1989260715):0.0312359673):0.1154768302):0.0311952864):0.0149160316,(((((((((Candida_al:0.1027755272,Saccharomy:0.1190206560):0.1333487870,Neurospora:0.1977309079):0.0522926266,Schizosacc:0.2019603227):0.0567441011,(Cryptococc:0.1948614959,Ustilago_m:0.1564451295):0.0775729694):0.0323959951,Glomus_int:0.1573670796):0.0194701292,Chytridiom:0.2228415254):0.0384370601,Encephalit:1.4622174644):0.0416231688,(((Drosophila:0.2160627753,(Mammalians:0.1080484094,Tunicates0:0.1739253014):0.0289624371):0.0346633757,Hydrozoa00:0.2058137032):0.0480963050,Monosiga_b:0.3020637584):0.0654894239):0.0380915725,(Dictyostel:0.3453588998,Mastigamoe:0.3844779231):0.0478795653):0.0129578395):1.7592083381,((Archaeoglo:0.5402784445,Methanococ:0.4088567459):0.0993669265,Pyrococcus:0.4058713829):0.1734405968):0.2193511807,Pyrobaculu:0.7507718047):0.1646616482,Sulfolobus:0.5404967897);" -prop_invariant :: Tree (PhyloLabel Int) -> Bool+prop_invariant :: Tree Phylo Int -> Bool prop_invariant t =   parseByteStringWith "Sample newick byte string" (newick Standard) (toNewick t)     ==-                   -- We have to convert Int to ByteString, but we cannot use-                   -- arbitrary instances with ByteStrings because many-                   -- characters are disallowed and break the Newick string.-       fmap (\l -> l { label = getName l }) t+    -- We have to convert Int to ByteString, but we cannot use+    -- arbitrary instances with ByteStrings because many+    -- characters are disallowed and break the Newick string.+    second getName t  spec :: Spec spec = describe "parseByteStringWith newick $ toNewickPhyloByteString" $ do   it "should be an invariant for a sample tree" $ do     let bs = toNewick samplePhyloByteStringTree-    parseByteStringWith "Newick string converted from tree object"-                        (newick Standard)-                        bs+    parseByteStringWith+      "Newick string converted from tree object"+      (newick Standard)+      bs       `shouldBe` samplePhyloByteStringTree-   it "should be an invariant for a sample tree" $ property prop_invariant
+ test/ELynx/Export/Tree/NexusSpec.hs view
@@ -0,0 +1,46 @@+{-# LANGUAGE OverloadedStrings #-}++-- |+-- Module      :  ELynx.Export.Tree.NexusSpec+-- Description :  Test export of trees in Nexus files+-- Copyright   :  (c) Dominik Schrempf 2020+-- License     :  GPL-3+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Tue Apr 28 18:08:14 2020.+module ELynx.Export.Tree.NexusSpec+  ( spec,+  )+where++import Data.ByteString.Lazy (ByteString)+import ELynx.Data.Tree+import ELynx.Export.Tree.Nexus+import ELynx.Import.Tree.Newick (NewickFormat (..))+import ELynx.Import.Tree.Nexus+import ELynx.Tools+import Test.Hspec++tree :: Tree Phylo ByteString+tree =+  Node+    (Phylo Nothing Nothing)+    ""+    [ Node+        (Phylo Nothing Nothing)+        ""+        [ Node (Phylo Nothing Nothing) "A" [],+          Node (Phylo Nothing Nothing) "B" []+        ],+      Node (Phylo Nothing Nothing) "C" []+    ]++spec :: Spec+spec = describe "toNexusTrees" $+  it "exports a nexus file with a TREES block" $+    do+      let ts = parseByteStringWith "NexusTrees" (nexusTrees Standard) (toNexusTrees [("tree1", tree)])+      head ts `shouldBe` ("tree1", tree)
test/ELynx/Import/Tree/NewickSpec.hs view
@@ -1,98 +1,83 @@-{- |-Module      :  ELynx.Import.Tree.NewickSpec-Copyright   :  (c) Dominik Schrempf 2020-License     :  GPL-3.0-or-later--Maintainer  :  dominik.schrempf@gmail.com-Stability   :  unstable-Portability :  portable--Creation date: Fri Jan 18 10:14:04 2019.---}+{-# LANGUAGE OverloadedStrings #-} +-- |+-- Module      :  ELynx.Import.Tree.NewickSpec+-- Copyright   :  (c) Dominik Schrempf 2020+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Fri Jan 18 10:14:04 2019. module ELynx.Import.Tree.NewickSpec-  ( spec+  ( spec,   ) where -import qualified Data.ByteString.Lazy.Char8    as L-import           Data.Tree-import           Test.Hspec-import           Test.Hspec.Megaparsec-import           Text.Megaparsec+import Data.ByteString.Lazy.Char8 (ByteString)+import qualified Data.ByteString.Lazy.Char8 as L+import ELynx.Data.Tree+import ELynx.Import.Tree.Newick+import ELynx.Tools+import Test.Hspec+import Test.Hspec.Megaparsec+import Text.Megaparsec -import           ELynx.Tools-import           ELynx.Data.Tree-import           ELynx.Import.Tree.Newick+-- sampleLabelByteString :: ByteString+-- sampleLabelByteString = "name:0.3" -sampleLabelByteString :: L.ByteString-sampleLabelByteString = L.pack "name:0.3"+-- sampleLeaf :: Tree Phylo ByteString+-- sampleLeaf = Node (Phylo (Just 0.3) Nothing) "name" [] -sampleLabel :: (PhyloLabel L.ByteString)-sampleLabel = PhyloLabel (L.pack "name") Nothing (Just 0.3)+-- sampleForestByteString :: ByteString+-- sampleForestByteString = "(l,l,(a,b))" -sampleForestByteString :: L.ByteString-sampleForestByteString = L.pack "(l,l,(a,b))"+noPL :: Phylo+noPL = Phylo Nothing Nothing -sampleForest :: [Tree (PhyloLabel L.ByteString)]-sampleForest =-  [ Node { rootLabel = PhyloLabel (L.pack "l") Nothing Nothing, subForest = [] }-  , Node { rootLabel = PhyloLabel (L.pack "l") Nothing Nothing, subForest = [] }-  , Node-    { rootLabel = PhyloLabel (L.pack "") Nothing Nothing-    , subForest = [ Node { rootLabel = PhyloLabel (L.pack "a") Nothing Nothing-                         , subForest = []-                         }-                  , Node { rootLabel = PhyloLabel (L.pack "b") Nothing Nothing-                         , subForest = []-                         }-                  ]-    }-  ]+-- sampleForest :: Forest Phylo ByteString+-- sampleForest =+--   [ Node noPL "l" [],+--     Node noPL "l" [],+--     Node+--       noPL+--       ""+--       [ Node noPL "a" [],+--         Node noPL "b" []+--       ]+--   ] -sampleNewickByteString1 :: L.ByteString+sampleNewickByteString1 :: ByteString sampleNewickByteString1 =   L.pack     "(Aeropyrum0:0.5478645225,(((((((((Arabidopsi:0.0701001024,Oryza_sati:0.0765988261):0.0309636193,Gymnosperm:0.0520325624):0.0338982245,Physcomitr:0.0768008916):0.0895714685,(Chlamydomo:0.1136227755,Dunaliella:0.1406347323):0.1117340620):0.0818876186,Rhodophyta:0.3405656487):0.0363527066,((((((Babesia_bo:0.1646969208,Theileria0:0.1519889486):0.1908081096,Plasmodium:0.3250696762):0.0637865908,(Toxoplasma:0.1153570425,Eimeria000:0.1671916078):0.0980136930):0.0518956330,Cryptospor:0.3175062809):0.1607708388,Ciliophora:0.5687502950):0.0624078848,(Phytophtho:0.2016424948,((Thalassios:0.1202730781,Phaeodacty:0.1290341329):0.1772775509,Phaeophyce:0.1989260715):0.0312359673):0.1154768302):0.0311952864):0.0149160316,(((((((((Candida_al:0.1027755272,Saccharomy:0.1190206560):0.1333487870,Neurospora:0.1977309079):0.0522926266,Schizosacc:0.2019603227):0.0567441011,(Cryptococc:0.1948614959,Ustilago_m:0.1564451295):0.0775729694):0.0323959951,Glomus_int:0.1573670796):0.0194701292,Chytridiom:0.2228415254):0.0384370601,Encephalit:1.4622174644):0.0416231688,(((Drosophila:0.2160627753,(Mammalians:0.1080484094,Tunicates0:0.1739253014):0.0289624371):0.0346633757,Hydrozoa00:0.2058137032):0.0480963050,Monosiga_b:0.3020637584):0.0654894239):0.0380915725,(Dictyostel:0.3453588998,Mastigamoe:0.3844779231):0.0478795653):0.0129578395):1.7592083381,((Archaeoglo:0.5402784445,Methanococ:0.4088567459):0.0993669265,Pyrococcus:0.4058713829):0.1734405968):0.2193511807,Pyrobaculu:0.7507718047):0.1646616482,Sulfolobus:0.5404967897);" -sampleNewickByteString2 :: L.ByteString+sampleNewickByteString2 :: ByteString sampleNewickByteString2 =   L.pack     "(Caenorhabd:0.0176707431,C0briggsae:0.0142817073,(Ancylostom:0.0711440844,(Pristionch:0.1301309005,((Brugia_mal:0.0757534325,Ascaris0su:0.0482660407)1:0.0563924634,(((Meloidogyn:0.1239621893,Heteroderi:0.0987968800)1:0.1136879428,Strongyloi:0.2483437292)1:0.0252467381,(Trichoceph:0.2985037612,((((((Coleoptera:0.0907850846,(Apis0melli:0.0754058285,Hemiptera0:0.1675359618)0.93:0.0085703192)1:0.0146980945,(Siphonapte:0.0556805916,Bombyx0mor:0.0968983509)1:0.0127867903)1:0.0167360185,((Drosophila:0.0492149086,Glossina0m:0.0534390467)1:0.0583462602,Anopheles0:0.0968919941)1:0.0431343553)1:0.0535616453,Crustacea0:0.2247268999)1:0.0252755187,Chelicerat:0.1537491558)1:0.0212497286,((Echinoderm:0.1803896615,(Cephalocho:0.1492264574,(Urochordat:0.2194747834,(Mammalia00:0.0393008407,Actinopter:0.0491700096):0.0858550024)1:0.0157515969)1:0.0132516777)1:0.0203423736,((((((((Neurospora:0.0721607581,Magnaporth:0.0814182810)1:0.0198940548,Gibberella:0.0858192964)1:0.0533872590,Eurotiomyc:0.1058840539)1:0.1266302603,(Candida0al:0.1349957509,Saccharomy:0.1553464572)1:0.1791344287)1:0.0529664967,Schizosacc:0.2550087905)1:0.0723650615,(Ustilago0m:0.2031812772,(Homobasidi:0.1473391802,Cryptococc:0.2070743149)1:0.0347868586)1:0.0790327507)1:0.0727415175,Glomales00:0.1779430068)1:0.0169066667,Chytridiom:0.3028920870)1:0.3311420273)1:0.0278566156)1:0.1049569161)1:0.1366217350)1:0.0171168289)1:0.0345725378)1:0.0542036935)1:0.0879337167)1;" --sampleNewickEmptyByteString :: L.ByteString+sampleNewickEmptyByteString :: ByteString sampleNewickEmptyByteString = L.pack "(,(,,),);" -sampleNewickEmpty :: Tree (PhyloLabel L.ByteString)-sampleNewickEmpty = Node-  { rootLabel = PhyloLabel (L.pack "") Nothing Nothing-  , subForest =-    [ Node { rootLabel = PhyloLabel (L.pack "") Nothing Nothing-           , subForest = []-           }-    , Node-      { rootLabel = PhyloLabel (L.pack "") Nothing Nothing-      , subForest =-        [ Node { rootLabel = PhyloLabel (L.pack "") Nothing Nothing-               , subForest = []-               }-        , Node { rootLabel = PhyloLabel (L.pack "") Nothing Nothing-               , subForest = []-               }-        , Node { rootLabel = PhyloLabel (L.pack "") Nothing Nothing-               , subForest = []-               }-        ]-      }-    , Node { rootLabel = PhyloLabel (L.pack "") Nothing Nothing-           , subForest = []-           }+sampleNewickEmpty :: Tree Phylo ByteString+sampleNewickEmpty =+  Node+    noPL+    ""+    [ Node noPL "" [],+      Node+        noPL+        ""+        [ Node noPL "" [],+          Node noPL "" [],+          Node noPL "" []+        ],+      Node noPL "" []     ]-  } -sampleTreeNewickRevbayes :: L.ByteString+sampleTreeNewickRevbayes :: ByteString sampleTreeNewickRevbayes =   L.pack "[&R](l[IDL]:0.3[KEYVALPAIRS],r[IDR]:0.4[KEYVALPARIS])[ID]:0.3;" @@ -101,55 +86,54 @@  spec :: Spec spec = do-  describe "branchLength" $ do-    it "parses a colon and a branch length"-      $             parse branchLength "" (L.pack ":13.2")-      `shouldParse` Just 13.2+  -- describe "branchLength" $ do+  --   it "parses a colon and a branch length" $+  --     parse branchLength "" ":13.2"+  --       `shouldParse` Just 13.2 -    it "returns Nothing if no branch length is given"-      $             parse branchLength "" (L.pack "")-      `shouldParse` Nothing+  --   it "returns Nothing if no branch length is given" $+  --     parse branchLength "" ""+  --       `shouldParse` Nothing -  describe "name" $ do-    it "parses a string of printable characters"-      $             parse name "" (L.pack "aName")-      `shouldParse` L.pack "aName"+  -- describe "name" $ do+  --   it "parses a string of printable characters" $+  --     parse name "" "aName"+  --       `shouldParse` "aName" -    it "parses blanks, colons, semicolons, parentheses, and sequare brackets"-      $             parse name "" (L.pack "aName bla")-      `shouldParse` L.pack "aName"+  --   it "parses blanks, colons, semicolons, parentheses, and sequare brackets" $+  --     parse name "" "aName bla"+  --       `shouldParse` "aName" -    it "allows empty names" $ parse name "" (L.pack "") `shouldParse` L.pack ""+  --   it "allows empty names" $ parse name "" "" `shouldParse` L.pack "" -  describe "node" $ do-    it "parses a tree node"-      $             parse node "" sampleLabelByteString-      `shouldParse` sampleLabel+  -- describe "node" $ do+  --   it "parses a tree node" $+  --     parse node "" sampleLabelByteString+  --       `shouldParse` sampleLabel -    it "parses tree nodes with empty names and branch lengths"-      $             parse node "" (L.pack "")-      `shouldParse` PhyloLabel (L.pack "") Nothing Nothing+  --   it "parses tree nodes with empty names and branch lengths" $+  --     parse node "" ""+  --       `shouldParse` PhyloLabelSoft "" Nothing Nothing -  describe "leaf"-    $             it "parses a leaf of a tree"-    $             parse leaf "" sampleLabelByteString-    `shouldParse` Node sampleLabel []+  -- describe "leaf" $+  --   it "parses a leaf of a tree" $+  --     parse leaf "" sampleLabelByteString+  --       `shouldParse` sampleLeaf -  describe "forest"-    $             it "parses a set of trees within brackets"-    $             parse forest "" sampleForestByteString-    `shouldParse` sampleForest+  -- describe "forestP" $+  --   it "parses a set of trees within brackets" $+  --     parse forestP "" sampleForestByteString+  --       `shouldParse` sampleForest    describe "newick" $ do-    it "parses a newick trees" $ do+    it "parses newick trees" $ do       parse (newick Standard) "" `shouldSucceedOn` sampleNewickByteString1       parse (newick Standard) "" `shouldSucceedOn` sampleNewickByteString2--    it "parses a weird newick tree without node labels nor branch lengths"-      $             parse (newick Standard) "" sampleNewickEmptyByteString-      `shouldParse` sampleNewickEmpty--  describe "newickRevBayes" $ it "parses newick trees in RevBayes format" $ do-    parse (newick RevBayes) "" `shouldSucceedOn` sampleTreeNewickRevbayes-    t2 <- parseFileWith (newick RevBayes) sampleNewickRevBayesFile-    length (leaves t2) `shouldBe` 102+    it "parses a weird newick tree without node labels nor branch lengths" $+      parse (newick Standard) "" sampleNewickEmptyByteString `shouldParse` sampleNewickEmpty+  describe "newickRevBayes" $+    it "parses newick trees in RevBayes format" $+      do+        parse (newick RevBayes) "" `shouldSucceedOn` sampleTreeNewickRevbayes+        t2 <- parseFileWith (newick RevBayes) sampleNewickRevBayesFile+        length (leaves t2) `shouldBe` 102
+ test/ELynx/Import/Tree/NexusSpec.hs view
@@ -0,0 +1,51 @@+{-# LANGUAGE OverloadedStrings #-}++-- |+-- Module      :  ELynx.Import.Tree.NexusSpec+-- Description :  Test import of trees in Nexus files+-- Copyright   :  (c) Dominik Schrempf 2020+-- License     :  GPL-3+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  unstable+-- Portability :  portable+--+-- Creation date: Tue Apr 28 18:08:14 2020.+module ELynx.Import.Tree.NexusSpec+  ( spec,+  )+where++import Data.ByteString.Lazy (ByteString)+import ELynx.Data.Tree+import ELynx.Import.Tree.Newick (NewickFormat (..))+import ELynx.Import.Tree.Nexus+import ELynx.Tools+import Test.Hspec++file :: FilePath+file = "data/SimpleTree.nex"++noPL :: Phylo+noPL = Phylo Nothing Nothing++res :: Tree Phylo ByteString+res =+  Node+    noPL+    ""+    [ Node+        noPL+        ""+        [ Node noPL "A" [],+          Node noPL "B" []+        ],+      Node noPL "C" []+    ]++spec :: Spec+spec = describe "trees" $+  it "parses a nexus file with a TREES block" $+    do+      ts <- parseFileWith (nexusTrees Standard) file+      head ts `shouldBe` ("tree1", res)