elynx-tree 0.1.0 → 0.2.1
raw patch · 36 files changed
+2076/−1307 lines, 36 filesdep −lifted-asyncdep −optparse-applicativedep −paralleldep ~basedep ~bytestringdep ~containers
Dependencies removed: lifted-async, optparse-applicative, parallel, transformers, vector
Dependency ranges changed: base, bytestring, containers, criterion, elynx-tools, hspec-megaparsec, math-functions, megaparsec, primitive, statistics
Files
- ChangeLog.md +8/−5
- README.md +375/−88
- bench/Bench.hs +6/−10
- elynx-tree.cabal +34/−53
- src/ELynx/Data/Tree.hs +63/−0
- src/ELynx/Data/Tree/Bipartition.hs +86/−67
- src/ELynx/Data/Tree/BranchSupportTree.hs +23/−22
- src/ELynx/Data/Tree/Distance.hs +102/−61
- src/ELynx/Data/Tree/EvoTree.hs +0/−44
- src/ELynx/Data/Tree/MeasurableTree.hs +48/−28
- src/ELynx/Data/Tree/Multipartition.hs +67/−63
- src/ELynx/Data/Tree/NamedTree.hs +7/−6
- src/ELynx/Data/Tree/PhyloTree.hs +10/−20
- src/ELynx/Data/Tree/SubSample.hs +64/−0
- src/ELynx/Data/Tree/Subset.hs +0/−105
- src/ELynx/Data/Tree/SumStat.hs +15/−12
- src/ELynx/Data/Tree/Tree.hs +135/−112
- src/ELynx/Distribution/BirthDeath.hs +34/−26
- src/ELynx/Distribution/BirthDeathCritical.hs +23/−18
- src/ELynx/Distribution/BirthDeathCriticalNoTime.hs +17/−15
- src/ELynx/Distribution/BirthDeathNearlyCritical.hs +44/−27
- src/ELynx/Distribution/CoalescentContinuous.hs +7/−5
- src/ELynx/Distribution/TimeOfOrigin.hs +33/−24
- src/ELynx/Distribution/TimeOfOriginNearCritical.hs +35/−26
- src/ELynx/Distribution/Types.hs +3/−2
- src/ELynx/Export/Tree/Newick.hs +29/−32
- src/ELynx/Import/Tree/Newick.hs +159/−55
- src/ELynx/Simulate/Coalescent.hs +32/−34
- src/ELynx/Simulate/PointProcess.hs +121/−85
- test/ELynx/Data/Tree/BipartitionSpec.hs +54/−49
- test/ELynx/Data/Tree/DistanceSpec.hs +182/−44
- test/ELynx/Data/Tree/MultipartitionSpec.hs +30/−26
- test/ELynx/Data/Tree/PhyloTreeArbitraryInstance.hs +33/−0
- test/ELynx/Data/Tree/TreeSpec.hs +73/−68
- test/ELynx/Export/Tree/NewickSpec.hs +27/−18
- test/ELynx/Import/Tree/NewickSpec.hs +97/−57
ChangeLog.md view
@@ -1,10 +1,13 @@ -# Table of Contents--1. [Changelog for ELynx](#org263dcac)+# Changelog for ELynx -<a id="org263dcac"></a>+## Version 0.2.1 -# Changelog for ELynx+- Validation and repetition of previous analyses is finally possible with the+ new `elynx` binary.+- A library `elynx-markov` for running Markov processes along phylogenetic trees+ has been split off `elynx-seq`. This library performs the computations when+ executing `slynx simulate ...`.+- Many other small improvements.
README.md view
@@ -2,30 +2,31 @@ # The ELynx Suite -Version: 0.1.0.+Version: 0.2.1. Reproducible evolution made easy. -The ELynx Suite is a Haskell library and a tool set for computational biology.-The goal of the ELynx Suite is reproducible research. Evolutionary sequences and-phylogenetic trees can be read, viewed, modified and simulated. Exact-specification of all options is necessary, and nothing is assumed about the data-(e.g., the type of the genetic code). The command line with all arguments is-consistently, and automatically logged. The work overhead in the beginning-usually pays off in the end.+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. -The Elynx Suite consists of three library packages and two executables providing+The Elynx Suite consists of four library packages and three executables providing a range of sub commands. The library packages are: +- **elynx-markov:** Simulate multi sequence alignments along phylogenetic trees. - **elynx-seq:** Handle evolutionary sequences and multi sequence alignments.-- **elynx-tree:** Handle phylogenetic trees. - **elynx-tools:** Tools for the provided executables.+- **elynx-tree:** Handle phylogenetic trees. The executables are: -- **slynx:** Analyze, modify, and simulate evolutionary sequences (FASTA format).-- **tlynx:** Analyze, modify, and simulate phylogenetic trees (Newick format).+- **slynx:** Analyze, modify, and simulate evolutionary sequences.+- **tlynx:** Analyze, modify, and simulate phylogenetic trees.+- **elynx:** Validate and redo past analyses. **ELynx is actively developed. We happily receive comments, ideas, feature requests, and pull requests!**@@ -67,10 +68,12 @@ slynx --help - ELynx Suite version 0.1.0. Developed by Dominik Schrempf. Compiled on January- 30, 2020, at 09:18 am, UTC.+ ELynx Suite version 0.2.1.+ Developed by Dominik Schrempf.+ Compiled on April 27, 2020, at 06:31 am, UTC. - Usage: slynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME] COMMAND+ Usage: slynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME] + [-f|--force] COMMAND Analyze, and simulate multi sequence alignments. Available options:@@ -80,21 +83,22 @@ Debug (default: Info) -o,--output-file-basename NAME Specify base name of output file+ -f,--force Ignore previous analysis and overwrite existing+ output files. Available commands: concatenate Concatenate sequences found in input files.- examine Examine sequences. If data is a multi sequence- alignment, additionally analyze columns.+ examine Examine sequences. If data is a multi sequence alignment, additionally analyze columns.+ filter-columns Filter columns of multi sequence alignments. filter-rows Filter rows (or sequences) found in input files.- filter-columns Filter columns of multi-sequence alignments. simulate Simulate multi sequence alignments. sub-sample Sub-sample columns from multi sequence alignments. translate Translate from DNA to Protein or DNAX to ProteinX. - File formats:+ Available sequence file formats: - FASTA - Alphabet types:+ Available alphabets: - DNA (nucleotides) - DNAX (nucleotides; including gaps) - DNAI (nucleotides; including gaps, and IUPAC codes)@@ -103,19 +107,20 @@ - ProteinS (amino acids; including gaps, and translation stops) - ProteinI (amino acids; including gaps, translation stops, and IUPAC codes) - The ELynx Suite- ---------------- A Haskell library and a tool set for computational biology. The goal of the- ELynx Suite is reproducible research. Evolutionary sequences and phylogenetic- trees can be read, viewed, modified and simulated. Exact specification of all- options is necessary, and nothing is assumed about the data (e.g., the type of- code). The command line with all arguments is consistently, and automatically- logged.+ 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 specific commands:+ Get help for sub commands: slynx examine --help @@ -125,13 +130,36 @@ slynx concatenate --help + ELynx Suite version 0.2.1.+ Developed by Dominik Schrempf.+ Compiled on April 27, 2020, at 06:31 am, UTC.+ Usage: slynx concatenate (-a|--alphabet NAME) INPUT-FILE Concatenate sequences found in input files. Available options:+ -h,--help Show this help text+ -V,--version Show version -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@@ -140,15 +168,37 @@ slynx examine --help - Usage: slynx examine (-a|--alphabet NAME) [INPUT-FILE] [--per-site]- Examine sequences. If data is a multi sequence alignment, additionally analyze- columns.+ ELynx Suite version 0.2.1.+ Developed by Dominik Schrempf.+ Compiled on April 27, 2020, at 06:31 am, UTC. + Usage: slynx examine (-a|--alphabet NAME) INPUT-FILE [--per-site]+ Examine sequences. If data is a multi sequence alignment, additionally analyze columns.+ Available options:+ -h,--help Show this help text+ -V,--version Show version -a,--alphabet NAME Specify alphabet type NAME 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@@ -157,12 +207,17 @@ slynx filter-rows --help - Usage: slynx filter-rows (-a|--alphabet NAME) [INPUT-FILE]- [--longer-than LENGTH] [--shorter-than LENGTH]- [--standard-characters]+ ELynx Suite version 0.2.1.+ Developed by Dominik Schrempf.+ Compiled on April 27, 2020, at 06:31 am, UTC.+ + Usage: slynx filter-rows (-a|--alphabet NAME) INPUT-FILE [--longer-than LENGTH] + [--shorter-than LENGTH] [--standard-characters] Filter rows (or sequences) found in input files. Available options:+ -h,--help Show this help text+ -V,--version Show version -a,--alphabet NAME Specify alphabet type NAME INPUT-FILE Read sequences from INPUT-FILE --longer-than LENGTH Only keep sequences longer than LENGTH@@ -170,21 +225,61 @@ --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 - Usage: slynx filter-columns (-a|--alphabet NAME) [INPUT-FILE]+ ELynx Suite version 0.2.1.+ Developed by Dominik Schrempf.+ Compiled on April 27, 2020, at 06:31 am, UTC.+ + Usage: slynx filter-columns (-a|--alphabet NAME) INPUT-FILE [--standard-chars DOUBLE]- Filter columns of multi-sequence alignments.+ Filter columns of multi sequence alignments. Available options:+ -h,--help Show this help text+ -V,--version Show version -a,--alphabet NAME Specify alphabet type NAME INPUT-FILE Read sequences from INPUT-FILE --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@@ -193,16 +288,22 @@ slynx simulate --help - Usage: slynx simulate (-t|--tree-file Name) [-s|--substitution-model MODEL]- [-m|--mixture-model MODEL] [-e|--edm-file NAME]- [-p|--siteprofile-files NAMES]- [-w|--mixture-model-weights "[DOUBLE,DOUBLE,...]"]+ ELynx Suite version 0.2.1.+ Developed by Dominik Schrempf.+ Compiled on April 27, 2020, at 06:31 am, UTC.+ + Usage: slynx simulate (-t|--tree-file Name) [-s|--substitution-model MODEL] + [-m|--mixture-model MODEL] [-e|--edm-file NAME] + [-p|--siteprofile-files NAMES] + [-w|--mixture-model-weights "[DOUBLE,DOUBLE,...]"] [-g|--gamma-rate-heterogeneity "(NCAT,SHAPE)"] (-l|--length NUMBER) [-S|--seed [INT]] Simulate multi sequence alignments. Available options:- -t,--tree-file Name Read trees from file NAME+ -h,--help Show this help text+ -V,--version Show version+ -t,--tree-file Name Read tree from Newick file NAME -s,--substitution-model MODEL Set the phylogenetic substitution model; available models are shown below (mutually exclusive with -m@@ -224,17 +325,24 @@ Substitution models: -s "MODEL[PARAMETER,PARAMETER,...]{STATIONARY_DISTRIBUTION}"- Supported DNA models: JC, HKY.+ Supported DNA models: JC, F81, HKY, GTR4. For example, -s HKY[KAPPA]{DOUBLE,DOUBLE,DOUBLE,DOUBLE}- Supported Protein models: Poisson, Poisson-Custom, LG, LG-Custom, WAG, WAG-Custom.+ -s GTR4[e_AC,e_AG,e_AT,e_CG,e_CT,e_GT]{DOUBLE,DOUBLE,DOUBLE,DOUBLE}+ where the 'e_XY' are the exchangeabilities from nucleotide X to Y.+ Supported Protein models: Poisson, Poisson-Custom, LG, LG-Custom, WAG, WAG-Custom, GTR20. MODEL-Custom means that only the exchangeabilities of MODEL are used, and a custom stationary distribution is provided. For example,+ -s LG -s LG-Custom{...}+ -s GTR20[e_AR,e_AN,...]{...}+ the 'e_XY' are the exchangeabilities from amino acid X to Y (alphabetical order).+ Notes: The F81 model for DNA is equivalent to the Poisson-Custom for proteins.+ The GTR4 model for DNA is equivalent to the GTR20 for proteins. Mixture models:- -m "MIXTURE(SUBSTITUTION_MODEL_1,SUBSTITUTION_MODEL_2)"+ -m "MIXTURE(SUBSTITUTION_MODEL_1,SUBSTITUTION_MODEL_2[PARAMETERS]{STATIONARY_DISTRIBUTION},...)" For example, -m "MIXTURE(JC,HKY[6.0]{0.3,0.2,0.2,0.3})" Mixture weights have to be provided with the -w option.@@ -249,6 +357,22 @@ LG exchangeabilities with stationary distributions given in FILE. -m "EDM(LG-Custom)" -e FILE 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@@ -257,12 +381,18 @@ slynx sub-sample --help - Usage: slynx sub-sample (-a|--alphabet NAME) [INPUT-FILE]+ ELynx Suite version 0.2.1.+ Developed by Dominik Schrempf.+ Compiled on April 27, 2020, at 06:31 am, UTC.+ + Usage: slynx sub-sample (-a|--alphabet NAME) INPUT-FILE (-n|--number-of-sites INT) (-m|--number-of-alignments INT) [-S|--seed [INT]] Sub-sample columns from multi sequence alignments. Available options:+ -h,--help Show this help text+ -V,--version Show version -a,--alphabet NAME Specify alphabet type NAME INPUT-FILE Read sequences from INPUT-FILE -n,--number-of-sites INT Number of sites randomly drawn with replacement@@ -272,8 +402,23 @@ integers with up to 256 elements (default: random) -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.+ 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@@ -282,17 +427,40 @@ slynx translate --help - Usage: slynx translate (-a|--alphabet NAME) [INPUT-FILE]- (-r|--reading-frame INT) (-u|--universal-code CODE)+ ELynx Suite version 0.2.1.+ Developed by Dominik Schrempf.+ Compiled on April 27, 2020, at 06:31 am, UTC.+ + Usage: slynx translate (-a|--alphabet NAME) INPUT-FILE (-r|--reading-frame INT)+ (-u|--universal-code CODE) Translate from DNA to Protein or DNAX to ProteinX. Available options:+ -h,--help Show this help text+ -V,--version Show version -a,--alphabet NAME Specify alphabet type NAME INPUT-FILE Read sequences from INPUT-FILE -r,--reading-frame INT Reading frame [0|1|2]. -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@@ -301,10 +469,12 @@ tlynx --help - ELynx Suite version 0.1.0. Developed by Dominik Schrempf. Compiled on January- 30, 2020, at 09:18 am, UTC.+ ELynx Suite version 0.2.1.+ Developed by Dominik Schrempf.+ Compiled on April 27, 2020, at 06:31 am, UTC. - Usage: tlynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME] COMMAND+ Usage: tlynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME] + [-f|--force] COMMAND Compare, examine, and simulate phylogenetic trees. Available options:@@ -314,35 +484,37 @@ Debug (default: Info) -o,--output-file-basename NAME Specify base name of output file+ -f,--force Ignore previous analysis and overwrite existing+ output files. Available commands:+ coalesce Simulate phylogenetic trees using the coalescent processes (see also the 'simulate' command for simulations using the birth and death process).+ compare Compare two phylogenetic trees (compute distances and branch-wise differences).+ connect Connect two phylogenetic trees in all ways (possibly honoring constraints). distance Compute distances between many phylogenetic trees. examine Compute summary statistics of phylogenetic trees.- simulate Simulate phylogenetic trees using birth and death- processes.- compare Compare two phylogenetic trees (compute distances and- branch-wise differences).- connect Connect two phylogenetic trees in all ways (possibly- honoring constraints).- shuffle Shuffle a phylogenetic tree (keep coalescent times,- but shuffle topology and leaves).+ 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). - File formats:- - Newick+ 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. - The ELynx Suite- ---------------- A Haskell library and a tool set for computational biology. The goal of the- ELynx Suite is reproducible research. Evolutionary sequences and phylogenetic- trees can be read, viewed, modified and simulated. Exact specification of all- options is necessary, and nothing is assumed about the data (e.g., the type of- code). The command line with all arguments is consistently, and automatically- logged.+ 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 specific commands:+ Get help for sub commands: slynx examine --help @@ -352,18 +524,43 @@ tlynx compare --help - Usage: tlynx compare [-n|--normalize] [-b|--bipartitions] [-i|--newick-iqtree]- NAME- Compare two phylogenetic trees (compute distances and branch-wise- differences).+ ELynx Suite version 0.2.1.+ Developed by Dominik Schrempf.+ Compiled on April 27, 2020, at 06:31 am, UTC. + Usage: tlynx compare [-n|--normalize] [-b|--bipartitions] [-t|--intersect] + [-f|--newick-format FORMAT] NAME+ Compare two phylogenetic trees (compute distances and branch-wise differences).+ Available options:+ -h,--help Show this help text+ -V,--version Show version -n,--normalize Normalize trees before comparison- -b,--bipartitions Print common and missing bipartitions- -i,--newick-iqtree Use IQ-TREE Newick format (internal node labels are- branch support values)+ -b,--bipartitions Print and plot common and missing bipartitions+ -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 -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@@ -372,14 +569,38 @@ tlynx examine --help - Usage: tlynx examine [INPUT-FILE] [-i|--newick-iqtree]+ ELynx Suite version 0.2.1.+ Developed by Dominik Schrempf.+ Compiled on April 27, 2020, at 06:31 am, UTC.+ + Usage: tlynx examine INPUT-FILE [-f|--newick-format FORMAT] Compute summary statistics of phylogenetic trees. Available options:+ -h,--help Show this help text+ -V,--version Show version INPUT-FILE Read trees from INPUT-FILE- -i,--newick-iqtree Use IQ-TREE Newick format (internal node labels are- branch support values)+ -f,--newick-format FORMAT+ Newick tree format; see 'tlynx+ --help' (default: Standard) -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@@ -388,13 +609,19 @@ tlynx simulate --help - Usage: tlynx simulate [-t|--nTrees INT] [-n|--nLeaves INT] [-H|--height DOUBLE]- [-M|--condition-on-mrca] [-l|--lambda DOUBLE]- [-m|--mu DOUBLE] [-r|--rho DOUBLE] [-u|--sub-sample]+ ELynx Suite version 0.2.1.+ Developed by Dominik Schrempf.+ Compiled on April 27, 2020, at 06:31 am, UTC.+ + Usage: tlynx simulate [-t|--nTrees INT] [-n|--nLeaves INT] [-H|--height DOUBLE] + [-M|--condition-on-mrca] [-l|--lambda DOUBLE] + [-m|--mu DOUBLE] [-r|--rho DOUBLE] [-u|--sub-sample] [-s|--summary-statistics] [-S|--seed [INT]]- Simulate phylogenetic trees using birth and death processes.+ Simulate phylogenetic trees using birth and death processes (see also the 'coalesce' command for simulations using the coalescent process). Available options:+ -h,--help Show this help text+ -V,--version Show version -t,--nTrees INT Number of trees (default: 10) -n,--nLeaves INT Number of leaves per tree (default: 5) -H,--height DOUBLE Fix tree height (no default)@@ -413,15 +640,75 @@ 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 -Documentation of the library can be found on [Hackage](https://hackage.haskell.org/):+Validate and (optionally) redo past ELynx analyses. + elynx --help++ ELynx Suite version 0.2.1.+ Developed by Dominik Schrempf.+ Compiled on April 27, 2020, at 06:31 am, UTC.+ + Usage: elynx COMMAND+ Validate and redo past ELynx analyses+ + Available options:+ -h,--help Show this help text+ -V,--version Show version+ + 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-tree](https://hackage.haskell.org/package/elynx-tree) - [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,8 +1,8 @@ {- | Module : Bench Description : Various benchmarks-Copyright : (c) Dominik Schrempf 2019-License : GPL-3+Copyright : (c) Dominik Schrempf 2020+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -13,13 +13,12 @@ -} import Criterion.Main-import qualified Data.ByteString.Lazy.Char8 as L+import qualified Data.ByteString.Lazy.Char8 as L import Data.Tree -import ELynx.Data.Tree.Bipartition-import ELynx.Data.Tree.PhyloTree+import ELynx.Data.Tree import ELynx.Import.Tree.Newick-import ELynx.Tools.InputOutput+import ELynx.Tools treeFileMany :: FilePath treeFileMany = "data/Many.trees"@@ -31,7 +30,4 @@ main = do ts <- getManyTrees defaultMain- [ bgroup "bipartition"- [ bench "manyTrees" $ whnf (map bipartitions) ts- ]- ]+ [bgroup "bipartition" [bench "manyTrees" $ whnf (map bipartitions) ts]]
elynx-tree.cabal view
@@ -1,9 +1,9 @@ cabal-version: 1.12 name: elynx-tree-version: 0.1.0+version: 0.2.1 license: GPL-3 license-file: LICENSE-copyright: Dominik Schrempf (2019)+copyright: Dominik Schrempf (2020) maintainer: dominik.schrempf@gmail.com author: Dominik Schrempf homepage: https://github.com/dschrempf/elynx#readme@@ -23,15 +23,15 @@ library exposed-modules:+ ELynx.Data.Tree ELynx.Data.Tree.Bipartition ELynx.Data.Tree.BranchSupportTree ELynx.Data.Tree.Distance- ELynx.Data.Tree.EvoTree ELynx.Data.Tree.MeasurableTree ELynx.Data.Tree.Multipartition ELynx.Data.Tree.NamedTree ELynx.Data.Tree.PhyloTree- ELynx.Data.Tree.Subset+ ELynx.Data.Tree.SubSample ELynx.Data.Tree.SumStat ELynx.Data.Tree.Tree ELynx.Distribution.BirthDeath@@ -50,24 +50,17 @@ other-modules: Paths_elynx_tree default-language: Haskell2010- ghc-options: -Wall -O2+ ghc-options: -Wall -fllvm build-depends:- QuickCheck >=2.13.2 && <2.14,- base >=4.12.0.0 && <4.13,- bytestring >=0.10.8.2 && <0.11,- containers >=0.6.0.1 && <0.7,- elynx-tools >=0.1.0 && <0.2,- lifted-async >=0.10.0.4 && <0.11,- math-functions >=0.3.2.1 && <0.4,- megaparsec >=7.0.5 && <7.1,+ 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.1 && <0.3,+ math-functions >=0.3.3.0 && <0.4,+ megaparsec >=8.0.0 && <8.1, mwc-random >=0.14.0.0 && <0.15,- optparse-applicative >=0.14.3.0 && <0.15,- parallel >=3.2.2.0 && <3.3,- primitive >=0.6.4.0 && <0.7,- quickcheck-instances >=0.3.22 && <0.4,- statistics >=0.15.0.0 && <0.16,- transformers >=0.5.6.2 && <0.6,- vector >=0.12.0.3 && <0.13+ primitive >=0.7.0.1 && <0.8,+ statistics >=0.15.2.0 && <0.16 test-suite tree-test type: exitcode-stdio-1.0@@ -77,33 +70,28 @@ 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 -O2 -eventlog -threaded -rtsopts- -with-rtsopts=-N+ ghc-options: -Wall -fllvm -threaded -rtsopts -with-rtsopts=-N build-depends: QuickCheck >=2.13.2 && <2.14,- base >=4.12.0.0 && <4.13,- bytestring >=0.10.8.2 && <0.11,- containers >=0.6.0.1 && <0.7,- elynx-tools >=0.1.0 && <0.2,+ 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.1 && <0.3, elynx-tree -any, hspec >=2.7.1 && <2.8,- hspec-megaparsec >=2.0.1 && <2.1,- lifted-async >=0.10.0.4 && <0.11,- math-functions >=0.3.2.1 && <0.4,- megaparsec >=7.0.5 && <7.1,+ 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,- optparse-applicative >=0.14.3.0 && <0.15,- parallel >=3.2.2.0 && <3.3,- primitive >=0.6.4.0 && <0.7,+ primitive >=0.7.0.1 && <0.8, quickcheck-instances >=0.3.22 && <0.4,- statistics >=0.15.0.0 && <0.16,- transformers >=0.5.6.2 && <0.6,- vector >=0.12.0.3 && <0.13+ statistics >=0.15.2.0 && <0.16 benchmark tree-bench type: exitcode-stdio-1.0@@ -112,23 +100,16 @@ other-modules: Paths_elynx_tree default-language: Haskell2010- ghc-options: -Wall -O2 -threaded -rtsopts -with-rtsopts=-N+ ghc-options: -Wall -fllvm -threaded -rtsopts -with-rtsopts=-N build-depends:- QuickCheck >=2.13.2 && <2.14,- base >=4.12.0.0 && <4.13,- bytestring >=0.10.8.2 && <0.11,- containers >=0.6.0.1 && <0.7,- criterion >=1.5.6.0 && <1.6,- elynx-tools >=0.1.0 && <0.2,+ 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.1 && <0.3, elynx-tree -any,- lifted-async >=0.10.0.4 && <0.11,- math-functions >=0.3.2.1 && <0.4,- megaparsec >=7.0.5 && <7.1,+ math-functions >=0.3.3.0 && <0.4,+ megaparsec >=8.0.0 && <8.1, mwc-random >=0.14.0.0 && <0.15,- optparse-applicative >=0.14.3.0 && <0.15,- parallel >=3.2.2.0 && <3.3,- primitive >=0.6.4.0 && <0.7,- quickcheck-instances >=0.3.22 && <0.4,- statistics >=0.15.0.0 && <0.16,- transformers >=0.5.6.2 && <0.6,- vector >=0.12.0.3 && <0.13+ primitive >=0.7.0.1 && <0.8,+ statistics >=0.15.2.0 && <0.16
+ src/ELynx/Data/Tree.hs view
@@ -0,0 +1,63 @@+{- |+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+ ( 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+ )+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++-- -- | An evolutionary label has some information about where the corresponding+-- -- node is on the tree, and if the node is 'extant', 'extinct', 'internal', or+-- -- 'external'. The latter two could also be determined from the tree. This could+-- -- be species, genes or individuals; probably more.+-- class EvoLabel n where+-- extant :: n -> Bool+-- extinct :: n -> Bool++-- internal :: n -> Bool+-- internal n = not $ extant n || extinct n++-- external :: n -> Bool+-- external = not . internal++-- -- -- | Glue branches together, so that one new tree emerges. It's root node is+-- -- -- new, the sub-forest has to be given (a list of trees).+-- -- glue :: (NodeType c)+-- -- => PhyloLabel a b c -- ^ New root node.+-- -- -> [PhyloTree a b c] -- ^ Sub-forest.+-- -- -> PhyloTree a b c+-- -- glue s@(PhyloLabel _ _ n) ts+-- -- | extant n = error "Root node cannot be of type 'Exant'."+-- -- | extinct n = error "Root node cannot be of type 'Extinct'."+-- -- | otherwise = Node s ts
src/ELynx/Data/Tree/Bipartition.hs view
@@ -1,8 +1,8 @@ {- | Module : ELynx.Data.Tree.Bipartition Description : Bipartitions on trees-Copyright : (c) Dominik Schrempf 2019-License : GPL-3+Copyright : (c) Dominik Schrempf 2020+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -26,30 +26,32 @@ module ELynx.Data.Tree.Bipartition ( -- * The 'Bipartition' data type.- Bipartition ()+ Bipartition() , bps , bp , bpmap , bphuman+ , sshow -- * Working with 'Bipartition's. , bipartition , bipartitions , bipartitionToBranchLength- , compatible- ) where+ , bpcompatible+ )+where -import qualified Data.Map as M+import Data.List ( intercalate )+import qualified Data.Map as M import Data.Maybe-import qualified Data.Set as S+import qualified Data.Set as S import Data.Tree -import ELynx.Data.Tree.Subset import ELynx.Data.Tree.Tree --- | Each branch of a tree partitions the leaves of the tree into two--- 'Subset's, 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 :: (Subset a, Subset a) -- ^ Tuple of partitions+-- | 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) @@ -64,20 +66,22 @@ -- | 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 ++ ")"+bphuman f (Bipartition (x, y)) = "(" ++ sshow f x ++ "|" ++ sshow f y ++ ")" --- | Create a bipartition from two 'S.Set's.-bp :: Ord a => Subset a -> Subset a -> Bipartition a-bp xs ys = if xs >= ys- then Bipartition (xs, ys)- else Bipartition (ys, xs)+-- | 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 +-- | 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)+ -- | 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 (smap f x) (smap f y)+bpmap f (Bipartition (x, y)) = bp (S.map f x) (S.map f y) -- | Create a bipartition from two 'S.Set's.-bpwith :: (Ord a, Ord b) => (a -> b) -> Subset a -> Subset a -> Bipartition b+bpwith :: (Ord a, Ord b) => (a -> b) -> S.Set a -> S.Set a -> Bipartition b bpwith f x y = bpmap f $ bp x y instance (Eq a) => Eq (Bipartition a) where@@ -88,83 +92,97 @@ -- Check if a bipartition is valid. For now, only checks if one set is empty. valid :: Bipartition a -> Bool-valid (Bipartition (xs, ys)) = not $ snull xs || snull ys+valid (Bipartition (xs, ys)) = not $ S.null xs || S.null ys -- | For a bifurcating root, get the bipartition induced by the root node. bipartition :: Ord a => Tree a -> Bipartition a-bipartition (Node _ [x, y]) = bp (sfromlist $ leaves x) (sfromlist $ leaves y)-bipartition _ = error "Root node is not bifurcating."+bipartition (Node _ [x, y]) =+ bp (S.fromList $ leaves x) (S.fromList $ leaves y)+bipartition _ = error "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 sempty pTree- else error "bipartitions: The tree contains duplicate leaves."- where ls = leaves t- pTree = partitionTree t+ then S.filter valid $ bipartitionsUnsafe S.empty pTree+ else error "bipartitions: The tree contains duplicate leaves."+ where+ ls = leaves t+ pTree = partitionTree t --- | See 'bipartitions', but do not check if leaves are unique.-bipartitionsUnsafe :: Ord a => Subset a -> Tree (Subset a) -> S.Set (Bipartition a)-bipartitionsUnsafe p (Node l [] ) = S.singleton $ bp p l+-- 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+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+ S.unions+ $ S.singleton (bp p ls)+ : [ bipartitionsUnsafe lvs x | (lvs, x) <- zip lsOthers xs ]+ where lsOthers = subForestGetSubsets p t -- | For a given rose 'Tree', remove all degree two nodes and reconnect the--- resulting disconnected pairs of branches and sum their branch lengths. For--- this operation, a combining binary function and a unity element is required,--- and so we need 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+-- 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+ -> (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, sempty) 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+ 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 '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, Subset a)+type Info c a = (c, S.Set a) -- | 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 ----- > (a, Subset a)-bipartitionToBranchLengthUnsafe :: (Ord a, Ord b, Monoid c)+-- > (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') [] ) =+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, 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)+ 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) -- | Determine compatibility between an bipartition and a subset. If both -- subsets of the bipartition share elements with the given subset, the@@ -172,8 +190,9 @@ -- 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'.-compatible :: (Show a, Ord a) => Bipartition a -> Subset a -> Bool+bpcompatible :: (Show a, Ord a) => Bipartition a -> S.Set a -> Bool -- compatible (Bipartition (l, r)) ss = sintersection l ss `sdisjoint` sintersection r ss-compatible (Bipartition (l, r)) ss = snull lOverlap || snull rOverlap- where lOverlap = sintersection l ss- rOverlap = 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
src/ELynx/Data/Tree/BranchSupportTree.hs view
@@ -1,8 +1,8 @@ {- | Module : ELynx.Data.Tree.BranchSupportTree Description : Node label with branch support-Copyright : (c) Dominik Schrempf 2019-License : GPL-3+Copyright : (c) Dominik Schrempf 2020+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -14,24 +14,26 @@ module ELynx.Data.Tree.BranchSupportTree ( BranchSupport- , BranchSupported (..)- , normalize+ , BranchSupported(..)+ , normalizeBranchSupport , collapse- ) where+ )+where import Data.List import Data.Tree --- TODO: Implementation of a proper data type. When parsing Newick files, we+-- 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? ----- data BranchSupport = BSNothing--- | BSInt Int--- | BSDouble Double deriving (Num)+-- 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'.@@ -45,32 +47,31 @@ setBranchSupport :: BranchSupport -> a -> a apply :: BranchSupported a => (Double -> Double) -> a -> a-apply f l = setBranchSupport (f <$> s) l- where s = getBranchSupport l+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.-normalize :: BranchSupported a => Tree a -> Tree a-normalize t = case mm of+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+accept _ Nothing = True+accept thresh (Just s) = s >= thresh -- | Collapse branches with support lower than given value. Note, branch length--- is ignored at the moment. Continue collapsing until a fix point is reached.+-- 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'+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+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,8 +1,8 @@ {- | Module : ELynx.Data.Tree.Distance Description : Compute distances between trees-Copyright : (c) Dominik Schrempf 2019-License : GPL-3+Copyright : (c) Dominik Schrempf 2020+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -11,8 +11,15 @@ Creation date: Thu Jun 13 17:15:54 2019. Various distance functions for phylogenetic trees (and binary trees in general).-All trees are assumed to be UNROOTED. +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@@ -24,22 +31,29 @@ , branchScoreWith , pairwise , adjacent- ) where+ )+where import Data.List-import qualified Data.Map as M+import qualified Data.Map as M import Data.Monoid-import qualified Data.Set as S+import qualified Data.Set as S import Data.Tree -import ELynx.Data.Tree.Bipartition (Bipartition,- bipartitionToBranchLength,- bipartitions)-import ELynx.Data.Tree.MeasurableTree (Measurable, getLen)-import ELynx.Data.Tree.Multipartition (Multipartition, compatible,- fromBipartition,- multipartitions)+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 ) -- Symmetric difference between two 'Set's. symmetricDifference :: Ord a => S.Set a -> S.Set a -> S.Set a@@ -48,25 +62,35 @@ -- | 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!+-- only. The tree is assumed to be UNROOTED! See+-- http://evolution.genetics.washington.edu/phylip/doc/treedist.html. -- -- XXX: Comparing a list of trees with this function recomputes bipartitions.-symmetricWith :: (Ord b) => (a -> b) -> Tree a -> Tree a -> Int-symmetricWith f t1 t2 = length $ symmetricDifference (bs t1) (bs t2)- where bs t = bipartitions $ fmap f t+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 'id' for comparisons.-symmetric :: (Ord a, Named a) => Tree a -> Tree a -> Int+-- | See 'symmetricWith', but with 'getName' for comparisons.+symmetric :: (Show a, Ord a, Named a) => Tree a -> Tree a -> Int symmetric = symmetricWith getName -countIncompatibilities :: (Ord a, Show a) => S.Set (Bipartition a) -> S.Set (Multipartition a) -> Int-countIncompatibilities bs ms = foldl' (\i b -> if any (compatible (fromBipartition b)) ms- then i- else i+1) 0 bs+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 -- | Number of incompatible splits. Similar to 'symmetricWith' but all--- bipartition induced by multifurcations are considered. For a detailed--- description of how the distance is calculated, see 'compatible'.+-- bipartitions induced by multifurcations are considered. For a detailed+-- description of how the distance is calculated, see+-- 'ELynx.Data.Tree.Bipartition.bpcompatible'. -- -- A multifurcation on a tree may (but not necessarily does) represent missing -- information about the order of bifurcations. In this case, it is interesting@@ -90,30 +114,34 @@ -- > AD|BC -- -- For an explanation of how compatibility of a bipartition with a--- multipartition is checked, see 'compatible'. Before using 'compatible',--- bipartitions are simply converted to multipartitions with two subsets.+-- multipartition is checked, see 'ELynx.Data.Tree.Multipartition.mpcompatible'.+-- Before using 'mpcompatible', bipartitions are simply converted to+-- multipartitions with two subsets. -- -- Only if a bipartition is not compatible with all induced multifurcations of -- the other tree, it is incompatible. -- -- 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 = countIncompatibilities putIncBs1 ms2 +- countIncompatibilities putIncBs2 ms1- where- -- Bipartitions.- bs t = bipartitions $ fmap f t- bs1 = bs t1- bs2 = bs 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+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 'id' for comparisons.+-- | See 'incompatibleSplitsWith', use 'getName' for comparisons. incompatibleSplits :: (Ord a, Named a) => Tree a -> Tree a -> Int incompatibleSplits = incompatibleSplitsWith getName @@ -121,36 +149,49 @@ -- 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.+-- nodes with a given function. Assumes that the trees are UNROOTED. See+-- http://evolution.genetics.washington.edu/phylip/doc/treedist.html. -- -- 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)+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 = 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+ -> 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 'id' for comparisons.+-- | 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 ]+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+ :: (Tree a -> Tree a -> b) -- ^ Distance function+ -> [Tree a] -- ^ Input trees+ -> [b] adjacent dist trs = [ dist x y | (x, y) <- zip trs (tail trs) ]
− src/ELynx/Data/Tree/EvoTree.hs
@@ -1,44 +0,0 @@-{- |-Module : ELynx.Data.Tree.EvoTree-Description : Evolutionary nodes-Copyright : (c) Dominik Schrempf 2019-License : GPL-3--Maintainer : dominik.schrempf@gmail.com-Stability : unstable-Portability : portable--Creation date: Thu Jan 17 14:19:26 2019.--XXX: This module is not used.---}--module ELynx.Data.Tree.EvoTree- ( EvoLabel (..)- ) where---- | An evolutionary label has some information about where the corresponding--- node is on the tree, and if the node is 'extant', 'extinct', 'internal', or--- 'external'. The latter two could also be determined from the tree. This could--- be species, genes or individuals; probably more.-class EvoLabel n where- extant :: n -> Bool- extinct :: n -> Bool-- internal :: n -> Bool- internal n = not $ extant n || extinct n-- external :: n -> Bool- external = not . internal---- -- | Glue branches together, so that one new tree emerges. It's root node is--- -- new, the sub-forest has to be given (a list of trees).--- glue :: (NodeType c)--- => PhyloLabel a b c -- ^ New root node.--- -> [PhyloTree a b c] -- ^ Sub-forest.--- -> PhyloTree a b c--- glue s@(PhyloLabel _ _ n) ts--- | extant n = error "Root node cannot be of type 'Exant'."--- | extinct n = error "Root node cannot be of type 'Extinct'."--- | otherwise = Node s ts
src/ELynx/Data/Tree/MeasurableTree.hs view
@@ -3,8 +3,8 @@ {- | Module : ELynx.Data.Tree.MeasurableTree Description : Functions on trees with branch lengths-Copyright : (c) Dominik Schrempf 2019-License : GPL-3+Copyright : (c) Dominik Schrempf 2020+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -16,7 +16,8 @@ module ELynx.Data.Tree.MeasurableTree- ( Measurable (..)+ ( Measurable(..)+ , extendBranchLength , distancesRootLeaves , distancesOriginLeaves , averageDistanceOriginLeaves@@ -26,18 +27,20 @@ , shortenStem , summarize , totalBranchLength- , normalize+ , normalizeBranchLength , prune , removeMultifurcations , ultrametric- ) where+ )+where -import qualified Data.ByteString.Lazy.Char8 as L+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.Equality (allNearlyEqual)+import ELynx.Tools -- | A 'Node' label with measurable and modifiable branch length to the parent. class Measurable a where@@ -46,10 +49,15 @@ -- | Set attached branch length. setLen :: Double -> a -> a --- | Elongate branch length. 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)@@ -58,7 +66,8 @@ -- 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]+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 -@@ -77,8 +86,9 @@ -- 'distancesOriginLeaves'. averageDistanceOriginLeaves :: (Measurable a) => Tree a -> Double averageDistanceOriginLeaves tr = sum ds / fromIntegral n- where ds = distancesOriginLeaves tr- n = length ds+ where+ ds = distancesOriginLeaves tr+ n = length ds -- | Height (max distance between origin and leaves) of a tree. Return 0 if the -- tree is empty.@@ -97,40 +107,50 @@ 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" $ map L.pack- [ "Leaves: " ++ show n ++ "."- , "Height: " ++ show h ++ "."- , "Average distance root to leaves: " ++ show h' ++ "."- , "Total branch length: " ++ show b ++ "." ]- where n = length . leaves $ t- h = height t- b = totalBranchLength t- h' = sum (distancesOriginLeaves t) / fromIntegral n+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.-normalize :: (Measurable a) => Tree a -> Tree a-normalize t = fmap (\n -> setLen (getLen n / s) n) t+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+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]+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.
src/ELynx/Data/Tree/Multipartition.hs view
@@ -1,8 +1,8 @@ {- | Module : ELynx.Data.Tree.Multipartition Description : Multipartitions on rose trees-Copyright : (c) Dominik Schrempf 2019-License : GPL-3+Copyright : (c) Dominik Schrempf 2020+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -16,30 +16,37 @@ -} module ELynx.Data.Tree.Multipartition- ( -- * The 'Multipartition' data type.- Multipartition ()+ (+ -- * The 'Multipartition' data type.+ Multipartition() , mps , mp , mpmap , mphuman , fromBipartition- -- * Working with 'Multipartition's.+ -- * Working with 'Multipartition's. , multipartitions , findSubset- , compatible- ) where+ , mpcompatible+ )+where -import Data.List (filter, find, foldl', intercalate)-import Data.Maybe (fromMaybe)-import qualified Data.Set as S+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)-import ELynx.Data.Tree.Subset+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 'Subset's, see 'ELynx.Data.Tree.Bipartition'. In a similar way, each+-- 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@@ -47,29 +54,26 @@ -- 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 (Subset a) -- ^ Set of partitions+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)) ++ ")"+mphuman f (Multipartition xs) =+ "(" ++ intercalate "|" (map (sshow f) (S.toList xs)) ++ ")" -- | Create a multipartition.-mp :: Ord a => [Subset a] -> Multipartition a-mp = mp' . filter (not . snull)+mp :: Ord a => [S.Set a] -> Multipartition a+mp = mp' . filter (not . S.null) -mp' :: Ord a => [Subset a] -> Multipartition a--- XXX: For now also allow multipartitions with no, one, or two elements.--- mp' [] = error "mp': Cannot create multipartition from empty list."--- mp' [_] = error "mp': Cannot create multipartition from list with one element."--- mp' [_, _] = error "mp': Cannot create multipartition from list with two elements."-mp' xs = Multipartition (S.fromList xs)+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 (smap f) xs+mpmap f (Multipartition xs) = Multipartition $ S.map (S.map f) xs instance (Eq a) => Eq (Multipartition a) where Multipartition xs == Multipartition ys = xs == ys@@ -79,47 +83,46 @@ -- | Convert bipartition to multipartition. fromBipartition :: Ord a => Bipartition a -> Multipartition a-fromBipartition bp = mp [l, r]- where (l, r) = bps bp+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 sempty (partitionTree t)- else error "multipartitions: The tree contains duplicate leaves."+ 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 => Subset a -> Tree (Subset 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 (sunion xs (rootLabel x)) y- r = multipartitionsUnsafe (sunion xs (rootLabel y)) x-multipartitionsUnsafe xs t@(Node _ ys ) = S.unions $- S.singleton (mp (xs : map rootLabel ys))- : zipWith multipartitionsUnsafe lvsOthers ys+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 multipartition containing a given element.-findSubset :: Ord a => a -> Multipartition a -> Subset a-findSubset l m = fromMaybe- -- Return the empty subset if nothing is found. This corresponds+-- | 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.- sempty- (find (smember l) ss)+ 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 (Subset a) -> a -> S.Set (Subset a)-addSubset m ss l = if not $ snull ss'- then ss' `S.insert` ss- else ss- where ss' = findSubset l m+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 -> Subset a -> S.Set (Subset a)+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@@ -134,7 +137,7 @@ -- -- 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 (Subset a)@. One could+-- 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.@@ -147,16 +150,17 @@ -- to 3. -- -- See also 'ELynx.Data.Tree.Bipartition.compatible'.-compatible :: (Ord a, Show a) => Multipartition a -> Multipartition a -> Bool--- compatible l r = traceShow lOverlaps $ traceShow rOverlaps $ and $-compatible 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+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/NamedTree.hs view
@@ -1,8 +1,8 @@ {- | Module : ELynx.Data.Tree.NamedTree Description : Trees with named nodes-Copyright : (c) Dominik Schrempf 2019-License : GPL-3+Copyright : (c) Dominik Schrempf 2020+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -13,11 +13,12 @@ -} module ELynx.Data.Tree.NamedTree- ( Named (..)- ) where+ ( Named(..)+ )+where -import qualified Data.ByteString.Lazy.Builder as L-import qualified Data.ByteString.Lazy.Char8 as L+import qualified Data.ByteString.Lazy.Builder as L+import qualified Data.ByteString.Lazy.Char8 as L -- | Data types with names. class Named a where
src/ELynx/Data/Tree/PhyloTree.hs view
@@ -3,8 +3,8 @@ {- | Module : ELynx.Data.Tree.PhyloTree Description : Phylogenetic trees-Copyright : (c) Dominik Schrempf 2019-License : GPL-3+Copyright : (c) Dominik Schrempf 2020+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -19,14 +19,15 @@ module ELynx.Data.Tree.PhyloTree- ( PhyloLabel (..)+ ( PhyloLabel(..) , removeBrInfo- ) where+ )+where import Data.Function-import Data.Maybe (fromMaybe)+import Data.Maybe ( fromMaybe ) import Data.Tree-import Test.QuickCheck hiding (label)+-- import Test.QuickCheck hiding ( label ) import ELynx.Data.Tree.BranchSupportTree import ELynx.Data.Tree.MeasurableTree@@ -47,26 +48,15 @@ instance Measurable (PhyloLabel a) where getLen = fromMaybe 0 . brLen setLen l x- | l >= 0 = x {brLen = Just l}+ | 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 Nothing l = l { brSup = Nothing } setBranchSupport (Just s) l- | s > 0 = l {brSup = Just s}+ | s > 0 = l { brSup = Just s } | otherwise = error "Branch support cannot be negative."---- 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) ) instance Named a => Named (PhyloLabel a) where getName = getName . label
+ src/ELynx/Data/Tree/SubSample.hs view
@@ -0,0 +1,64 @@+{- |+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/Subset.hs
@@ -1,105 +0,0 @@-{-# LANGUAGE GeneralizedNewtypeDeriving #-}--{- |-Module : ELynx.Data.Tree.Subset-Description : A subset is a set of elements-Copyright : (c) Dominik Schrempf 2019-License : GPL-3--Maintainer : dominik.schrempf@gmail.com-Stability : unstable-Portability : portable--Creation date: Fri Dec 13 11:02:43 2019.---}--module ELynx.Data.Tree.Subset- ( Subset ()- , sfromset- , sfromlist- , smap- , snull- , sempty- , ssingleton- , sunion- , sunions- , sdifference- , sintersection- , sdisjoint- , smember- , sshow- ) where--import Data.List (intercalate)-import qualified Data.Set as S---- | A 'Subset' is a set of elements of type a. For example, on phylogenetic--- trees, a 'Subset' is a set of leaves. In this case, a 'Subset' is induced,--- for example, by a node on the (rooted) tree. The 'Subset's of leaf nodes are--- singletons. The 'Subset's of the root node is the set of all leaves.--- 'Subset's are the building blocks of partitions. Each branch on the tree--- induces a bipartition, or a pair of 'Subset's, see--- 'ELynx.Data.Tree.Bipartition'. Multifurcations induce multipartitions, see--- 'ELynx.Data.Tree.Multipartition'.------ Internally, a subset is just an 'S.Set', since the order of elements--- within the subset is not important, but the uniqueness of elements is.-newtype Subset a = SS {subset :: S.Set a}- deriving (Show, Read, Eq, Ord, Semigroup, Monoid, Foldable)---- | Create a subset from a set.-sfromset :: S.Set a -> Subset a-sfromset = SS---- | Create a subset from a list. Throws an error if duplicate elements are--- present in the list.-sfromlist :: Ord a => [a] -> Subset a-sfromlist l = if S.size s == length l- then sfromset s- else error "pfromlist: List contains duplicate elements."- where s = S.fromList l---- | Map a function over all elements in a subset.-smap :: Ord b => (a -> b) -> Subset a -> Subset b-smap f = SS . S.map f . subset---- | Is the subset empty?-snull :: Subset a -> Bool-snull = S.null . subset---- | The empty subset.-sempty :: Subset a-sempty = SS S.empty---- | A subset with one element.-ssingleton :: a -> Subset a-ssingleton = SS . S.singleton---- | Unite two subsets.-sunion :: Ord a => Subset a -> Subset a -> Subset a-sunion p q = SS $ subset q `S.union` subset p---- | Unite a list of subsets.-sunions :: Ord a => [Subset a] -> Subset a-sunions = SS . S.unions . map subset---- | Difference of two subsets.-sdifference :: Ord a => Subset a -> Subset a -> Subset a-sdifference p q = SS $ subset p S.\\ subset q---- | Intersection of two subsets.-sintersection :: Ord a => Subset a -> Subset a -> Subset a-sintersection p q = SS $ S.intersection (subset p) (subset q)---- | Are two subsets disjoint?-sdisjoint :: Ord a => Subset a -> Subset a -> Bool-sdisjoint p q = S.disjoint (subset p) (subset q)---- | Check if an element is member of a subset.-smember :: Ord a => a -> Subset a -> Bool-smember x = S.member x . subset---- | Show the elements of a subset in a human readable way.-sshow :: (a -> String) -> Subset a -> String-sshow f = intercalate "," . map f . S.toList . subset
src/ELynx/Data/Tree/SumStat.hs view
@@ -2,7 +2,7 @@ Module : ELynx.Data.Tree.SumStat Description : Summary statistics for phylogenetic trees Copyright : (c) Dominik Schrempf 2018-License : GPL-3+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -10,6 +10,9 @@ 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@@ -17,11 +20,11 @@ , NChildSumStat , toNChildSumStat , formatNChildSumStat- ) where+ )+where -import qualified Data.ByteString.Builder as L-import qualified Data.ByteString.Lazy.Char8 as L-import Data.Monoid ((<>))+import qualified Data.ByteString.Builder as L+import qualified Data.ByteString.Lazy.Char8 as L import Data.Tree import ELynx.Data.Tree.MeasurableTree @@ -41,17 +44,17 @@ -- nLeaves2 branchLength2 -- .... formatNChildSumStat :: NChildSumStat -> L.ByteString-formatNChildSumStat s = L.toLazyByteString . mconcat $ map formatNChildSumStatLine s+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'+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+ where+ nChSS = map toNChildSumStat ts+ sumCh = sum $ map (snd . head) nChSS
src/ELynx/Data/Tree/Tree.hs view
@@ -1,8 +1,8 @@ {- | Module : ELynx.Data.Tree.Tree Description : Functions related to phylogenetic trees-Copyright : (c) Dominik Schrempf 2019-License : GPL-3+Copyright : (c) Dominik Schrempf 2020+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -10,6 +10,10 @@ 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@@ -33,6 +37,11 @@ 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. @@ -43,31 +52,30 @@ ( singleton , degree , leaves- , subTree- , subSample- , nSubSamples , pruneWith+ , dropLeafWith+ , intersectWith , merge , tZipWith , partitionTree , subForestGetSubsets+ , subTree , bifurcating , roots , connect , clades- ) where+ )+where -import Control.Monad-import Control.Monad.Primitive-import Data.Maybe-import qualified Data.Sequence as Seq-import qualified Data.Set as Set+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-import System.Random.MWC -import ELynx.Data.Tree.Subset-import ELynx.Tools.Random -- | The simplest tree. Usually an extant leaf. singleton :: a -> Tree a@@ -80,65 +88,65 @@ -- | Get leaves of tree. leaves :: Tree a -> [a] leaves (Node l []) = [l]-leaves (Node _ f) = concatMap leaves f+leaves (Node _ f ) = concatMap leaves f --- -- | Check if ancestor and daughters of first tree are a subset of the ancestor--- -- and daughters of the second tree. Useful to test if, e.g., speciations agree.--- rootNodesAgreeWith :: (Ord c) => (a -> c) -> Tree a -> (b -> c) -> Tree b -> Bool--- rootNodesAgreeWith f s g t =--- f (rootLabel s) == g (rootLabel t) &&--- S.fromList sDs `S.isSubSetOf` S.fromList tDs--- where sDs = map (f . rootLabel) (subForest s)--- tDs = map (g . rootLabel) (subForest t)+-- | 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) --- | 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+-- | 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 --- XXX: If module gets too big, move the sampling functions into their own--- module.--- | Extract a random sub tree 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 we have--- fast sub-sampling as well as lookup and don't have to recompute them 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 '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) --- | 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+-- | 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 --- | Prune or remove degree 2 inner nodes. The information stored in a pruned--- node can be used to change the daughter node. To discard this information,--- use, @pruneWith const tree@, otherwise @pruneWith (\daughter parent ->--- combined) tree@.-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)+-- 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 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 :)).+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 @@ -146,34 +154,48 @@ -- 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))+ 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 (Subset a)-partitionTree (Node l []) = Node (ssingleton l) []-partitionTree (Node _ xs) = Node (sunions $ map rootLabel xs') xs'+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)- => Subset a -- ^ Complementary partition at the stem- -> Tree (Subset a) -- ^ Tree with partition nodes- -> [Subset a]+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 = [ sunions $ lvs- : take i lvsChildren ++ drop (i+1) lvsChildren- | i <- [0 .. (nChildren - 1)] ]- lvsOthers = map (sunion lvs) lvsOtherChildren+ 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 and does not include degree two nodes. I--- know, one should use a proper data structure to encode bifurcating trees, but--- I don't have enough time for this now.+-- | 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@@ -197,35 +219,38 @@ 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 _ [] ) = [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."+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."+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 in all possible ways.+-- | 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.@@ -233,13 +258,11 @@ -- 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]+connect n l r = [ Node n [x, y] | x <- roots l, y <- roots r ] -- | Get clades induced by multifurcations.------ XXX: Probably introduce a new module defining a Clade.-clades :: Ord a => Tree a -> [Subset a]-clades (Node _ [] ) = []-clades (Node _ [x]) = clades x+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 = sfromlist (leaves t) : concatMap clades (subForest t)+clades t = S.fromList (leaves t) : concatMap clades (subForest t)
src/ELynx/Distribution/BirthDeath.hs view
@@ -5,7 +5,7 @@ Module : ELynx.Distribution.BirthDeath Description : Birth and death distribution Copyright : (c) Dominik Schrempf 2018- License : GPL-3+ License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -26,11 +26,14 @@ , cumulative , density , quantile- ) where+ )+where -import Data.Data (Data, Typeable)-import GHC.Generics (Generic)-import qualified Statistics.Distribution as D+import Data.Data ( Data+ , Typeable+ )+import GHC.Generics ( Generic )+import qualified Statistics.Distribution as D import ELynx.Distribution.Types @@ -43,17 +46,17 @@ } deriving (Eq, Typeable, Data, Generic) instance D.Distribution BirthDeathDistribution where- cumulative = cumulative+ 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@@ -61,24 +64,29 @@ -- | 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 [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)) instance D.ContGen BirthDeathDistribution where genContVar = D.genContinuous
src/ELynx/Distribution/BirthDeathCritical.hs view
@@ -5,7 +5,7 @@ Module : ELynx.Distribution.BirthDeathCritical Description : Birth and death distribution Copyright : (c) Dominik Schrempf 2018- License : GPL-3+ License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -27,11 +27,14 @@ , cumulative , density , quantile- ) where+ )+where -import Data.Data (Data, Typeable)-import GHC.Generics (Generic)-import qualified Statistics.Distribution as D+import Data.Data ( Data+ , Typeable+ )+import GHC.Generics ( Generic )+import qualified Statistics.Distribution as D import ELynx.Distribution.Types @@ -43,14 +46,13 @@ } deriving (Eq, Typeable, Data, Generic) instance D.Distribution BirthDeathCriticalDistribution where- cumulative = cumulative+ 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@@ -58,19 +60,22 @@ -- | 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
@@ -5,7 +5,7 @@ Module : ELynx.Distribution.BirthDeathCriticalNoTime Description : Birth and death distribution Copyright : (c) Dominik Schrempf 2018- License : GPL-3+ License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -27,11 +27,14 @@ , cumulative , density , quantile- ) where+ )+where -import Data.Data (Data, Typeable)-import GHC.Generics (Generic)-import qualified Statistics.Distribution as D+import Data.Data ( Data+ , Typeable+ )+import GHC.Generics ( Generic )+import qualified Statistics.Distribution as D import ELynx.Distribution.Types @@ -42,13 +45,12 @@ } deriving (Eq, Typeable, Data, Generic) instance D.Distribution BirthDeathCriticalNoTimeDistribution where- cumulative = cumulative+ 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@@ -57,17 +59,17 @@ -- | 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
@@ -5,7 +5,7 @@ Module : ELynx.Distribution.BirthDeathNearlyCritical Description : Birth and death distribution Copyright : (c) Dominik Schrempf 2018- License : GPL-3+ License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -29,11 +29,14 @@ , cumulative , density , quantile- ) where+ )+where -import Data.Data (Data, Typeable)-import GHC.Generics (Generic)-import qualified Statistics.Distribution as D+import Data.Data ( Data+ , Typeable+ )+import GHC.Generics ( Generic )+import qualified Statistics.Distribution as D import ELynx.Distribution.Types @@ -46,16 +49,16 @@ } deriving (Eq, Typeable, Data, Generic) instance D.Distribution BirthDeathNearlyCriticalDistribution where- cumulative = cumulative+ 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@@ -63,27 +66,41 @@ -- | 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 ++ "."+ | p >= 0 && p <= 1+ = res+ | otherwise+ = error+ $ "PointProcess.quantile: p must be in [0,1] range. Got: "+ ++ show p+ ++ "." where- den = l*(-3.0 + 2.0*t*(-1.0+p)*l)+m- t1 = (2.0 + t*(l - 4.0*p*l + m)) / den- t2Nom = 4.0 + t*(l*(4.0 + t*l + 8.0*p*(1.0 + t*l)) + 2.0*(2.0 + t*l - 4.0*p*(1.0 + t*l))*m + t*m*m)- t2 = t2Nom / (den**2)- res = 0.5 * (t1 + sqrt t2)+ den = l * (-3.0 + 2.0 * t * (-1.0 + p) * l) + m+ t1 = (2.0 + t * (l - 4.0 * p * l + m)) / den+ t2Nom =+ 4.0+ + t+ * ( l+ * (4.0 + t * l + 8.0 * p * (1.0 + t * l))+ + 2.0+ * (2.0 + t * l - 4.0 * p * (1.0 + t * l))+ * m+ + t+ * m+ * m+ )+ t2 = t2Nom / (den ** 2)+ res = 0.5 * (t1 + sqrt t2) instance D.ContGen BirthDeathNearlyCriticalDistribution where genContVar = D.genContinuous
src/ELynx/Distribution/CoalescentContinuous.hs view
@@ -2,7 +2,7 @@ Module : ELynx.Distribution.CoalescentContinuous Description : Distribution of coalescent times Copyright : (c) Dominik Schrempf 2018-License : GPL-3+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -14,13 +14,15 @@ module ELynx.Distribution.CoalescentContinuous ( coalescentDistributionCont- ) where+ )+where -import Numeric.SpecFunctions (choose)+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+ :: Int -- ^ Sample size.+ -> ExponentialDistribution coalescentDistributionCont n = exponential (choose n 2)
src/ELynx/Distribution/TimeOfOrigin.hs view
@@ -5,7 +5,7 @@ Module : ELynx.Distribution.TimeOfOrigin Description : Distribution of time of origin for birth and death trees Copyright : (c) Dominik Schrempf 2018- License : GPL-3+ License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -26,11 +26,14 @@ , cumulative , density , quantile- ) where+ )+where -import Data.Data (Data, Typeable)-import GHC.Generics (Generic)-import qualified Statistics.Distribution as D+import Data.Data ( Data+ , Typeable+ )+import GHC.Generics ( Generic )+import qualified Statistics.Distribution as D import ELynx.Distribution.Types @@ -43,15 +46,15 @@ } deriving (Eq, Typeable, Data, Generic) instance D.Distribution TimeOfOriginDistribution where- cumulative = cumulative+ 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@@ -61,24 +64,30 @@ density :: TimeOfOriginDistribution -> Time -> Double density (TOD nn l m) x | x < 0 = 0- | otherwise = n * l**n * d**2 * t1**(n-1.0) * ex / (t2**(n+1.0))- where d = l - m- n = fromIntegral nn- ex = exp(-d*x)- t1 = 1.0 - ex- t2 = l - m*ex+ | otherwise = n * l ** n * d ** 2 * t1 ** (n - 1.0) * ex / (t2 ** (n + 1.0))+ where+ d = l - m+ n = fromIntegral nn+ ex = exp (-d * x)+ t1 = 1.0 - ex+ t2 = l - m * ex -- | The inverted cumulative probability distribution 'cumulative'. See also -- 'D.ContDistr'. quantile :: TimeOfOriginDistribution -> Double -> Time quantile (TOD n' l m) p- | p >= 0 && p <= 1 = -1.0/d * log(t1/t2)- | otherwise =- error $ "PointProcess.quantile: p must be in [0,1] range. Got: " ++ show p ++ "."- where d = l - m- n = fromIntegral n'- t1 = l*(1.0-p**(1.0/n))- t2 = l - p**(1.0/n)*m+ | 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
@@ -5,7 +5,7 @@ Module : ELynx.Distribution.TimeOfOriginNearCritical Description : Distribution of time of origin for birth and death trees Copyright : (c) Dominik Schrempf 2018- License : GPL-3+ License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -26,11 +26,14 @@ , cumulative , density , quantile- ) where+ )+where -import Data.Data (Data, Typeable)-import GHC.Generics (Generic)-import qualified Statistics.Distribution as D+import Data.Data ( Data+ , Typeable+ )+import GHC.Generics ( Generic )+import qualified Statistics.Distribution as D import ELynx.Distribution.Types @@ -43,17 +46,17 @@ } deriving (Eq, Typeable, Data, Generic) instance D.Distribution TimeOfOriginNearCriticalDistribution where- cumulative = cumulative+ 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@@ -61,24 +64,30 @@ -- | 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
@@ -2,7 +2,7 @@ Module : ELynx.Distribution.Types Description : Data types for distributions on trees Copyright : (c) Dominik Schrempf 2018-License : GPL-3+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -16,7 +16,8 @@ module ELynx.Distribution.Types ( Time , Rate- ) where+ )+where -- | Branch lengths are measured in time. type Time = Double
src/ELynx/Export/Tree/Newick.hs view
@@ -1,8 +1,8 @@ {- | Module : ELynx.Export.Tree.Newick Description : Export tree objects to Newick format-Copyright : (c) Dominik Schrempf 2019-License : GPL-3+Copyright : (c) Dominik Schrempf 2020+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -10,7 +10,8 @@ Creation date: Thu Jan 17 13:51:47 2019. -Parts of the code are from https://hackage.haskell.org/package/BiobaseNewick.+Some functions are inspired by+[Biobase.Newick.Import](https://hackage.haskell.org/package/BiobaseNewick). See nomenclature in 'ELynx.Data.Tree.Tree'. @@ -18,40 +19,36 @@ module ELynx.Export.Tree.Newick ( toNewick- -- , toNewickPhyloIntTree- -- , toNewickPhyloByteStringTree- ) where+ )+where -import qualified Data.ByteString.Lazy.Builder as L-import qualified Data.ByteString.Lazy.Char8 as L-import Data.List (intersperse)--- import Data.Maybe+import qualified Data.ByteString.Lazy.Builder as L+import qualified Data.ByteString.Lazy.Char8 as L+import Data.List ( intersperse ) import Data.Tree --- import ELynx.Data.Tree.BranchSupportTree--- import ELynx.Data.Tree.MeasurableTree-import ELynx.Data.Tree.NamedTree-import ELynx.Data.Tree.PhyloTree-import ELynx.Tools.ByteString (c2w)+import ELynx.Data.Tree+import ELynx.Tools -- | 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 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
src/ELynx/Import/Tree/Newick.hs view
@@ -1,8 +1,10 @@+{-# LANGUAGE DeriveGeneric #-}+ {- | Module : ELynx.Import.Tree.Newick Description : Import Newick trees-Copyright : (c) Dominik Schrempf 2019-License : GPL-3+Copyright : (c) Dominik Schrempf 2020+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -10,18 +12,23 @@ Creation date: Thu Jan 17 14:56:27 2019. -Code partly taken from Biobase.Newick.Import.+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.+ -} module ELynx.Import.Tree.Newick ( Parser -- * Newick tree format+ , NewickFormat(..)+ , description , newick , oneNewick , manyNewick@@ -30,69 +37,108 @@ , node , name , branchLength- -- * Newick tree format with branch support as node names (e.g., used by IQ-TREE)- , newickIqTree- , oneNewickIqTree- , manyNewickIqTree- ) where+ )+where -import qualified Data.ByteString.Lazy as L+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 Text.Megaparsec.Byte.Lexer ( decimal+ , float+ ) -import ELynx.Data.Tree.PhyloTree-import ELynx.Tools.ByteString (c2w)+import ELynx.Data.Tree+import ELynx.Tools -- | Shortcut. type Parser = Parsec Void L.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+ deriving (Eq, Show, Read, Bounded, Enum, Generic)++instance FromJSON NewickFormat++instance ToJSON NewickFormat++-- | Short description of the supported Newick formats.+description :: NewickFormat -> String+description Standard =+ "Standard: Branch support values are stored in square brackets after branch lengths."+description IqTree+ = "IqTree: Branch support values are stored as node names after the closing bracket of forests."+description RevBayes+ = "RevBayes Key-value pairs is provided in square brackets after node names as well as branch lengths. XXX: Key value pairs are IGNORED at the moment."+ -- | Parse a single Newick tree. Also succeeds when more trees follow.-newick :: Parser (Tree (PhyloLabel L.ByteString))-newick = tree <* char (c2w ';') <?> "newick"+newick :: NewickFormat -> Parser (Tree (PhyloLabel L.ByteString))+newick Standard = newickStandard+newick IqTree = newickIqTree+newick RevBayes = newickRevBayes -- | Parse a single Newick tree. Fails when end of file is not reached.-oneNewick :: Parser (Tree (PhyloLabel L.ByteString))-oneNewick = newick <* space <* eof <?> "oneNewick"+oneNewick :: NewickFormat -> Parser (Tree (PhyloLabel L.ByteString))+oneNewick Standard = oneNewickStandard+oneNewick IqTree = oneNewickIqTree+oneNewick RevBayes = oneNewickRevBayes -- | Parse many Newick trees until end of file.-manyNewick :: Parser [Tree (PhyloLabel L.ByteString)]-manyNewick = some (newick <* space) <* eof <?> "manyNewick"+manyNewick :: NewickFormat -> Parser [Tree (PhyloLabel L.ByteString)]+manyNewick Standard = manyNewickStandard+manyNewick IqTree = manyNewickIqTree+manyNewick RevBayes = manyNewickRevBayes +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. Fails when end of file is not reached.+oneNewickStandard :: Parser (Tree (PhyloLabel L.ByteString))+oneNewickStandard = newickStandard <* eof <?> "oneNewick"++-- | Parse many Newick trees until end of file.+manyNewickStandard :: Parser [Tree (PhyloLabel L.ByteString)]+manyNewickStandard = some newickStandard <* eof <?> "manyNewick"+ tree :: Parser (Tree (PhyloLabel L.ByteString))-tree = space *> (branched <|> leaf) <?> "tree"+tree = branched <|> leaf <?> "tree" branched :: Parser (Tree (PhyloLabel L.ByteString)) branched = do f <- forest- n <- node- <?> "branched"+ 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 = do- _ <- char (c2w '(')- f <- tree `sepBy1` char (c2w ',')- _ <- char (c2w ')')- <?> "forest"- return f+forest = between (w '(') (w ')') (tree `sepBy1` w ',') <?> "forest" +-- TODO: Why try? branchSupport :: Parser (Maybe Double) branchSupport = optional $ do- _ <- try $ char (c2w '[')+ _ <- try $ w '[' s <- try float <|> try decimalAsDouble- _ <- try $ char (c2w ']')+ _ <- try $ w ']' return s -- | A 'leaf' is a 'node' without children. leaf :: Parser (Tree (PhyloLabel L.ByteString)) leaf = do- n <- node- <?> "leaf"+ n <- node <?> "leaf" return $ Node n [] -- | A 'node' has a name and a 'branchLength'.@@ -100,8 +146,7 @@ node = do n <- name b <- branchLength- s <- branchSupport- <?> "node"+ s <- branchSupport <?> "node" return $ PhyloLabel n s b checkNameCharacter :: Word8 -> Bool@@ -114,7 +159,7 @@ -- | Branch length. branchLength :: Parser (Maybe Double)-branchLength = (optional $ char (c2w ':') *> branchLengthGiven) <?> "branchLength"+branchLength = optional (w ':' *> branchLengthGiven) <?> "branchLength" branchLengthGiven :: Parser Double branchLengthGiven = try float <|> decimalAsDouble@@ -123,44 +168,45 @@ decimalAsDouble = fromIntegral <$> (decimal :: Parser Int) ----------------------------------------------------------------------------------- IQ-TREE STUFF.+-- IQ-TREE. --- | IQ-TREE stores the branch support as node names after the closing bracket of a forest.+-- 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 = treeIqTree <* char (c2w ';') <?> "newickIqTree"+newickIqTree = space *> treeIqTree <* w ';' <* space <?> "newickIqTree" --- | IQ-TREE stores the branch support as node names after the closing bracket of a forest.+-- See 'newickIqTree'. Parse a single Newick tree. Fails when end of file is not+-- reached. oneNewickIqTree :: Parser (Tree (PhyloLabel L.ByteString))-oneNewickIqTree = newickIqTree <* space <* eof <?> "oneNewickIqTree"+oneNewickIqTree = newickIqTree <* eof <?> "oneNewickIqTree" --- | IQ-TREE stores the branch support as node names after the closing bracket of a forest.+-- See 'newickIqTree'. Parse many Newick trees until end of file. manyNewickIqTree :: Parser [Tree (PhyloLabel L.ByteString)]-manyNewickIqTree = some (newickIqTree <* space) <* eof <?> "manyNewickIqTree"+manyNewickIqTree = some newickIqTree <* eof <?> "manyNewickIqTree" -- IQ-TREE stores the branch support as node names after the closing bracket of a forest. treeIqTree :: Parser (Tree (PhyloLabel L.ByteString))-treeIqTree = space *> (branchedIqTree <|> leaf) <?> "treeIqTree"---- IQ-TREE stores the branch support as node names after the closing bracket of a forest.-forestIqTree :: Parser [Tree (PhyloLabel L.ByteString)]-forestIqTree = do- _ <- char (c2w '(')- f <- treeIqTree `sepBy1` char (c2w ',')- _ <- char (c2w ')')- <?> "forestIqTree"- return f+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 f <- forestIqTree s <- branchSupportIqTree- n <- nodeIqTree- <?> "branchedIqTree"- let n' = n {brSup = s}+ n <- nodeIqTree <?> "branchedIqTree"+ let n' = n { brSup = s } return $ Node 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+ _ <- w '('+ f <- treeIqTree `sepBy1` w ','+ _ <- w ')' <?> "forestIqTree"+ 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 @@ -168,6 +214,64 @@ nodeIqTree :: Parser (PhyloLabel L.ByteString) nodeIqTree = do n <- name+ b <- branchLength <?> "nodeIqTree"+ return $ PhyloLabel n Nothing b+++--------------------------------------------------------------------------------+-- RevBayes.++-- RevBayes uses square brackets and key-value pairs to define information+-- about nodes and branches. Parse a single Newick tree. Also succeeds when more+-- trees follow.+--+-- XXX: Key value pairs are IGNORED at the moment.+newickRevBayes :: Parser (Tree (PhyloLabel L.ByteString))+newickRevBayes =+ space *> 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 = newickRevBayes <* eof <?> "oneNewickRevBayes"++-- See 'newickRevBayes'. Parse many Newick trees until end of file.+manyNewickRevBayes :: Parser [Tree (PhyloLabel L.ByteString)]+manyNewickRevBayes = some newickRevBayes <* eof <?> "manyNewickRevBayes"++treeRevBayes :: Parser (Tree (PhyloLabel L.ByteString))+treeRevBayes = branchedRevBayes <|> leafRevBayes <?> "treeRevBayes"++branchedRevBayes :: Parser (Tree (PhyloLabel L.ByteString))+branchedRevBayes = do+ f <- forestRevBayes+ n <- nodeRevBayes <?> "branchedRevBayes"+ return $ Node n f++forestRevBayes :: Parser [Tree (PhyloLabel L.ByteString)]+forestRevBayes = do+ _ <- w '('+ f <- treeRevBayes `sepBy1` w ','+ _ <- w ')' <?> "forestRevBayes"+ return f++nodeRevBayes :: Parser (PhyloLabel L.ByteString)+nodeRevBayes = do+ n <- name+ _ <- optional brackets b <- branchLength- <?> "nodeIqTree"+ _ <- optional brackets <?> "nodeRevBayes" return $ PhyloLabel n Nothing b++leafRevBayes :: Parser (Tree (PhyloLabel L.ByteString))+leafRevBayes = do+ n <- nodeRevBayes <?> "leafRevBayes"+ return $ Node n []++-- Drop anything between brackets.+brackets :: Parser ()+brackets = do+ _ <- between (w '[')+ (w ']')+ (takeWhileP (Just "allCharsButBracketEnd") (/= c2w ']'))+ return ()
src/ELynx/Simulate/Coalescent.hs view
@@ -2,7 +2,7 @@ Module : ELynx.Simulate.Coalescent Description : Generate coalescent trees Copyright : (c) Dominik Schrempf 2018-License : GPL-3+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -10,55 +10,53 @@ Creation date: Wed May 16 13:13:11 2018. -TODO: This is unused, is it?- -} module ELynx.Simulate.Coalescent ( simulate- ) where+ )+where import Control.Monad.Primitive import Data.Tree import Statistics.Distribution import System.Random.MWC -import ELynx.Data.Tree.MeasurableTree-import ELynx.Data.Tree.PhyloTree-import ELynx.Data.Tree.Tree+import ELynx.Data.Tree import ELynx.Distribution.CoalescentContinuous -- | 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)+ => Int -- ^ Number of leaves.+ -> Gen (PrimState m)+ -> m (Tree (PhyloLabel Int)) simulate n = simulate' n 0 trs- where trs = [ singleton (PhyloLabel i Nothing Nothing) | i <- [0..n-1] ]+ where trs = [ singleton (PhyloLabel i Nothing Nothing) | i <- [0 .. n - 1] ] -simulate' :: (PrimMonad m)- => Int- -> Int- -> [Tree (PhyloLabel Int)]- -> Gen (PrimState m)- -> m (Tree (PhyloLabel Int))+simulate'+ :: (PrimMonad m)+ => Int+ -> Int+ -> [Tree (PhyloLabel Int)]+ -> Gen (PrimState m)+ -> m (Tree (PhyloLabel 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+ | n <= 0 = error "Cannot construct trees without leaves."+ | n == 1 && length trs /= 1 = error "Too many trees provided."+ | n == 1 && length trs == 1 = return $ head trs+ | otherwise = do -- Indices of the leaves to join will be i-1 and i.- i <- uniformR (1, n-1) g- -- The time of the coalescent event.- t <- genContVar (coalescentDistributionCont n) g- let trs' = map (lengthenStem 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]- -- 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+ 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+ -- Join the two chosen trees.+ tm = Node (PhyloLabel a Nothing Nothing) [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
@@ -4,7 +4,7 @@ Module : ELynx.Simulate.PointProcess Description : Point process and functions Copyright : (c) Dominik Schrempf 2018- License : GPL-3+ License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -27,18 +27,18 @@ , toReconstructedTree , simulateReconstructedTree , simulateNReconstructedTrees- ) where+ )+where import Control.Monad import Control.Monad.Primitive-import Data.List (mapAccumL)+import Data.List ( mapAccumL ) import Data.Tree-import qualified Statistics.Distribution as D (genContVar)+import qualified Statistics.Distribution as D+ ( genContVar ) import System.Random.MWC -import ELynx.Data.Tree.MeasurableTree-import ELynx.Data.Tree.PhyloTree-import ELynx.Data.Tree.Tree+import ELynx.Data.Tree import ELynx.Distribution.BirthDeath import ELynx.Distribution.BirthDeathCritical import ELynx.Distribution.BirthDeathCriticalNoTime@@ -46,8 +46,7 @@ import ELynx.Distribution.TimeOfOrigin import ELynx.Distribution.TimeOfOriginNearCritical import ELynx.Distribution.Types-import ELynx.Tools.Equality-import ELynx.Tools.List+import ELynx.Tools epsNearCriticalPointProcess :: Double epsNearCriticalPointProcess = 1e-5@@ -73,38 +72,49 @@ -- | 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)+ => 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) -- 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."+ 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+ 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+ 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+ 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."+ | 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.@@ -113,36 +123,37 @@ -- 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+ (m =~= l) && not c = do+ !vs <- replicateM (n - 1) (D.genContVar (BDCD t l) g)+ return $ PointProcess [0 .. (n - 1)] vs t | (abs (m - l) <= epsNearCriticalPointProcess) && not c = do- !vs <- replicateM (n-1) (D.genContVar (BDNCD t l m) g)- return $ PointProcess [0..(n-1)] vs t+ !vs <- replicateM (n - 1) (D.genContVar (BDNCD t l m) g)+ return $ PointProcess [0 .. (n - 1)] vs t | not c = do- !vs <- replicateM (n-1) (D.genContVar (BDD t l m) g)- return $ PointProcess [0..(n-1)] vs t+ !vs <- replicateM (n - 1) (D.genContVar (BDD t l m) g)+ return $ PointProcess [0 .. (n - 1)] vs t | (m =~= l) && c = do- !vs <- replicateM (n-2) (D.genContVar (BDCD t l) g)- vs' <- randomInsert t vs g- return $ PointProcess [0..(n-1)] vs' t+ !vs <- replicateM (n - 2) (D.genContVar (BDCD t l) g)+ vs' <- randomInsertList t vs g+ return $ PointProcess [0 .. (n - 1)] vs' t | (abs (m - l) <= epsNearCriticalPointProcess) && c = do- !vs <- replicateM (n-2) (D.genContVar (BDNCD t l m) g)- vs' <- randomInsert t vs g- return $ PointProcess [0..(n-1)] vs' t+ !vs <- replicateM (n - 2) (D.genContVar (BDNCD t l m) g)+ vs' <- randomInsertList t vs g+ return $ PointProcess [0 .. (n - 1)] vs' t | c = do- !vs <- replicateM (n-2) (D.genContVar (BDD t l m) g)- vs' <- randomInsert t vs g- return $ PointProcess [0..(n-1)] vs' t+ !vs <- replicateM (n - 2) (D.genContVar (BDD t l m) g)+ vs' <- randomInsertList t vs g+ return $ PointProcess [0 .. (n - 1)] vs' t | otherwise = error "simulate: Fell through guard, this should never happen." -- | Sort the values of a point process and their indices to be (the indices -- that they will have while creating the tree). sort :: (Ord b) => PointProcess a b -> ([b], [Int]) sort (PointProcess _ vs _) = (vsSorted, isSorted)- where vsIsSorted = sortWithIndices 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]@@ -151,8 +162,7 @@ -- 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)+fAcc is i = (i : is, i') where i' = i - length (filter (< i) is) -- | See 'simulateReconstructedTree', but n times. simulateNReconstructedTrees@@ -179,7 +189,8 @@ -> Rate -- ^ Death rate -> Gen (PrimState m) -- ^ Generator (see 'System.Random.MWC') -> m (Tree (PhyloLabel Int))-simulateReconstructedTree n t l m g = toReconstructedTree 0 <$> simulate n t l m g+simulateReconstructedTree n t l m g =+ toReconstructedTree 0 <$> simulate n t l m g -- | Convert a point process to a reconstructed tree. See Lemma 2.2. @@ -191,45 +202,70 @@ -- 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 (PhyloLabel 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+ 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 -- 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.+ -> [Tree (PhyloLabel a)] -- Leaves with accumulated root branch lengths.+ -> [Double] -- Accumulated heights of the leaves.+ -> Tree (PhyloLabel a)+toReconstructedTree' [] [] _ trs _ = head trs 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+ !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
test/ELynx/Data/Tree/BipartitionSpec.hs view
@@ -2,8 +2,8 @@ {- | Module : ELynx.Data.Tree.BipartitionSpec-Copyright : (c) Dominik Schrempf 2019-License : GPL-3+Copyright : (c) Dominik Schrempf 2020+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -14,71 +14,76 @@ -} module ELynx.Data.Tree.BipartitionSpec- (spec- ) where+ ( spec+ )+where -import qualified Data.ByteString.Lazy as L-import qualified Data.Map as M+import qualified Data.ByteString.Lazy as L+import qualified Data.Map as M import Data.Monoid-import qualified Data.Set as S+import qualified Data.Set as S import Data.Tree import Test.Hspec -import ELynx.Data.Tree.Bipartition-import ELynx.Data.Tree.MeasurableTree-import ELynx.Data.Tree.PhyloTree-import ELynx.Data.Tree.Subset+import ELynx.Data.Tree import ELynx.Import.Tree.Newick-import ELynx.Tools.InputOutput+import ELynx.Tools -sfrom :: [L.ByteString] -> Subset L.ByteString-sfrom = sfromset . S.fromList+sfrom :: [L.ByteString] -> S.Set L.ByteString+sfrom = S.fromList treeFileSimple :: FilePath treeFileSimple = "data/TreeDist.trees" getSimpleTrees :: IO [Tree (PhyloLabel L.ByteString)]-getSimpleTrees = parseFileWith manyNewick treeFileSimple+getSimpleTrees = parseFileWith (manyNewick 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 = 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 })+ ] 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 = 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"])+ ] 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 = 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"])+ ] spec :: Spec spec = do- describe "bipartitions" $- it "calculates correct bipartitions for sample trees" $ do- simpleTrees <- map removeBrInfo <$> 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- simpleTrees <- getSimpleTrees- bipartitionToBranchLength label (Sum . getLen) (simpleTrees !! 2)- `shouldBe` bipartitionToBranchAnswer+ describe "bipartitions"+ $ it "calculates correct bipartitions for sample trees"+ $ do+ simpleTrees <- map removeBrInfo <$> 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+ simpleTrees <- getSimpleTrees+ bipartitionToBranchLength label (Sum . getLen) (simpleTrees !! 2)+ `shouldBe` bipartitionToBranchAnswer
test/ELynx/Data/Tree/DistanceSpec.hs view
@@ -1,9 +1,10 @@ {-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE OverloadedStrings #-} {- | Module : ELynx.Data.Tree.DistanceSpec-Copyright : (c) Dominik Schrempf 2019-License : GPL-3+Copyright : (c) Dominik Schrempf 2020+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -14,73 +15,190 @@ -} module ELynx.Data.Tree.DistanceSpec- (spec- ) where+ ( spec+ )+where -import qualified Data.ByteString.Lazy.Char8 as L+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.Distance-import ELynx.Data.Tree.NamedTree-import ELynx.Data.Tree.PhyloTree+import Test.QuickCheck.Instances.Containers+ ( )+import ELynx.Data.Tree+import ELynx.Data.Tree.PhyloTreeArbitraryInstance+ ( ) import ELynx.Import.Tree.Newick-import ELynx.Tools.Equality-import ELynx.Tools.InputOutput+import ELynx.Tools treeFileSimple :: FilePath treeFileSimple = "data/TreeDist.trees" getSimpleTrees :: IO [Tree (PhyloLabel L.ByteString)]-getSimpleTrees = parseFileWith manyNewick treeFileSimple+getSimpleTrees = parseFileWith (manyNewick Standard) treeFileSimple treeFileMany :: FilePath treeFileMany = "data/Many.trees" getManyTrees :: IO [Tree (PhyloLabel L.ByteString)]-getManyTrees = parseFileWith manyNewick treeFileMany+getManyTrees = parseFileWith (manyNewick 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+ :: (Num b, Eq b) => (Tree a -> Tree a -> b) -> Tree a -> Bool prop_dist_same_tree distanceMeasure t = distanceMeasure t t == 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 = 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 = 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 = Node+ ' '+ [ Node ' ' [Node 'A' [], Node 'C' []]+ , Node ' ' [Node 'B' [], Node ' ' [Node 'D' [], Node 'E' []]]+ ] +incSplitTree1a :: Tree (PhyloLabel L.ByteString)+incSplitTree1a = parseByteStringWith "" (oneNewick IqTree) "((a,b)0.7,(c,d));"++incSplitTree1b :: Tree (PhyloLabel L.ByteString)+incSplitTree1b = parseByteStringWith "" (oneNewick IqTree) "((a,b)0.7,c,d);"++incSplitTree2 :: Tree (PhyloLabel L.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));"++incSplitTree4 :: Tree (PhyloLabel L.ByteString)+incSplitTree4 = parseByteStringWith "" (oneNewick IqTree) "(((a,c),b),(d,e));"+ spec :: Spec spec = do describe "symmetric" $ do@@ -93,9 +211,13 @@ -- 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)+ it "is zero for a collection of random trees"+ $ property+ $ prop_dist_same_tree+ (symmetric :: Tree (PhyloLabel Double)+ -> Tree (PhyloLabel Double)+ -> Int+ ) describe "incompatibleSplit" $ do it "calculates correct distances for sample trees" $ do@@ -103,10 +225,22 @@ incompatibleSplits bifurcatingComp multifurcating `shouldBe` 0 incompatibleSplits bifurcatingIncomp multifurcating `shouldBe` 2 incompatibleSplits multifurcating bifurcatingIncomp `shouldBe` 2- it "is zero for a collection of random trees" $- property $ prop_dist_same_tree- (incompatibleSplits :: Tree (PhyloLabel Double) -> Tree (PhyloLabel Double) -> Int)-+ 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+ ) describe "branchScore" $ do it "calculates correct distances for sample trees" $ do@@ -114,6 +248,10 @@ -- 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)+ it "is zero for a collection of random trees"+ $ property+ $ prop_dist_same_tree+ (branchScore :: Tree (PhyloLabel Double)+ -> Tree (PhyloLabel Double)+ -> Double+ )
test/ELynx/Data/Tree/MultipartitionSpec.hs view
@@ -2,8 +2,8 @@ {- | Module : ELynx.Data.Tree.MultipartitionSpec-Copyright : (c) Dominik Schrempf 2019-License : GPL-3+Copyright : (c) Dominik Schrempf 2020+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -14,46 +14,50 @@ -} module ELynx.Data.Tree.MultipartitionSpec- (spec- ) where+ ( spec+ )+where -import Data.Set (Set, empty, fromList)+import Data.Set ( Set+ , empty+ , fromList+ ) import Data.Tree import Test.Hspec import Test.QuickCheck-import Test.QuickCheck.Instances.Containers ()+import Test.QuickCheck.Instances.Containers+ ( ) -import ELynx.Data.Tree.BranchSupportTree-import ELynx.Data.Tree.MeasurableTree-import ELynx.Data.Tree.Multipartition (Multipartition, mp,- multipartitions)-import ELynx.Data.Tree.NamedTree-import ELynx.Data.Tree.PhyloTree-import ELynx.Data.Tree.Subset (sfromlist)+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 [sfromlist [1], sfromlist [3], sfromlist [4,5,6]]- , mp [sfromlist [1,3], sfromlist [4], sfromlist [5], sfromlist [6]] ]+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 [sfromlist [1], sfromlist [2], sfromlist [3,4], sfromlist [5]] ]+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+ :: (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)+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/PhyloTreeArbitraryInstance.hs view
@@ -0,0 +1,33 @@+{- |+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/TreeSpec.hs view
@@ -2,8 +2,8 @@ {- | Module : ELynx.Data.Tree.TreeSpec-Copyright : (c) Dominik Schrempf 2019-License : GPL-3+Copyright : (c) Dominik Schrempf 2020+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -14,32 +14,28 @@ -} module ELynx.Data.Tree.TreeSpec- (spec- ) where+ ( spec+ )+where -import qualified Data.ByteString.Lazy.Char8 as L+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 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.Bipartition (bipartition, compatible)-import ELynx.Data.Tree.MeasurableTree (height, prune)-import ELynx.Data.Tree.NamedTree (getName)-import ELynx.Data.Tree.PhyloTree (PhyloLabel (PhyloLabel),- label)-import ELynx.Data.Tree.Subset (Subset, smap)-import ELynx.Data.Tree.Tree (bifurcating, clades,- connect, leaves,- pruneWith, roots,- subTree)-import ELynx.Import.Tree.Newick (manyNewick, newick,- oneNewick)-import ELynx.Tools.InputOutput (parseByteStringWith,- parseFileWith)+import ELynx.Data.Tree+import ELynx.Import.Tree.Newick+ hiding ( node )+import ELynx.Tools node :: Int -> Tree Int node n = Node n []@@ -54,53 +50,61 @@ 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);"+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 sampleTreeBS+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- -- XXX: Skip not bifurcating trees. This is ugly, I know.- | not $ bifurcating t = True- | length (leaves t) < 3 = length (roots t) == 1- | otherwise = length (roots t) == 2 * length (leaves t) - 3+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- -- XXX: Skip not bifurcating trees. This is ugly, I know.- | not (bifurcating l) || not (bifurcating r) = True+ | 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 = Subset a+type Constraint a = S.Set a compatibleAll :: (Show a, Ord a) => Tree a -> [Constraint a] -> Bool-compatibleAll (Node _ [l, r]) cs = all (compatible (bipartition l)) cs &&- all (compatible (bipartition r)) cs+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 (smap f) cs)+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+ 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+ 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@@ -109,29 +113,30 @@ 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+ 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" [] ] ]+ [ 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))+ modifyMaxSize (* 100)+ $ it "returns the correct number of rooted trees for arbitrary trees"+ $ property (prop_roots :: (Tree Int -> Bool)) - 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)+ -- 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 "data/ConnectA.tree"- b <- parseFileWith oneNewick "data/ConnectB.tree"- c <- parseFileWith manyNewick "data/ConnectConstraints.tree"- let ts = connect (PhyloLabel "" Nothing (Just 1.0)) a b- cs = concatMap clades c :: [Constraint (PhyloLabel L.ByteString)]+ 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` 63 length ts' `shouldBe` 15
test/ELynx/Export/Tree/NewickSpec.hs view
@@ -1,7 +1,7 @@ {- | Module : ELynx.Export.Tree.NewickSpec-Copyright : (c) Dominik Schrempf 2019-License : GPL-3+Copyright : (c) Dominik Schrempf 2020+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -12,41 +12,50 @@ -} module ELynx.Export.Tree.NewickSpec- (spec- ) where+ ( spec+ )+where -import qualified Data.ByteString.Lazy.Char8 as L+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 Test.QuickCheck hiding ( label )+import Test.QuickCheck.Instances.Containers+ ( ) -import ELynx.Data.Tree.NamedTree-import ELynx.Data.Tree.PhyloTree+import ELynx.Data.Tree+import ELynx.Data.Tree.PhyloTreeArbitraryInstance+ ( ) import ELynx.Export.Tree.Newick import ELynx.Import.Tree.Newick-import ELynx.Tools.InputOutput+import ELynx.Tools samplePhyloByteStringTree :: Tree (PhyloLabel L.ByteString)-samplePhyloByteStringTree =- parseByteStringWith "Sample newick byte string" newick sampleNewickByteString1+samplePhyloByteStringTree = parseByteStringWith "Sample newick byte string"+ (newick Standard)+ sampleNewickByteString1 sampleNewickByteString1 :: L.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);"+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 t = parseByteStringWith "Sample newick byte string" newick (toNewick t) ==+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+ fmap (\l -> l { label = getName l }) 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 bs `shouldBe` samplePhyloByteStringTree+ (newick Standard)+ bs+ `shouldBe` samplePhyloByteStringTree - it "should be an invariant for a sample tree" $- property prop_invariant+ it "should be an invariant for a sample tree" $ property prop_invariant
test/ELynx/Import/Tree/NewickSpec.hs view
@@ -1,7 +1,7 @@ {- | Module : ELynx.Import.Tree.NewickSpec-Copyright : (c) Dominik Schrempf 2019-License : GPL-3+Copyright : (c) Dominik Schrempf 2020+License : GPL-3.0-or-later Maintainer : dominik.schrempf@gmail.com Stability : unstable@@ -12,16 +12,20 @@ -} module ELynx.Import.Tree.NewickSpec- (spec) where+ ( spec+ )+where -import qualified Data.ByteString.Lazy.Char8 as L+import qualified Data.ByteString.Lazy.Char8 as L import Data.Tree-import ELynx.Data.Tree.PhyloTree-import ELynx.Import.Tree.Newick import Test.Hspec import Test.Hspec.Megaparsec import Text.Megaparsec +import ELynx.Tools+import ELynx.Data.Tree+import ELynx.Import.Tree.Newick+ sampleLabelByteString :: L.ByteString sampleLabelByteString = L.pack "name:0.3" @@ -33,83 +37,119 @@ 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 = []}- ]- }+ [ 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 = []+ }+ ]+ } ] sampleNewickByteString1 :: L.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);"+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 = 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;"+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 = 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 = 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 = []+ }+ ]+ } +sampleTreeNewickRevbayes :: L.ByteString+sampleTreeNewickRevbayes =+ L.pack "[&R](l[IDL]:0.3[KEYVALPAIRS],r[IDR]:0.4[KEYVALPARIS])[ID]:0.3;"++sampleNewickRevBayesFile :: String+sampleNewickRevBayesFile = "data/NewickRevBayes.tree"+ 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+ it "parses a colon and a branch length"+ $ parse branchLength "" (L.pack ":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 "" (L.pack "")+ `shouldParse` Nothing describe "name" $ do- it "parses a string of printable characters" $- parse name "" (L.pack "aName") `shouldParse` L.pack "aName"+ it "parses a string of printable characters"+ $ parse name "" (L.pack "aName")+ `shouldParse` L.pack "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 "" (L.pack "aName bla")+ `shouldParse` L.pack "aName" - it "allows empty names" $- parse name "" (L.pack "") `shouldParse` L.pack ""+ it "allows empty names" $ parse name "" (L.pack "") `shouldParse` L.pack "" describe "node" $ do- it "parses a tree node" $- parse node "" sampleLabelByteString `shouldParse` sampleLabel+ 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 "" (L.pack "")+ `shouldParse` PhyloLabel (L.pack "") 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` Node sampleLabel [] - describe "forest" $- it "parses a set of trees within brackets" $- parse forest "" sampleForestByteString `shouldParse` sampleForest+ describe "forest"+ $ it "parses a set of trees within brackets"+ $ parse forest "" sampleForestByteString+ `shouldParse` sampleForest describe "newick" $ do it "parses a newick trees" $ do- parse newick "" `shouldSucceedOn` sampleNewickByteString1- parse newick "" `shouldSucceedOn` sampleNewickByteString2+ parse (newick Standard) "" `shouldSucceedOn` sampleNewickByteString1+ parse (newick Standard) "" `shouldSucceedOn` sampleNewickByteString2 - it "parses a weird newick tree without node labels nor branch lengths" $- parse newick "" sampleNewickEmptyByteString `shouldParse` sampleNewickEmpty+ 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