diff --git a/ChangeLog.md b/ChangeLog.md
--- a/ChangeLog.md
+++ b/ChangeLog.md
@@ -1,10 +1,13 @@
-# Table of Contents
-
-1. [Changelog for ELynx](#org263dcac)
+# Changelog for ELynx
-
+## 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.
diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -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)
diff --git a/bench/Bench.hs b/bench/Bench.hs
--- a/bench/Bench.hs
+++ b/bench/Bench.hs
@@ -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]]
diff --git a/elynx-tree.cabal b/elynx-tree.cabal
--- a/elynx-tree.cabal
+++ b/elynx-tree.cabal
@@ -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
diff --git a/src/ELynx/Data/Tree.hs b/src/ELynx/Data/Tree.hs
new file mode 100644
--- /dev/null
+++ b/src/ELynx/Data/Tree.hs
@@ -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
diff --git a/src/ELynx/Data/Tree/Bipartition.hs b/src/ELynx/Data/Tree/Bipartition.hs
--- a/src/ELynx/Data/Tree/Bipartition.hs
+++ b/src/ELynx/Data/Tree/Bipartition.hs
@@ -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
diff --git a/src/ELynx/Data/Tree/BranchSupportTree.hs b/src/ELynx/Data/Tree/BranchSupportTree.hs
--- a/src/ELynx/Data/Tree/BranchSupportTree.hs
+++ b/src/ELynx/Data/Tree/BranchSupportTree.hs
@@ -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
diff --git a/src/ELynx/Data/Tree/Distance.hs b/src/ELynx/Data/Tree/Distance.hs
--- a/src/ELynx/Data/Tree/Distance.hs
+++ b/src/ELynx/Data/Tree/Distance.hs
@@ -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) ]
diff --git a/src/ELynx/Data/Tree/EvoTree.hs b/src/ELynx/Data/Tree/EvoTree.hs
deleted file mode 100644
--- a/src/ELynx/Data/Tree/EvoTree.hs
+++ /dev/null
@@ -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
diff --git a/src/ELynx/Data/Tree/MeasurableTree.hs b/src/ELynx/Data/Tree/MeasurableTree.hs
--- a/src/ELynx/Data/Tree/MeasurableTree.hs
+++ b/src/ELynx/Data/Tree/MeasurableTree.hs
@@ -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.
diff --git a/src/ELynx/Data/Tree/Multipartition.hs b/src/ELynx/Data/Tree/Multipartition.hs
--- a/src/ELynx/Data/Tree/Multipartition.hs
+++ b/src/ELynx/Data/Tree/Multipartition.hs
@@ -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
diff --git a/src/ELynx/Data/Tree/NamedTree.hs b/src/ELynx/Data/Tree/NamedTree.hs
--- a/src/ELynx/Data/Tree/NamedTree.hs
+++ b/src/ELynx/Data/Tree/NamedTree.hs
@@ -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
diff --git a/src/ELynx/Data/Tree/PhyloTree.hs b/src/ELynx/Data/Tree/PhyloTree.hs
--- a/src/ELynx/Data/Tree/PhyloTree.hs
+++ b/src/ELynx/Data/Tree/PhyloTree.hs
@@ -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
diff --git a/src/ELynx/Data/Tree/SubSample.hs b/src/ELynx/Data/Tree/SubSample.hs
new file mode 100644
--- /dev/null
+++ b/src/ELynx/Data/Tree/SubSample.hs
@@ -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
+
diff --git a/src/ELynx/Data/Tree/Subset.hs b/src/ELynx/Data/Tree/Subset.hs
deleted file mode 100644
--- a/src/ELynx/Data/Tree/Subset.hs
+++ /dev/null
@@ -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
diff --git a/src/ELynx/Data/Tree/SumStat.hs b/src/ELynx/Data/Tree/SumStat.hs
--- a/src/ELynx/Data/Tree/SumStat.hs
+++ b/src/ELynx/Data/Tree/SumStat.hs
@@ -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
diff --git a/src/ELynx/Data/Tree/Tree.hs b/src/ELynx/Data/Tree/Tree.hs
--- a/src/ELynx/Data/Tree/Tree.hs
+++ b/src/ELynx/Data/Tree/Tree.hs
@@ -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)
diff --git a/src/ELynx/Distribution/BirthDeath.hs b/src/ELynx/Distribution/BirthDeath.hs
--- a/src/ELynx/Distribution/BirthDeath.hs
+++ b/src/ELynx/Distribution/BirthDeath.hs
@@ -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
diff --git a/src/ELynx/Distribution/BirthDeathCritical.hs b/src/ELynx/Distribution/BirthDeathCritical.hs
--- a/src/ELynx/Distribution/BirthDeathCritical.hs
+++ b/src/ELynx/Distribution/BirthDeathCritical.hs
@@ -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
diff --git a/src/ELynx/Distribution/BirthDeathCriticalNoTime.hs b/src/ELynx/Distribution/BirthDeathCriticalNoTime.hs
--- a/src/ELynx/Distribution/BirthDeathCriticalNoTime.hs
+++ b/src/ELynx/Distribution/BirthDeathCriticalNoTime.hs
@@ -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
diff --git a/src/ELynx/Distribution/BirthDeathNearlyCritical.hs b/src/ELynx/Distribution/BirthDeathNearlyCritical.hs
--- a/src/ELynx/Distribution/BirthDeathNearlyCritical.hs
+++ b/src/ELynx/Distribution/BirthDeathNearlyCritical.hs
@@ -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
diff --git a/src/ELynx/Distribution/CoalescentContinuous.hs b/src/ELynx/Distribution/CoalescentContinuous.hs
--- a/src/ELynx/Distribution/CoalescentContinuous.hs
+++ b/src/ELynx/Distribution/CoalescentContinuous.hs
@@ -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)
diff --git a/src/ELynx/Distribution/TimeOfOrigin.hs b/src/ELynx/Distribution/TimeOfOrigin.hs
--- a/src/ELynx/Distribution/TimeOfOrigin.hs
+++ b/src/ELynx/Distribution/TimeOfOrigin.hs
@@ -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
diff --git a/src/ELynx/Distribution/TimeOfOriginNearCritical.hs b/src/ELynx/Distribution/TimeOfOriginNearCritical.hs
--- a/src/ELynx/Distribution/TimeOfOriginNearCritical.hs
+++ b/src/ELynx/Distribution/TimeOfOriginNearCritical.hs
@@ -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
diff --git a/src/ELynx/Distribution/Types.hs b/src/ELynx/Distribution/Types.hs
--- a/src/ELynx/Distribution/Types.hs
+++ b/src/ELynx/Distribution/Types.hs
@@ -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
diff --git a/src/ELynx/Export/Tree/Newick.hs b/src/ELynx/Export/Tree/Newick.hs
--- a/src/ELynx/Export/Tree/Newick.hs
+++ b/src/ELynx/Export/Tree/Newick.hs
@@ -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
diff --git a/src/ELynx/Import/Tree/Newick.hs b/src/ELynx/Import/Tree/Newick.hs
--- a/src/ELynx/Import/Tree/Newick.hs
+++ b/src/ELynx/Import/Tree/Newick.hs
@@ -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 ()
diff --git a/src/ELynx/Simulate/Coalescent.hs b/src/ELynx/Simulate/Coalescent.hs
--- a/src/ELynx/Simulate/Coalescent.hs
+++ b/src/ELynx/Simulate/Coalescent.hs
@@ -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
diff --git a/src/ELynx/Simulate/PointProcess.hs b/src/ELynx/Simulate/PointProcess.hs
--- a/src/ELynx/Simulate/PointProcess.hs
+++ b/src/ELynx/Simulate/PointProcess.hs
@@ -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 ( 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
diff --git a/test/ELynx/Data/Tree/BipartitionSpec.hs b/test/ELynx/Data/Tree/BipartitionSpec.hs
--- a/test/ELynx/Data/Tree/BipartitionSpec.hs
+++ b/test/ELynx/Data/Tree/BipartitionSpec.hs
@@ -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
diff --git a/test/ELynx/Data/Tree/DistanceSpec.hs b/test/ELynx/Data/Tree/DistanceSpec.hs
--- a/test/ELynx/Data/Tree/DistanceSpec.hs
+++ b/test/ELynx/Data/Tree/DistanceSpec.hs
@@ -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
+ )
diff --git a/test/ELynx/Data/Tree/MultipartitionSpec.hs b/test/ELynx/Data/Tree/MultipartitionSpec.hs
--- a/test/ELynx/Data/Tree/MultipartitionSpec.hs
+++ b/test/ELynx/Data/Tree/MultipartitionSpec.hs
@@ -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)
diff --git a/test/ELynx/Data/Tree/PhyloTreeArbitraryInstance.hs b/test/ELynx/Data/Tree/PhyloTreeArbitraryInstance.hs
new file mode 100644
--- /dev/null
+++ b/test/ELynx/Data/Tree/PhyloTreeArbitraryInstance.hs
@@ -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))
diff --git a/test/ELynx/Data/Tree/TreeSpec.hs b/test/ELynx/Data/Tree/TreeSpec.hs
--- a/test/ELynx/Data/Tree/TreeSpec.hs
+++ b/test/ELynx/Data/Tree/TreeSpec.hs
@@ -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
diff --git a/test/ELynx/Export/Tree/NewickSpec.hs b/test/ELynx/Export/Tree/NewickSpec.hs
--- a/test/ELynx/Export/Tree/NewickSpec.hs
+++ b/test/ELynx/Export/Tree/NewickSpec.hs
@@ -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
diff --git a/test/ELynx/Import/Tree/NewickSpec.hs b/test/ELynx/Import/Tree/NewickSpec.hs
--- a/test/ELynx/Import/Tree/NewickSpec.hs
+++ b/test/ELynx/Import/Tree/NewickSpec.hs
@@ -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