[med-svn] [r-cran-phangorn] 01/08: New upstream version 2.1.1
Andreas Tille
tille at debian.org
Wed Dec 7 08:46:53 UTC 2016
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository r-cran-phangorn.
commit 6c5685461593a2c0d17e27e9c29fcd75f69a9e6c
Author: Andreas Tille <tille at debian.org>
Date: Wed Dec 7 09:25:46 2016 +0100
New upstream version 2.1.1
---
DESCRIPTION | 22 +-
MD5 | 124 ++--
NAMESPACE | 26 +-
NEWS | 36 +
R/Densi.R | 231 ------
R/RcppExports.R | 11 +
R/ancestral_pml.R | 32 +-
R/bootstrap.R | 55 +-
R/distSeq.R | 97 +--
R/distTree.R | 104 ++-
R/fitch.R | 3 +-
R/networx.R | 120 +++-
R/parsimony.R | 46 +-
R/phyDat.R | 388 ++++++-----
R/phylo.R | 36 +-
R/sankoff.R | 2 -
R/simSeq.R | 17 +-
R/superTree.R | 135 ++++
R/treeManipulation.R | 65 +-
R/treedist.R | 154 +++-
R/zzz.R | 5 +
README.md | 6 +-
build/vignette.rds | Bin 375 -> 376 bytes
inst/doc/Ancestral.pdf | Bin 260891 -> 245485 bytes
inst/doc/IntertwiningTreesAndNetworks.html | 6 +-
inst/doc/Networx.Rmd | 2 +-
inst/doc/Networx.html | 16 +-
inst/doc/Trees.pdf | Bin 167890 -> 167906 bytes
inst/doc/phangorn-specials.pdf | Bin 180115 -> 179955 bytes
man/SOWH.test.Rd | 2 +-
man/ancestral.pml.Rd | 2 +-
man/as.splits.Rd | 9 +-
man/densiTree.Rd | 7 +-
man/dist.p.Rd | 2 +-
man/hadamard.Rd | 14 +-
man/maxCladeCred.Rd | 2 +-
man/midpoint.Rd | 2 +-
man/modelTest.Rd | 4 +-
man/nni.Rd | 2 +-
man/parsimony.Rd | 5 +-
man/phyDat.Rd | 19 +-
man/plot.networx.Rd | 2 +-
man/pml.Rd | 4 +-
man/pml.fit.Rd | 31 +-
man/simSeq.Rd | 5 +-
man/superTree.Rd | 22 +-
man/treedist.Rd | 20 +-
src/RcppExports.cpp | 30 +
src/dist.c | 9 +-
src/fitch.c | 18 +-
src/ml.c | 51 +-
src/phangorn.c | 25 +-
src/phangorn_help.cpp | 48 ++
src/sankoff.c | 11 +-
src/sprdist.c | 1047 ++++++++++++++++++++++++++++
tests/testthat/test_dist_tree.R | 22 +
tests/testthat/test_distances.R | 22 +-
tests/testthat/test_hadamard.R | 21 +
tests/testthat/test_parsimony.R | 34 +-
tests/testthat/test_phyDat.R | 78 ++-
tests/testthat/test_pml.R | 93 ++-
tests/testthat/test_splits.R | 39 ++
tests/testthat/test_superTree.R | 26 +
tests/testthat/test_treeManipulation.R | 31 +-
tests/testthat/test_treedist.R | 27 +
vignettes/Networx.Rmd | 2 +-
vignettes/phangorn.bib | 79 +++
67 files changed, 2681 insertions(+), 925 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 944dc6b..b4eaf19 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,18 +1,23 @@
Package: phangorn
Title: Phylogenetic Analysis in R
-Version: 2.0.4
-Date: 2016-06-20
+Version: 2.1.1
+Date: 2016-12-03
Authors at R: c(person("Klaus", "Schliep", email="klaus.schliep at gmail.com", role = c("aut", "cre")),
person("Emmanuel", "Paradis", role = c("aut")),
+ person("Leonardo", "de Oliveira Martins", role = c("aut")),
person("Alastair", "Potts", role = c("aut")),
+ person("Tim W.", "White", role = c("aut")),
+ person("Cyrill", "Stachniss", role = c("ctb")),
person("Michelle", "Kendall", email="m.kendall at imperial.ac.uk", role = c("ctb")))
Maintainer: Klaus Schliep <klaus.schliep at gmail.com>
Description: Phylogenetic analysis in R: Estimation of phylogenetic
trees and networks using Maximum Likelihood, Maximum Parsimony,
distance methods and Hadamard conjugation.
-Depends: R (>= 3.0.0), ape (>= 3.5)
-Imports: quadprog, igraph (>= 1.0), Matrix, parallel, nnls, methods,
- utils, stats, graphics, grDevices, fastmatch
+Depends: R (>= 3.2.0), ape (>= 4.0)
+Imports: quadprog, igraph (>= 1.0), Matrix, parallel, methods, utils,
+ stats, graphics, grDevices, fastmatch, magrittr, Rcpp (>=
+ 0.12.0)
+LinkingTo: Rcpp
Suggests: testthat, seqLogo, seqinr, xtable, flashClust, rgl, knitr,
rmarkdown, Biostrings
ByteCompile: TRUE
@@ -21,9 +26,12 @@ VignetteBuilder: utils, knitr
URL: https://github.com/KlausVigo/phangorn
Repository: CRAN
NeedsCompilation: yes
-Packaged: 2016-06-20 15:44:07 UTC; klaus
+Packaged: 2016-12-03 22:25:47 UTC; klaus
Author: Klaus Schliep [aut, cre],
Emmanuel Paradis [aut],
+ Leonardo de Oliveira Martins [aut],
Alastair Potts [aut],
+ Tim W. White [aut],
+ Cyrill Stachniss [ctb],
Michelle Kendall [ctb]
-Date/Publication: 2016-06-21 08:32:40
+Date/Publication: 2016-12-04 11:42:04
diff --git a/MD5 b/MD5
index 9156cd4..804bd2e 100644
--- a/MD5
+++ b/MD5
@@ -1,33 +1,35 @@
-1b369c1893556a9402077452b8e55651 *DESCRIPTION
-2552be19064d4d91f8600cfc228052bd *NAMESPACE
-d5f86a348bb7b8b442738f553682d61d *NEWS
+1caf854c8e3d0058f87b4e52606bf9a6 *DESCRIPTION
+f8dbf6fef0c3d19291595e95e407f2eb *NAMESPACE
+bbe4b979bbcd18715eee2072e9305ed8 *NEWS
9441c4e1029cac2a635d71e9acad4ffb *R/Coalescent.R
-63bb6ba4ff0d5fa9c07ab8f7e014faca *R/Densi.R
+2abb1811fe55334d88817d7c7ad2984d *R/Densi.R
+b743c4769d78295ac28702f376be88f7 *R/RcppExports.R
0708ebbd9a9a0ba33f932d172d01ecc6 *R/SH.R
f916be823b9a838d49145695b0a37aa1 *R/SOWH.R
-8b8f9aeefc014d7ff3d0cd51077f4650 *R/ancestral_pml.R
-f0570fdfa97da0f51b063b88ce85057d *R/bootstrap.R
+3cfc639a09a0e3c003929e7ff015df46 *R/ancestral_pml.R
+fd1201d71dddabc5bedf25753bf0d317 *R/bootstrap.R
089d06695316d8734e5647d50eaa2b12 *R/cladePar.R
b308a91826087b5dfcbc2075a0cc3a90 *R/clanistic.R
62d2375bf3497590ad11af0de6e77efa *R/dist.p.R
-857b3f1c2ad9d97648c81b20105b4d2a *R/distSeq.R
-fbf4c9f43a5cb52edf3779f125115cd9 *R/distTree.R
-cff63aa886a22f913ad0c01b53ff9fb8 *R/fitch.R
+64b92983be4d54d166713ba0e69cd6a9 *R/distSeq.R
+f13145e51b142a4124e83ecdf947a9d2 *R/distTree.R
+31c6e8e64ddaaef935571c2b0a136d69 *R/fitch.R
3c0ccfb95517188a1d8f4059b7f192d2 *R/hadamard.R
42fe926e668d07ef6aee042b4a69de95 *R/modelTest.R
c4c58c9e35b58e21176f482fd6150002 *R/neighborNet.R
-2426acbef78a4ea50e95c20df9b64b5e *R/networx.R
-f411f2c21b87776125e0d3c510e17b26 *R/parsimony.R
-fc739a28612271d61f81d221e01a0005 *R/phyDat.R
-974b4bdd5db6a33a750f6bce183f47ed *R/phylo.R
-4832725b6e807c67398b7c25621161d8 *R/sankoff.R
-eec6086aac90ebcb9403402682295c9b *R/simSeq.R
+76da6b713d332fcce747148396563a1b *R/networx.R
+3223d2ae933dadce81bc77821cc99cde *R/parsimony.R
+c90c106b1ea0dc8dc6f088b3fd51b280 *R/phyDat.R
+7faa9011d21ceb145671ad0f32a54d00 *R/phylo.R
+b33943cef8181fd573939c9e1a003fdb *R/sankoff.R
+1e46863e6ff14de9597b86c4dedbd02d *R/simSeq.R
+9f64c372e0892738ccb107683440af5d *R/superTree.R
2b4d6884b4ee2b92e617cea76f79d3da *R/sysdata.rda
-979aea186a4b879b66420ed49a25285d *R/treeManipulation.R
-39801dcd237ea83998e2ed1e82fb59b8 *R/treedist.R
-c925e4cd5ce95bf91968fc909d6eef76 *R/zzz.R
-058f733944f5c56b40167cb2f3dd29bc *README.md
-2222bf60c7d441a496a74fe403b025f4 *build/vignette.rds
+40a44fc429cb60a112307c22f8f0de8f *R/treeManipulation.R
+fcffc8940c7bee194489d3b8180bced0 *R/treedist.R
+e07b91ecaa99e7e6023fed5d9760705e *R/zzz.R
+7191ca543e47e3eaf4ec6fb397d2e079 *README.md
+4a4c3e4965d40f619e29c1707cb472f0 *build/vignette.rds
4a92b07bcd170e85bd237e460acafa93 *data/Laurasiatherian.RData
4b269c2e640293341b9b1c04d4dd7f4e *data/chloroplast.RData
19f84425e8caaf9f605490cdff17dd81 *data/yeast.RData
@@ -35,19 +37,19 @@ b4d3fa3f40aae3a68846d1f28274e9a0 *inst/CITATION
61b3722df2ead748e0db8595428412a1 *inst/README
56efa08410f500fc3b4851b71058ad27 *inst/doc/Ancestral.R
bf3e208e4fceec83ea3d522adb1b0272 *inst/doc/Ancestral.Rnw
-74d40684ccc50a3970cdfb50d2ad72d8 *inst/doc/Ancestral.pdf
+d3291689993a930c8352226395dfedf1 *inst/doc/Ancestral.pdf
0862fe6de86381ba4ad9c8d1b4b02198 *inst/doc/IntertwiningTreesAndNetworks.R
af3eb8b8f4e4512cb5e22c8cec8538eb *inst/doc/IntertwiningTreesAndNetworks.Rmd
-7a58319240ce93cc06c335f3237b7a0d *inst/doc/IntertwiningTreesAndNetworks.html
+c274e3f8a6d69e55ac9b75c2972d140d *inst/doc/IntertwiningTreesAndNetworks.html
d24a83a2803108ede87eff4de3c13ae8 *inst/doc/Networx.R
-4d2150a5608862c1eff7000e4096e44c *inst/doc/Networx.Rmd
-9f24fa87e64cd7ef1056fa34b0d5d0d7 *inst/doc/Networx.html
+fe641c0b4c7aac7aabf24c47ed6cde40 *inst/doc/Networx.Rmd
+640f7d533c3af0b9b62e6762e33cec3c *inst/doc/Networx.html
1effe50d23f90af16c248d52ddedda58 *inst/doc/Trees.R
5aa283308c069783e018a708a19c6439 *inst/doc/Trees.Rnw
-6119d31f15debb53dda3436a6acbdcb6 *inst/doc/Trees.pdf
+8b7350a9103d1cf8226c306c70a0601e *inst/doc/Trees.pdf
ef0fbc15af66dff9fb288de1abbf8901 *inst/doc/phangorn-specials.R
25becba1a0c2b64b7166f0b5eb72c4f2 *inst/doc/phangorn-specials.Rnw
-8aab7335a82edfd290c74b4a6c911efc *inst/doc/phangorn-specials.pdf
+3ae1044de42534384f6d31000b699d1f *inst/doc/phangorn-specials.pdf
3009f9da02a198cb557d435cc5ad8c7f *inst/extdata/Blosum62.dat
72f496d4c6e937ffe744d25bd8891313 *inst/extdata/Dayhoff.dat
5aa357dab0e72b27023e4302bc37dbad *inst/extdata/FLU.dat
@@ -76,11 +78,11 @@ e185a8c9c7aed17a3e216cb78d517fd4 *inst/extdata/trees/woodmouse.mrbayes.nex.run2.
bda5669b71de975e1309d909c495b71f *man/Laurasiatherian.Rd
4d6fb151c28a597eadc0d4ec2c59a504 *man/NJ.Rd
5dd84777963b553e8c89bddc951b17c8 *man/SH.test.Rd
-b3418f878911164dd974046f0492c79c *man/SOWH.test.Rd
+095dc6d3e4c11d8ae6a83967a88dfb4c *man/SOWH.test.Rd
0ad5edc3e25f1b6bb10c816e8b866320 *man/addConfidences.Rd
7e03e8af0658fc2b19d1e51b0fad7a26 *man/allTrees.Rd
-c7080f03a162f87e0c700028bd699516 *man/ancestral.pml.Rd
-9f8c3cdeffbbbdc184680e0b797cfca2 *man/as.splits.Rd
+8a9aaa3662de9d68b72a0790cb255866 *man/ancestral.pml.Rd
+074100ce52115557988ff5f316420999 *man/as.splits.Rd
3b53888e5f37fd54344286a9175d8dc2 *man/bab.Rd
7c2aead89db45adc53a75fbecb93e90e *man/bootstrap.pml.Rd
22ad0ef60c2edd8b3d003f170f2fa15a *man/chloroplast.Rd
@@ -88,58 +90,64 @@ c7080f03a162f87e0c700028bd699516 *man/ancestral.pml.Rd
5befa0cf9c4edc213566b3f60f8838ed *man/consensusNet.Rd
55baf4b98181f03c2c1489a98d504e2d *man/cophenetic.networx.Rd
c090d7ec936738daaea76e99f1f9bab0 *man/delta.score.Rd
-5868068eb5f8be4302c0c3caadd25943 *man/densiTree.Rd
+50af51e1ddb4295a8ee2791e49237f00 *man/densiTree.Rd
f44e43abbbf7707821627c8e88f6e9fa *man/designTree.Rd
b888d59a2b1fe266b098544b8ac86eed *man/dfactorial.Rd
1e97d5bbdae42152778a84aff952a8cd *man/dist.hamming.Rd
-c240de2acabd4c5ccec785a09d21608c *man/dist.p.Rd
+5505391fb32d03f9a3c2b8a14ebf03de *man/dist.p.Rd
789d097683def6234f407e69cc3ee9b2 *man/distanceHadamard.Rd
72abceeb247f22b3da6560df0c20468a *man/getClans.Rd
-c1cacff95f20b574c03126d3f8e24106 *man/hadamard.Rd
+619845889ee2d3b241af8a114eac9330 *man/hadamard.Rd
f8cf4238c40481ee6c09cc935fc07f12 *man/lento.Rd
-734851dec63f4aab42447cd0e31b4641 *man/maxCladeCred.Rd
-a4d9197a0880a1d6b12d2e9371d5cd53 *man/midpoint.Rd
-a7d146d4074b23c6ceb162f602f2a669 *man/modelTest.Rd
+aec58935c6251ba62808ea67d130a7f2 *man/maxCladeCred.Rd
+e18f9490c68933d8f946f3d3d3993f25 *man/midpoint.Rd
+30c8fe2f9863788d48b2846d9a303905 *man/modelTest.Rd
7c85721291853f713eac8637c358f7a5 *man/neighborNet.Rd
-4dcf6f133fbbc6c0e3e43413d33ea3f7 *man/nni.Rd
-974c5df5eec7ea710deffd027a37e43e *man/parsimony.Rd
+14f5a70080bc79967555dfcbe839c459 *man/nni.Rd
+ad48146c276212ea82c236bbdabc923f *man/parsimony.Rd
ad095b45d2e89f45bc9938857c0bbe6d *man/phangorn-package.Rd
-40aabd33c9dc2db426fd36fe2cc61cf5 *man/phyDat.Rd
-fa86e2334aabf74dcdc7fd15ce78bbb4 *man/plot.networx.Rd
-ddc71ef3ac07e482ef2c97bff7cb1e97 *man/pml.Rd
-7aeb801f22ae8c8f712052e432ebf564 *man/pml.fit.Rd
+cf9e0177d4e92a1576d5f5e7d158cb4c *man/phyDat.Rd
+8be8be7777fcec2c8ede18bead01240c *man/plot.networx.Rd
+9d9a2c930b17961fef4845dd6c52991e *man/pml.Rd
+f30af37ac4943192c369d32a712751ba *man/pml.fit.Rd
9506d15a81467d18f40bb4741b3d9d28 *man/pmlCluster.Rd
8b06bdee57c510a690ea04b40bac4844 *man/pmlMix.Rd
ff67d17fa29bf88cea645905334c5ecc *man/pmlPart.Rd
bfc83a71d7de4e1c0b994b8b9c853439 *man/read.aa.Rd
-ae8bb6597edcdb96eeb86b5367cc4d5f *man/simSeq.Rd
+6d5fd22b2c6dbe501b0ea18f7b96cf85 *man/simSeq.Rd
2e7848b9bac018bde56302f70066a49a *man/splitsNetwork.Rd
-2abee82ad55769e4bf99f0a2f0dc102e *man/superTree.Rd
-345fb2ad0be13f59ee5de23ccc3b60ae *man/treedist.Rd
+1276a142c217c44f2e0b1277a600c25d *man/superTree.Rd
+ba64cab5b9f597473c7f465a9f4eca34 *man/treedist.Rd
64399c9fb7fb25d103c25aa925ab7a10 *man/upgma.Rd
b97649fe8b91a0632a1ded89a6f43125 *man/yeast.Rd
9a8d913f59a18a676121ea930cbe958f *src/Makevars
-c73230a47eefd569a5a63e6b2e43b736 *src/dist.c
-fc069fca9f2d8216cd6f7f458af13721 *src/fitch.c
-36ec67f90b9d4d79ddca24440c56b099 *src/ml.c
-e79d2a91b347507dccb722295db682ef *src/phangorn.c
+3178e5cc701bd341b1f3fa3626f1e0b5 *src/RcppExports.cpp
+80c91abd1d67687878dc151c7440aaf9 *src/dist.c
+4244d0890739f29f204ad5745970f151 *src/fitch.c
+b36357bbe2bc513c2216e905c30f173b *src/ml.c
+1c1ea4701d57f7b819a659fdc586aa63 *src/phangorn.c
+13be023855d3786e428a079593ed4177 *src/phangorn_help.cpp
d8e5c864aa0b122ce426a67479d3a610 *src/read_aa.c
-ebbc6aa310a27020920560e78f71ad72 *src/sankoff.c
+636bf81aacbb5ac0e06d0ae063a0ce75 *src/sankoff.c
+2910bd3e2e2383f7bac46180d62243a4 *src/sprdist.c
7d1eaef831fb2372b367f5be1b0a5790 *tests/testthat.R
-85915294ad9354b55ae10dd2f00ea595 *tests/testthat/test_distances.R
-ee0d53cf5ddd4840ec6649eef91b2df7 *tests/testthat/test_parsimony.R
-3b6b3635995a131b5beb515ae2897e28 *tests/testthat/test_phyDat.R
-9e246e3c61392cb3d4fcd76669e237e3 *tests/testthat/test_pml.R
-feb64ac20ca2b3749fe89fd5cb4d3863 *tests/testthat/test_splits.R
-9d0145ac2a13bcb4f1f4099e25d8cb45 *tests/testthat/test_treeManipulation.R
-b6a732e1ab1ab285f66aca5a22720039 *tests/testthat/test_treedist.R
+97d94788752f552d019cdaafcf7a66bc *tests/testthat/test_dist_tree.R
+92512e2f264d4a11382ddd89e2a515d1 *tests/testthat/test_distances.R
+7f6d10a2be16a7c04c1591ed79d7e02e *tests/testthat/test_hadamard.R
+59cb415b0af45af3bd1494226260791e *tests/testthat/test_parsimony.R
+06a0631572d94f45d220b55847f9f661 *tests/testthat/test_phyDat.R
+870574752144d6ca09f3c3a784abe284 *tests/testthat/test_pml.R
+7688c6569491421417861398008a0232 *tests/testthat/test_splits.R
+20ec41bfea311dd762e984f834a26948 *tests/testthat/test_superTree.R
+b46633e27e0dd867a89d4c40e4eabfb3 *tests/testthat/test_treeManipulation.R
+98d3d4091f9bacf2771d037df005daaf *tests/testthat/test_treedist.R
bf3e208e4fceec83ea3d522adb1b0272 *vignettes/Ancestral.Rnw
af3eb8b8f4e4512cb5e22c8cec8538eb *vignettes/IntertwiningTreesAndNetworks.Rmd
-4d2150a5608862c1eff7000e4096e44c *vignettes/Networx.Rmd
+fe641c0b4c7aac7aabf24c47ed6cde40 *vignettes/Networx.Rmd
a8250afd3b341e6c3f889e9454d5bc5f *vignettes/Trees.RData
5aa283308c069783e018a708a19c6439 *vignettes/Trees.Rnw
6af5f4d4c2e93469cc19d46a97ab5d0f *vignettes/exdna.txt
1ca8b1d97a011a9d0ecd9630d548dfb3 *vignettes/movie.gif
25becba1a0c2b64b7166f0b5eb72c4f2 *vignettes/phangorn-specials.Rnw
-d8f77ef2f47f4ca67f5e3e85ae952f12 *vignettes/phangorn.bib
+16fa639fbbf08ae193a36f714197a2fa *vignettes/phangorn.bib
d3069d1eff9e70bed655b8962bf4ee2b *vignettes/primates.dna
diff --git a/NAMESPACE b/NAMESPACE
index 76a1ded..27af5dd 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -4,21 +4,22 @@ importFrom(ape, as.DNAbin, as.phylo, plot.phylo, as.prop.part, nj,old2new.phylo,
new2old.phylo, cophenetic.phylo, is.rooted, unroot, root, is.binary.tree, as.alignment,
di2multi, multi2di, .uncompressTipLabel, .compressTipLabel, prop.part, Ntip,
postprocess.prop.part, speciesTree, plotPhyloCoor, cladogram.plot, phylogram.plot,
- node.depth.edgelength, drop.tip, stree, rtree, is.ultrametric, .PlotPhyloEnv,
+ node.depth.edgelength, drop.tip, stree, rtree, rcoal, rmtree, is.ultrametric, .PlotPhyloEnv,
nodelabels, edgelabels, BOTHlabels, read.nexus.data, write.nexus.data,
read.dna, write.dna, reorder.phylo, dist.nodes, collapse.singles, rotate, prop.clades,
- fastme.bal, fastme.ols, Nnode, image.DNAbin, ONEwise)
-# as.AAbin, as.phyDat.AAbin, image.AAbin
+ fastme.bal, fastme.ols, Nnode, image.DNAbin, ONEwise, as.AAbin, as.phyDat.AAbin, image.AAbin)
+# is.binary,
importFrom(igraph, graph, graph_from_adjacency_matrix, layout_with_kk, plot.igraph,
shortest_paths, set.edge.attribute, set_edge_attr, decompose, components, groups,
- vcount, graph.edgelist, graph_from_edgelist, topo_sort)
+ vcount, graph.edgelist, graph_from_edgelist, topo_sort, layout_nicely)
#layout.kamada.kawai, graph.adjacency, get.shortest.paths, topological.sort,
-
+importFrom(Rcpp,sourceCpp)
+importFrom(magrittr,"%>%")
importFrom(Matrix, Matrix, sparseMatrix, spMatrix, crossprod, solve) #, %*%
importFrom(fastmatch, fmatch)
-importFrom(nnls, nnls)
+#importFrom(nnls, nnls)
#importFrom(parallel, mclapply, detectCores)
import(parallel)
importFrom(quadprog, solve.QP, solve.QP.compact)
@@ -40,22 +41,23 @@ importFrom(methods, hasArg)
export(pml, optim.pml, pml.control, parsimony, optim.parsimony, pratchet, NJ, UNJ, PNJ,
phyDat, cbind.phyDat, read.phyDat, write.phyDat, as.Matrix, as.phyDat, as.splits, as.networx,
- dfactorial, ldfactorial, hadamard, nni, allSitePattern, allSplits, fhm,
+ dfactorial, ldfactorial, hadamard, nni, allSitePattern, allSplits, allCircularSplits, fhm,
distanceHadamard, treedist, sankoff, fitch, h2st, h4st, dist.logDet, dist.hamming,
dist.ml, dist.p, upgma, wpgma, write.nexus.splits, read.nexus.splits,
read.nexus.networx, write.nexus.networx, write.splits, pmlPart,
pmlCluster, pmlMix, pmlPen, read.aa, allTrees, designTree, designSplits, nnls.tree,
nnls.splits, nnls.phylo, nnls.networx, neighborNet, pmlPart2multiPhylo,
splitsNetwork, simSeq, SH.test, bootstrap.pml, bootstrap.phyDat, RF.dist, wRF.dist,
- KF.dist, path.dist, rNNI,
+ KF.dist, path.dist, sprdist, SPR.dist, rNNI,
rSPR, plotBS, Ancestors, Descendants, mrca.phylo, Children, Siblings, pace, modelTest,
lento, compatible, acgt2ry, ancestral.pars, ancestral.pml, CI, RI, getClans, getSlices,
getClips, getDiversity, midpoint, pruneTree, acctran, getRoot, plotAnc, consensusNet,
bab, random.addition, diversity, baseFreq, densiTree, superTree, coalSpeciesTree,
pml.fit, pml.init, pml.free, edQt, lli, cladePar, addConfidences, createLabel,
countCycles, SOWH.test, plot.networx, presenceAbsence, as.networx.splits, addTrivialSplits,
- phyDat2alignment, AICc, readDist, writeDist, discrete.gamma, matchSplits,
- removeUndeterminedSites, phyDat2MultipleAlignment, delta.score, maxCladeCred)
+ phyDat2alignment, AICc, readDist, writeDist, discrete.gamma, matchSplits,
+ distinct.splits, c.splits, dna2codon, codon2dna,
+ removeUndeterminedSites, phyDat2MultipleAlignment, delta.score, maxCladeCred, as.MultipleAlignment)
S3method('[', splits)
S3method(addConfidences, phylo)
@@ -77,11 +79,13 @@ S3method(as.networx, phylo)
#S3method(as.igraph, networx)
S3method(as.phyDat, DNAbin)
S3method(as.phyDat, alignment)
+S3method(as.phyDat, character)
S3method(as.phyDat, data.frame)
+S3method(as.phyDat, factor)
S3method(as.phyDat, matrix)
S3method(as.phyDat, MultipleAlignment)
S3method(as.MultipleAlignment, phyDat)
-#S3method(as.AAbin, phyDat) # momentan raus
+S3method(as.AAbin, phyDat)
S3method(as.DNAbin, phyDat)
S3method(as.phylo, splits)
S3method(as.splits, phylo)
diff --git a/NEWS b/NEWS
index aa4130a..92cee36 100644
--- a/NEWS
+++ b/NEWS
@@ -1,3 +1,39 @@
+ CHANGES in PHANGORN VERSION 2.1.0
+
+NEW FEATURES
+
+ o new functions to compute the (approximate) SPR distance (sprdist, SPR.dist)
+
+ contributed by Leonardo de Oliveira Martins.
+
+ o super tree methods based on NNI and SPR distances
+
+BUG FIXES
+
+ o fixed bug in KF.dist
+
+OTHER CHANGES
+
+ o improvements to as.networx. It often now produces networks with less edges
+
+ resulting in much nicer plots
+
+ o plot.networx does take a different layout algorithm
+
+ o as.data.frame.phyDat and as.character.phyDat return amino acids
+
+ now in upper cases
+
+ o more unit tests
+
+ o improved cbind.phyDat, faster and more flexible
+
+ o phangorn now requires ape 4.0
+
+ o phangorn now imports Rcpp, but not nnls any more
+
+
+
CHANGES in PHANGORN VERSION 2.0.4
NEW FEATURES
diff --git a/R/Densi.R b/R/Densi.R
index aa71547..db53164 100644
--- a/R/Densi.R
+++ b/R/Densi.R
@@ -18,237 +18,6 @@ getAges <- function(x){
}
-# now more memoryefficient
-# from phytools code by Liam Revell with a few changes
-my.supertree<-function(trees,method=c("pratchet","optim.parsimony"), trace=0, ...){
- # set method
- method<-method[1]
- # some minor error checking
- if(!inherits(trees,"multiPhylo")) stop("trees must be object of class 'multiPhylo.'")
-
- labels <- lapply(trees, function(x)sort(x$tip.label))
- ulabels <- unique(labels)
- lul <- length(ulabels)
- # compute matrix representation phylogenies
- X<-vector("list", lul) # list of bipartitions
- characters<-0 # number of characters
- weights <- NULL
- species<-trees[[1]]$tip.label
- for(i in 1:lul){
- pos <- match(labels, ulabels[i])
- ind <- which(!is.na(pos))
- temp<-prop.part(trees[ind]) # find all bipartitions
- # create matrix representation of trees[[i]] in X[[i]]
- X[[i]]<-matrix(0,nrow=length(trees[[ind[1]]]$tip.label),ncol=length(temp)-1)
- for(j in 1:ncol(X[[i]])) X[[i]][c(temp[[j+1]]),j]<-1
- rownames(X[[i]])<-attr(temp,"labels") # label rows
-# if(i==1) species<-trees[[ind[1]]]$tip.label
-# else
- species<-union(species,trees[[ind[1]]]$tip.label) # accumulate labels
- characters<-characters+ncol(X[[i]]) # count characters
- weights <- c(weights, attr(temp, "number")[-1])
- }
- XX<-matrix(data="?",nrow=length(species),ncol=characters,dimnames=list(species))
- j<-1
- for(i in 1:length(X)){
- # copy each of X into supermatrix XX
- XX[rownames(X[[i]]),c(j:((j-1)+ncol(X[[i]])))]<-X[[i]][1:nrow(X[[i]]),1:ncol(X[[i]])]
- j<-j+ncol(X[[i]])
- }
- # compute contrast matrix for phangorn
- contrast<-matrix(data=c(1,0,0,1,1,1),3,2,dimnames=list(c("0","1","?"),c("0","1")),byrow=TRUE)
- # convert XX to phyDat object
- XX<-phyDat(XX,type="USER",contrast=contrast, compress=FALSE)
- attr(XX, "weight") <- weights
- # estimate supertree
- if(method=="pratchet"){
- if(hasArg(start)){
- start<-list(...)$start
- if(inherits(start,"phylo")){
- supertree<-pratchet(XX,all=TRUE, trace=0, ...)
- } else {
- if(start=="NJ") start<-NJ(dist.hamming(XX))
- else if(start=="random") start<-rtree(n=length(XX),tip.label=names(XX))
- else {
- warning("do not recognize that option for start; using random starting tree")
- tree<-rtree(n=length(XX),tip.label=names(XX))
- }
- args<-list(...)
- args$start<-start
- args$data<-XX
- args$all<-TRUE
- supertree<-do.call(pratchet,args)
- }
- } else supertree<-pratchet(XX,all=TRUE, trace=0, ...)
- if(inherits(supertree,"phylo"))
- if(trace>0)message(paste("The MRP supertree, optimized via pratchet(),\nhas a parsimony score of ",
- attr(supertree,"pscore")," (minimum ",characters,")",sep=""))
- else if(inherits(supertree,"multiPhylo"))
- if(trace>0)message(paste("pratchet() found ",length(supertree)," supertrees\nwith a parsimony score of ",
- attr(supertree[[1]],"pscore")," (minimum ",characters,")",sep=""))
- } else if(method=="optim.parsimony"){
- if(hasArg(start)){
- start<-list(...)$start
- if(inherits(start,"phylo")){
- supertree<-optim.parsimony(tree=start,data=XX, trace=0, ...)
- } else {
- if(start=="NJ") start<-NJ(dist.hamming(XX))
- else if(start=="random") start<-rtree(n=length(XX),tip.label=names(XX))
- else {
- warning("do not recognize that option for tree; using random starting tree")
- start<-rtree(n=length(XX),tip.label=names(XX))
- }
- supertree<-optim.parsimony(tree=start,data=XX,...)
- }
- } else {
- if(trace>0)message("no input starting tree or option for optim.parsimony; using random addition tree")
- start<-random.addition(XX) # rtree(n=length(XX),tip.label=names(XX))
- supertree<-optim.parsimony(tree=start,data=XX, trace=0, ...)
- }
- if(inherits(supertree,"phylo"))
- if(trace>0)message(paste("The MRP supertree, optimized via optim.parsimony(),\nhas a parsimony score of ",
- attr(supertree,"pscore")," (minimum ",characters,")",sep=""))
- else if(inherits(supertree,"multiPhylo"))
- if(trace>0)message(paste("optim.parsimony() found ",length(supertree)," supertrees\nwith a parsimony score of ",
- attr(supertree[[1]],"pscore")," (minimum ",characters,")",sep=""))
- }
- return(supertree)
-}
-
-
-my.supertree.Old<-function(trees,method=c("pratchet","optim.parsimony"), trace=0, ...){
- # set method
- method<-method[1]
- # some minor error checking
- if(!inherits(trees,"multiPhylo")) stop("trees must be object of class 'multiPhylo.'")
- # compute matrix representation phylogenies
- X<-list() # list of bipartitions
- characters<-0 # number of characters
- for(i in 1:length(trees)){
- temp<-prop.part(trees[[i]]) # find all bipartitions
- # create matrix representation of trees[[i]] in X[[i]]
- X[[i]]<-matrix(0,nrow=length(trees[[i]]$tip.label),ncol=length(temp)-1)
- for(j in 1:ncol(X[[i]])) X[[i]][c(temp[[j+1]]),j]<-1
- rownames(X[[i]])<-attr(temp,"labels") # label rows
- if(i==1) species<-trees[[i]]$tip.label
- else species<-union(species,trees[[i]]$tip.label) # accumulate labels
- characters<-characters+ncol(X[[i]]) # count characters
- }
- XX<-matrix(data="?",nrow=length(species),ncol=characters,dimnames=list(species))
- j<-1
- for(i in 1:length(X)){
- # copy each of X into supermatrix XX
- XX[rownames(X[[i]]),c(j:((j-1)+ncol(X[[i]])))]<-X[[i]][1:nrow(X[[i]]),1:ncol(X[[i]])]
- j<-j+ncol(X[[i]])
- }
- # compute contrast matrix for phangorn
- contrast<-matrix(data=c(1,0,0,1,1,1),3,2,dimnames=list(c("0","1","?"),c("0","1")),byrow=TRUE)
- # convert XX to phyDat object
- XX<-phyDat(XX,type="USER",contrast=contrast)
- # estimate supertree
- if(method=="pratchet"){
- if(hasArg(start)){
- start<-list(...)$start
- if(inherits(start,"phylo")){
- supertree<-pratchet(XX,all=TRUE, trace=0, ...)
- } else {
- if(start=="NJ") start<-NJ(dist.hamming(XX))
- else if(start=="random") start<-rtree(n=length(XX),tip.label=names(XX))
- else {
- warning("do not recognize that option for start; using random starting tree")
- tree<-rtree(n=length(XX),tip.label=names(XX))
- }
- args<-list(...)
- args$start<-start
- args$data<-XX
- args$all<-TRUE
- supertree<-do.call(pratchet,args)
- }
- } else supertree<-pratchet(XX,all=TRUE, trace=0, ...)
- if(inherits(supertree,"phylo"))
- if(trace>0)message(paste("The MRP supertree, optimized via pratchet(),\nhas a parsimony score of ",
- attr(supertree,"pscore")," (minimum ",characters,")",sep=""))
- else if(inherits(supertree,"multiPhylo"))
- if(trace>0)message(paste("pratchet() found ",length(supertree)," supertrees\nwith a parsimony score of ",
- attr(supertree[[1]],"pscore")," (minimum ",characters,")",sep=""))
- } else if(method=="optim.parsimony"){
- if(hasArg(start)){
- start<-list(...)$start
- if(inherits(start,"phylo")){
- supertree<-optim.parsimony(tree=start,data=XX, trace=0, ...)
- } else {
- if(start=="NJ") start<-NJ(dist.hamming(XX))
- else if(start=="random") start<-rtree(n=length(XX),tip.label=names(XX))
- else {
- warning("do not recognize that option for tree; using random starting tree")
- start<-rtree(n=length(XX),tip.label=names(XX))
- }
- supertree<-optim.parsimony(tree=start,data=XX,...)
- }
- } else {
- if(trace>0)message("no input starting tree or option for optim.parsimony; using random addition tree")
- start<-random.addition(XX) # rtree(n=length(XX),tip.label=names(XX))
- supertree<-optim.parsimony(tree=start,data=XX, trace=0, ...)
- }
- if(inherits(supertree,"phylo"))
- if(trace>0)message(paste("The MRP supertree, optimized via optim.parsimony(),\nhas a parsimony score of ",
- attr(supertree,"pscore")," (minimum ",characters,")",sep=""))
- else if(inherits(supertree,"multiPhylo"))
- if(trace>0)message(paste("optim.parsimony() found ",length(supertree)," supertrees\nwith a parsimony score of ",
- attr(supertree[[1]],"pscore")," (minimum ",characters,")",sep=""))
- }
- return(supertree)
-}
-
-
-# we want a rooted supertree
-superTree = function(tree, method="pratchet", rooted=TRUE, ...){
- fun = function(x){
- x=reorder(x, "postorder")
- nTips = length(x$tip.label)
- x$edge[x$edge>nTips] = x$edge[x$edge>nTips] + 2L
- l=nrow(x$edge)
- oldroot = x$edge[l,1L]
- x$edge=rbind(x$edge,matrix(c(rep(nTips+2,2),oldroot,nTips+1),2L,2L))
- x$edge.length=c(x$edge.length, 100, 100)
- x$tip.label=c(x$tip.label, "ZZZ")
- x$Nnode=x$Nnode+1L
- x
- }
- if(!is.null(attr(tree, "TipLabel")))tree = .uncompressTipLabel(tree)
- tree = unclass(tree)
- if(rooted) tree = lapply(tree, fun)
- class(tree)="multiPhylo"
- res = my.supertree(tree, method=method, ...)
- if(rooted){
- if(inherits(res,"multiPhylo")){
- res = lapply(res, root, "ZZZ")
- res = lapply(res, drop.tip, "ZZZ")
- class(res) = "multiPhylo"
- }
- else{
- res = root(res, "ZZZ")
- res = drop.tip(res, "ZZZ")
- }
- }
- if(inherits(res,"multiPhylo")){
- fun = function(x){
- x$edge.length <- rep(.1, nrow(x$edge))
- x
- }
- res <- lapply(res, fun)
- res <- lapply(res, reorder, "postorder")
- class(res) = "multiPhylo"
- }
- else{
- res$edge.length = rep(.1, nrow(res$edge))
- res <- reorder(res, "postorder")
- }
- res
-}
-
-
-
densiTree <- function(x, type="cladogram", alpha=1/length(x), consensus=NULL, optim=FALSE, scaleX=FALSE, col=1, width=1, cex=.8, ...) {
if(!inherits(x,"multiPhylo"))stop("x must be of class multiPhylo")
compressed <- ifelse(is.null(attr(x, "TipLabel")), FALSE, TRUE)
diff --git a/R/RcppExports.R b/R/RcppExports.R
new file mode 100644
index 0000000..d540bb3
--- /dev/null
+++ b/R/RcppExports.R
@@ -0,0 +1,11 @@
+# Generated by using Rcpp::compileAttributes() -> do not edit by hand
+# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+allDescCPP <- function(orig, nTips) {
+ .Call('phangorn_allDescCPP', PACKAGE = 'phangorn', orig, nTips)
+}
+
+allChildrenCPP <- function(orig) {
+ .Call('phangorn_allChildrenCPP', PACKAGE = 'phangorn', orig)
+}
+
diff --git a/R/ancestral_pml.R b/R/ancestral_pml.R
index b7331a3..dac7d69 100644
--- a/R/ancestral_pml.R
+++ b/R/ancestral_pml.R
@@ -1,11 +1,11 @@
#
# ancestral sequences ML
#
-ancestral.pml <- function (object, type=c("ml", "bayes"))
+ancestral.pml <- function (object, type=c("marginal", "ml", "bayes"))
{
call <- match.call()
type <- match.arg(type)
- pt <- match.arg(type, c("ml", "bayes"))
+ pt <- match.arg(type, c("marginal", "joint", "ml", "bayes"))
tree = object$tree
INV <- object$INV
@@ -15,7 +15,7 @@ ancestral.pml <- function (object, type=c("ml", "bayes"))
if (is.null(attr(tree, "order")) || attr(tree, "order") ==
"cladewise")
tree <- reorder(tree, "postorder")
- q = length(tree$tip.label)
+ nTips = length(tree$tip.label)
node <- tree$edge[, 1]
edge <- tree$edge[, 2]
m = length(edge) + 1 # max(edge)
@@ -50,7 +50,7 @@ ancestral.pml <- function (object, type=c("ml", "bayes"))
mNodes = as.integer(max(node) + 1)
contrast = attr(data, "contrast")
nco = as.integer(dim(contrast)[1])
- for(i in 1:l)dat[i,(q + 1):m] <- .Call("LogLik2", data, P[i,], nr, nc, node, edge, nTips, mNodes, contrast, nco, PACKAGE = "phangorn")
+ for(i in 1:l)dat[i,(nTips + 1):m] <- .Call("LogLik2", data, P[i,], nr, nc, node, edge, nTips, mNodes, contrast, nco, PACKAGE = "phangorn")
parent <- tree$edge[, 1]
child <- tree$edge[, 2]
@@ -70,7 +70,7 @@ ancestral.pml <- function (object, type=c("ml", "bayes"))
for(i in 1:l){
tmp = tmp + w[i] * dat[[i, j]]
}
- if (pt == "bayes") tmp = tmp * rep(bf, each=nr)
+ if ( (pt == "bayes") || (pt == "marginal")) tmp = tmp * rep(bf, each=nr)
tmp = tmp / rowSums(tmp)
result[[j]] = tmp
}
@@ -80,10 +80,30 @@ ancestral.pml <- function (object, type=c("ml", "bayes"))
}
+#joint_reconstruction <- function(object){
+#
+#}
+
+
+ancestral2phyDat <- function(x) {
+ eps <- 1.0e-5
+ contr <- attr(x, "contrast")
+ # a bit too complicated
+ ind1 <- which( apply(contr, 1, function(x)sum(x > eps)) == 1L)
+ ind2 <- which( contr[ind1, ] > eps, arr.ind = TRUE)
+ pos <- ind2[match(as.integer(1L:ncol(contr)), ind2[,2]),1]
+ # only first hit
+ res <- lapply(x, function(x, pos) pos[apply(x, 1, which.max)], pos)
+ attributes(res) <- attributes(x)
+ return(res)
+}
+
+
fast.tree = function(tree, node){
parent = c(node, Ancestors(tree, node))
children = Descendants(tree, parent, 'children')
- l = sapply(children, length)
+ l = lengths(children)
+# l = sapply(children, length)
edge = cbind(rep(parent, l), unlist(children))
obj = list(edge=edge, Nnode=sum(l>0), tip.label=as.character(edge[is.na(match(edge[,2], edge[,1])),2]))
class(obj) = 'phylo'
diff --git a/R/bootstrap.R b/R/bootstrap.R
index f472f57..7991378 100644
--- a/R/bootstrap.R
+++ b/R/bootstrap.R
@@ -1,4 +1,4 @@
-bootstrap.pml = function (x, bs = 100, trees = TRUE, multicore=FALSE, mc.cores = NULL, ...)
+bootstrap.pml <- function (x, bs = 100, trees = TRUE, multicore=FALSE, mc.cores = NULL, ...)
{
if(multicore && is.null(mc.cores)){
mc.cores <- detectCores()
@@ -81,11 +81,11 @@ bootstrap.phyDat <- function (x, FUN, bs = 100, multicore=FALSE, mc.cores = NULL
}
-matchEdges = function(tree1, tree2){
- bp1 = bip(tree1)
- bp2 = bip(tree2)
- l = length(tree1$tip.label)
- fn = function(x, y){
+matchEdges <- function(tree1, tree2){
+ bp1 <- bip(tree1)
+ bp2 <- bip(tree2)
+ l <- length(tree1$tip.label)
+ fn <- function(x, y){
if(x[1]==1)return(x)
else return(y[-x])
}
@@ -120,10 +120,11 @@ plotBS <- function (tree, BStrees, type = "unrooted", bs.col = "black",
else plot(tree, type = type, ...)
if(hasArg(BStrees)){
- BStrees <- .uncompressTipLabel(BStrees)
- if(any(unlist(lapply(BStrees, is.rooted)))){
- BStrees <- lapply(BStrees, unroot)
- }
+ BStrees <- .uncompressTipLabel(BStrees) # check if needed
+ if(any(is.rooted(BStrees)))BStrees <- unroot(BStrees)
+# if(any(unlist(lapply(BStrees, is.rooted)))){
+# BStrees <- lapply(BStrees, unroot)
+# }
x = prop.clades(tree, BStrees)
x = round((x/length(BStrees)) * 100)
tree$node.label = x
@@ -163,12 +164,10 @@ plotBS <- function (tree, BStrees, type = "unrooted", bs.col = "black",
maxCladeCred <- function(x, tree=TRUE, part=NULL, rooted=TRUE){
if(inherits(x, "phylo")) x <- c(x)
- if(!rooted){
- x <- lapply(x, unroot)
- class(x) <- "multiPhylo"
- x <- .compressTipLabel(x)
- }
- if(is.null(part))pp <- prop.part(x)
+ if (is.null(part)){
+ if (!rooted) pp <- unroot(x) %>% prop.part
+ else pp <- prop.part(x)
+ }
else pp <- part
pplabel <- attr(pp, "labels")
if(!rooted)pp <- oneWise(pp)
@@ -180,6 +179,7 @@ maxCladeCred <- function(x, tree=TRUE, part=NULL, rooted=TRUE){
res <- numeric(l)
for(i in 1:l){
tmp <- checkLabels(x[[i]], pplabel)
+ if(!rooted)tmp <- unroot(tmp)
ppi <- prop.part(tmp) # trees[[i]]
if(!rooted)ppi <- oneWise(ppi)
indi <- fmatch(ppi, pp)
@@ -200,12 +200,13 @@ mcc <- maxCladeCred
cladeMatrix <- function(x, rooted=FALSE){
- if(!rooted){
- x <- .uncompressTipLabel(x)
- x <- lapply(x, unroot)
- class(x) <- "multiPhylo"
- x <- .compressTipLabel(x)
- }
+ if(!rooted) x <- unroot(x)
+# if(!rooted){
+# x <- .uncompressTipLabel(x)
+# x <- lapply(x, unroot)
+# class(x) <- "multiPhylo"
+# x <- .compressTipLabel(x)
+# }
pp <- prop.part(x)
pplabel <- attr(pp, "labels")
if(!rooted)pp <- oneWise(pp)
@@ -232,7 +233,11 @@ cladeMatrix <- function(x, rooted=FALSE){
}
-moving_average <- function(x, window=50){
- cx <- c(0, cumsum(x))
- (cx[(window+1):length(cx)] - cx[1:(length(cx)-window)])/(window)
+moving_average <- function(obj, window=50){
+ fun <- function(x){
+ cx <- c(0, cumsum(x))
+ (cx[(window+1):length(cx)] - cx[1:(length(cx)-window)])/(window)
+ }
+ res <- apply(obj$X, 1, fun)
+ rownames(res) <- c()
}
diff --git a/R/distSeq.R b/R/distSeq.R
index 0c955dc..8269d5a 100644
--- a/R/distSeq.R
+++ b/R/distSeq.R
@@ -118,9 +118,11 @@ dist.ml <- function (x, model = "JC69", exclude = "none", bf = NULL, Q = NULL, k
}
tmp = (contrast %*% eig[[2]])[ind1, ] * (contrast %*% (t(eig[[3]]) * bf))[ind2, ]
tmp2 = vector("list", k)
- wdiag = .Call("PWI", as.integer(1:n), as.integer(1:n), as.integer(n),
- as.integer(n), rep(1, n), as.integer(li), PACKAGE = "phangorn")
- wdiag = which(wdiag > 0)
+# wdiag = .Call("PWI", as.integer(1:n), as.integer(1:n), as.integer(n),
+# as.integer(n), rep(1, n), as.integer(li), PACKAGE = "phangorn")
+# wdiag = which(wdiag > 0)
+
+ wshared <- which(rowSums(contrast[ind1, ] * contrast[ind2, ]) > 0)
tmp2 = vector("list", k)
for (i in 1:(l - 1)) {
@@ -131,7 +133,7 @@ dist.ml <- function (x, model = "JC69", exclude = "none", bf = NULL, Q = NULL, k
w0[index] = 0.0
ind = w0 > 0
- old.el <- 1 - (sum(w0[wdiag])/sum(w0))
+ old.el <- 1 - (sum(w0[wshared])/sum(w0))
if (old.el > eps)
old.el <- 10
else old.el <- fun(old.el)
@@ -141,7 +143,7 @@ dist.ml <- function (x, model = "JC69", exclude = "none", bf = NULL, Q = NULL, k
for(lk in 1:k) tmp2[[lk]] <- tmp[ind, , drop = FALSE]
# FS0 verwenden!!!
res <- .Call("FS5", eig, nc, as.double(old.el), w, g, tmp2, as.integer(k), as.integer(sum(ind)),
- bf, w0[ind], ll.0[ind], PACKAGE = "phangorn")
+ w0[ind], ll.0, PACKAGE = "phangorn") #bf,
d[pos] <- res[1] # res[[1]]
v[pos] <- res[2] # res[[2]]
pos = pos + 1
@@ -159,91 +161,6 @@ dist.ml <- function (x, model = "JC69", exclude = "none", bf = NULL, Q = NULL, k
return(d)
}
-
-dist.mlOld <- function (x, model = "JC69", exclude = "none", bf = NULL, Q = NULL, ...)
-{
- if (!inherits(x,"phyDat"))
- stop("x has to be element of class phyDat")
- l = length(x)
- d = numeric((l * (l - 1))/2)
- v = numeric((l * (l - 1))/2)
- contrast <- attr(x, "contrast")
- nc <- as.integer(attr(x, "nc"))
- nr <- as.integer(attr(x, "nr"))
- con = rowSums(contrast > 0) < 2
- if (exclude == "all") {
- index = con[x[[1]]]
- for (i in 2:l) index = index & con[x[[i]]]
- index = which(index)
- x = subset(x, , index)
- }
-# model <- match.arg(model, c("JC69", "WAG", "JTT", "LG", "Dayhoff", "cpREV", "mtmam", "mtArt", "MtZoa", "mtREV24"))
- model <- match.arg(model, c("JC69", .aamodels))
-# if (!is.na(match(model, c("WAG", "JTT", "LG", "Dayhoff", "cpREV", "mtmam", "mtArt", "MtZoa", "mtREV24"))))
- if (!is.na(match(model, .aamodels)))
- getModelAA(model, bf = is.null(bf), Q = is.null(Q))
- if (is.null(bf))
- bf <- rep(1/nc, nc)
- if (is.null(Q))
- Q <- rep(1, (nc - 1) * nc/2L)
-
- bf = as.double(bf)
- eig <- edQt(Q = Q, bf = bf)
- k = 1
- w = as.double(1)
- g = as.double(1)
- fun <- function(s) -(nc - 1)/nc * log(1 - nc/(nc - 1) * s)
- eps <- (nc - 1)/nc
- n = as.integer(dim(contrast)[1])
- ind1 = rep(1:n, n:1)
- ind2 = unlist(lapply(n:1, function(x) seq_len(x) + n - x))
- li <- as.integer(length(ind1))
- weight = as.double(attr(x, "weight"))
- ll.0 = as.double(weight * 0)
- if (exclude == "pairwise") {
- index = con[ind1] & con[ind2]
- index = which(!index)
- }
- tmp = (contrast %*% eig[[2]])[ind1, ] * (contrast %*% (t(eig[[3]]) * bf))[ind2, ]
- tmp2 = vector("list", k)
- wdiag = .Call("PWI", as.integer(1:n), as.integer(1:n), as.integer(n),
- as.integer(n), rep(1, n), as.integer(li), PACKAGE = "phangorn")
- wdiag = which(wdiag > 0)
- for (i in 1:(l - 1)) {
- for (j in (i + 1):l) {
- w0 = .Call("PWI", as.integer(x[[i]]), as.integer(x[[j]]),
- nr, n, weight, li, PACKAGE = "phangorn")
- if (exclude == "pairwise")
- w0[index] = 0.0
- ind = w0 > 0
-
- old.el <- 1 - (sum(w0[wdiag])/sum(w0))
- if (old.el > eps)
- old.el <- 10
- else old.el <- fun(old.el)
- # sind = sum(ind)
- # tmp2 = vector("list", k)
- tmp2[[1]] <- tmp[ind, , drop = FALSE]
- # FS0 verwenden!!!
- res <- .Call("FS5", eig, nc, as.double(old.el), w, g, tmp2, 1L, as.integer(sum(ind)),
- bf, w0[ind], ll.0, PACKAGE = "phangorn")
- d[k] <- res[1] # res[[1]]
- v[k] <- res[2] # res[[2]]
- k = k + 1
- }
- }
- attr(d, "Size") <- l
- if (is.list(x))
- attr(d, "Labels") <- names(x)
- else attr(d, "Labels") <- colnames(x)
- attr(d, "Diag") <- FALSE
- attr(d, "Upper") <- FALSE
- attr(d, "call") <- match.call()
- attr(d, "variance") <- v
- class(d) <- "dist"
- return(d)
-}
-
dist.logDet = function (x)
{
diff --git a/R/distTree.R b/R/distTree.R
index 2667649..ca5d097 100644
--- a/R/distTree.R
+++ b/R/distTree.R
@@ -185,7 +185,7 @@ UNJ <- function(x)
result = list(edge = rbind(c(nam[2], nam[1]), edge),
edge.length=edge.length, tip.label = labels, Nnode=m)
class(result) <- "phylo"
- reorder(result)
+ reorder(result, "postorder")
}
@@ -323,7 +323,7 @@ designUltra <- function (tree, sparse=TRUE)
v=numeric( n * (n - 1)/2)
m = 1L
for (i in 1:length(leri)) {
- if (!is.null(leri[[i]])) {
+ if (length(leri[[i]])>1) {
if(length(leri[[i]])==2)ind = getIndex(bp[[leri[[i]][1] ]], bp[[leri[[i]][2] ]], n)
else {
ind=NULL
@@ -352,7 +352,7 @@ designUltra <- function (tree, sparse=TRUE)
designUnrooted2 <- function (tree, sparse=TRUE)
{
- if (is.null(attr(tree, "order")) || attr(tree, "order") == "cladewise")
+ if (is.null(attr(tree, "order")) || attr(tree, "order") != "postorder")
tree = reorder(tree, "postorder")
leri = allChildren(tree)
bp = bip(tree)
@@ -368,8 +368,7 @@ designUnrooted2 <- function (tree, sparse=TRUE)
p=1L
m = 1L
for (i in 1:length(leri)) {
- if (!is.null(leri[[i]])) {
-
+ if (length(leri[[i]]) > 1) {
if(length(leri[[i]])==2){
ind = getIndex(bp[[leri[[i]][1] ]], bp[[leri[[i]][2] ]], n)
ytmp = rep(bp[[leri[[i]][1] ]], each = length(bp[[leri[[i]][2] ]]))
@@ -456,11 +455,11 @@ nnls.tree <- function(dm, tree, rooted=FALSE, trace=1){
betahat = bhat[tree$edge[,1]] - bhat[tree$edge[,2]]
if(!any(betahat<0)){
- if(!rooted){
+# if(!rooted){
RSS = sum((y-(X%*%betahattmp))^2)
if(trace)print(paste("RSS:", RSS))
attr(tree, "RSS") = RSS
- }
+# }
tree$edge.length = betahat
return(tree)
}
@@ -468,23 +467,42 @@ nnls.tree <- function(dm, tree, rooted=FALSE, trace=1){
# non-negative LS
n = dim(X)[2]
l = nrow(tree$edge)
- Amat = matrix(0, n, l)
+
lab = attr(X, "nodes")
# vielleicht solve.QP.compact
ind1 = match(tree$edge[,1], lab)
- Amat[cbind(ind1, 1:l)] = 1
ind2 = match(tree$edge[,2], lab)
- Amat[cbind(ind2, 1:l)] = -1
- betahat <- quadprog::solve.QP(as.matrix(Dmat),as.vector(dvec),Amat)$sol
+# Amat = matrix(0, n, l)
+# Amat[cbind(ind1, 1:l)] = 1
+# Amat[cbind(ind2, 1:l)] = -1
+# betahat <- quadprog::solve.QP(as.matrix(Dmat),as.vector(dvec),Amat)$sol
+
+ Amat <- matrix(0, 2, l)
+ Amat[1,] <- 1
+ Amat[2,] <- -1
+
+ Aind <- matrix(0L, 3, l)
+ Aind[1,] <- 2L
+ Aind[2,] <- as.integer(ind1)
+ Aind[3,] <- as.integer(ind2)
+
+ if(any(is.na(Aind))){
+ na_ind <- which(is.na(Aind), arr.ind = TRUE)
+ Aind[is.na(Aind)] <- 0L
+ for(i in 1:nrow(na_ind)) Aind[1, na_ind[i, 2]] <- Aind[1, na_ind[i, 2]] - 1L
+ }
+
+ betahat <- quadprog::solve.QP.compact(as.matrix(Dmat),as.vector(dvec),Amat, Aind)$sol
+
# quadratic programing solving
- if(!rooted){
+# if(!rooted){
RSS = sum((y-(X%*%betahat))^2)
if(trace)print(paste("RSS:", RSS))
attr(tree, "RSS") = RSS
- }
+# }
bhat = numeric(max(tree$edge))
bhat[as.integer(lab)] = betahat
betahat = bhat[tree$edge[,1]] - bhat[tree$edge[,2]]
@@ -506,40 +524,7 @@ nnls.splits <- function(x, dm, trace=0){
y = dm[lower.tri(dm)]
x = SHORTwise(x, k)
- l <- sapply(x, length)
- if(any(l==0)) x = x[-which(l==0)]
-
- X = splits2design(x)
-
- if(any(is.na(y))){
- ind = which(is.na(y))
- X = X[-ind,,drop=FALSE]
- y= y[-ind]
- }
- X = as.matrix(X)
- n = dim(X)[2]
- int = sapply(x, length)
- Amat = diag(n) # (int)
- betahat <- nnls(X, y)
- ind = (betahat$x > 1e-8) | int==1
- x = x[ind]
- RSS <- betahat$deviance
- attr(x, "weights") = betahat$x[ind]
- if(trace)print(paste("RSS:", RSS))
- attr(x, "RSS") = RSS
- x
-}
-
-
-nnls.splitsOld <- function(x, dm, trace=0){
- labels=attr(x, "labels")
- dm = as.matrix(dm)
- k = dim(dm)[1]
- dm = dm[labels,labels]
- y = dm[lower.tri(dm)]
-
- x = SHORTwise(x, k)
- l <- sapply(x, length)
+ l <- lengths(x)
if(any(l==0)) x = x[-which(l==0)]
X = splits2design(x)
@@ -562,11 +547,17 @@ nnls.splitsOld <- function(x, dm, trace=0){
return(x)
}
n = dim(X)[2]
+ int <- lengths(x)
+
+# Amat = diag(n) # (int)
+# betahat <- quadprog::solve.QP(as.matrix(Dmat),as.vector(dvec),Amat)$sol # quadratic programing solving
+# faster version
+ Amat <- matrix(1, 1, n)
+ Aind <- matrix(0L, 2L, n)
+ Aind[1,] <- 1L
+ Aind[2,] <- as.integer(1L:n)
+ betahat <- quadprog::solve.QP.compact(as.matrix(Dmat),as.vector(dvec),Amat, Aind)$sol
- int = sapply(x, length)
- # int = as.numeric(int==1)# (int>1)
- Amat = diag(n) # (int)
- betahat <- quadprog::solve.QP(as.matrix(Dmat),as.vector(dvec),Amat)$sol # quadratic programing solving
RSS = sum((y-(X%*%betahat))^2)
ind = (betahat > 1e-8) | int==1
x = x[ind]
@@ -601,7 +592,7 @@ designSplits <- function (x, splits = "all", ...)
if (is.na(splits)) stop("invalid splits method")
if (splits == -1) stop("ambiguous splits method")
if(splits==1) X <- designAll(x)
- if(splits==2) X <- designStar(x)
+ if(splits==2) X <- designStar(x, ...)
return(X)
}
@@ -625,9 +616,12 @@ designAll <- function(n, add.split=FALSE){
}
-designStar = function(n){
- res=NULL
- for(i in 1:(n-1)) res = rbind(res,cbind(matrix(0,(n-i),i-1),1,diag(n-i)))
+# faster sparse version
+designStar = function(n, sparse=TRUE){
+# res=NULL
+# for(i in 1:(n-1)) res = rbind(res,cbind(matrix(0,(n-i),i-1),1,diag(n-i)))
+ res <- stree(n) %>% as.splits %>% splits2design
+ if(!sparse) return(as.matrix(res))
res
}
diff --git a/R/fitch.R b/R/fitch.R
index c3ee924..c47fc18 100644
--- a/R/fitch.R
+++ b/R/fitch.R
@@ -503,7 +503,7 @@ bab_old <- function (data, tree = NULL, trace = 1, ...)
ind = (1:L[a])[score<=bound]
for(i in 1:length(ind))trees[[a+1]][[i]] <- .Call("AddOne", tmpTree, as.integer(inord[a+1L]), as.integer(ind[i]), as.integer(L[a]), as.integer(M[a]), PACKAGE="phangorn")
l = length(ind)
- os = order(score[ind], decreasing=TRUE)
+ os = order(score[ind], decreasing=TRUE)
PSC = rbind(PSC, cbind(rep(a+1, l), os, score[ind][os] ))
}
else{
@@ -535,6 +535,7 @@ bab_old <- function (data, tree = NULL, trace = 1, ...)
bab <- function (data, tree = NULL, trace = 1, ...)
{
+ if(!is.null(tree)) data <- subset(data, tree$tip.label)
o = order(attr(data, "weight"), decreasing = TRUE)
data = subset(data, , o)
nr <- attr(data, "nr")
diff --git a/R/networx.R b/R/networx.R
index 0d24075..c86f58b 100644
--- a/R/networx.R
+++ b/R/networx.R
@@ -28,7 +28,8 @@ as.Matrix.splits <- function(x, ...){
labels = attr(x, "labels")
l = length(x)
j = unlist(x)
- i = rep(1:l, sapply(x, length))
+# i = rep(1:l, sapply(x, length))
+ i = rep(1:l, lengths(x))
sparseMatrix(i,j, x = rep(1L, length(i)), dimnames = list(NULL, labels)) # included x und labels
}
@@ -167,6 +168,7 @@ countCycles <- function(splits, tree=NULL, ord=NULL){
c.splits <- function (..., recursive=FALSE)
{
x <- list(...)
+ if (length(x) == 1 && !inherits(x[[1]], "splits")) x <- x[[1]]
n <- length(x)
match.names <- function(a, b) {
if (any(!(a %in% b)))
@@ -179,11 +181,25 @@ c.splits <- function (..., recursive=FALSE)
cycle <- attr(x[[1]], "cycle")
for (i in 2:n) {
match.names(labels, attr(x[[i]], "labels"))
+ x[[i]] <- changeOrder(x[[i]], labels)
}
- res = structure(NextMethod("c"), class=c("splits", "prop.part"))
- attr(res, "labels") = labels
- attr(res, "weight") = as.vector(sapply(x, attr, "weight"))
- attr(res, "cycle") = cycle
+ w <- as.vector(sapply(x, attr, "weights"))
+ x <- lapply(x, unclass)
+ res <- structure(do.call("c", x), class=c("splits", "prop.part"))
+ names(res) <- NULL
+# res <- structure(NextMethod("c"), class=c("splits", "prop.part"))
+ attr(res, "labels") <- labels
+ attr(res, "weights") <- w
+ attr(res, "cycle") <- cycle
+ res
+}
+
+
+distinct.splits <- function(...){
+ tmp <- c(...)
+ res <- unique(tmp)
+ attributes(res) <- attributes(tmp)
+ attr(res, "weights") <- tabulate(match(tmp, res))
res
}
@@ -219,8 +235,9 @@ as.splits.multiPhylo <- function(x, ...){
# firstTip = x[[1]]$tip[1]
# x = lapply(x, root, firstTip) # old trick
lx <- length(x)
- x <- lapply(x, unroot)
- class(x) <- "multiPhylo"
+ x <- unroot(x)
+# lapply(x, unroot)
+# class(x) <- "multiPhylo"
splits <- prop.part(x)
class(splits)='list'
weights = attr(splits, 'number')
@@ -422,7 +439,8 @@ splitsNetwork <- function(dm, splits=NULL, gamma=.1, lambda=1e-6, weight=NULL){
k = dim(dm)[1]
if(!is.null(splits)){
- tmp = which(sapply(splits, length)==k)
+# tmp = which(sapply(splits, length)==k)
+ tmp = which(lengths(splits)==k)
splits = splits[-tmp]
lab = attr(splits, "labels")
dm = dm[lab, lab]
@@ -436,7 +454,8 @@ splitsNetwork <- function(dm, splits=NULL, gamma=.1, lambda=1e-6, weight=NULL){
y = dm[lower.tri(dm)]
if(is.null(splits))ind = c(2^(0:(k-2)),2^(k-1)-1)
- else ind = which(sapply(splits, length)==1)
+ else ind = which(lengths(splits)==1)
+# else ind = which(sapply(splits, length)==1)
# y2 = lm(y~X[,ind]-1)$res
n = dim(X)[2]
@@ -456,6 +475,7 @@ splitsNetwork <- function(dm, splits=NULL, gamma=.1, lambda=1e-6, weight=NULL){
Amat <- cbind(ind1,diag(n))
bvec <- c(gamma, rep(0,n))
+# needs quadprog::solve.QP.compact
solution <- quadprog::solve.QP(Dmat,dvec,Amat,bvec=bvec, meq=1)$sol
ind2 <- which(solution>1e-8)
@@ -472,6 +492,8 @@ splitsNetwork <- function(dm, splits=NULL, gamma=.1, lambda=1e-6, weight=NULL){
Amat2 <- diag(n2)
bvec2 <- rep(0, n2)
+ # needs quadprog::solve.QP.compact
+ # bvec2 not used
solution2 <- quadprog::solve.QP(Dmat, dvec, Amat2)$sol
RSS1 = sum((y-X[,ind2]%*%solution[ind2])^2)
@@ -512,9 +534,11 @@ allCircularSplits <- function(k, labels=NULL){
tmp = (1L:y)+x
tmp %% (k+1L) + tmp %/% (k+1L)
}
- for(i in 2:l){
- res[(ind+1):(ind+k)] <- lapply(0L:(k-1L), fun, i)
- ind <- ind+k
+ if(k>4L){
+ for(i in 2:l){
+ res[(ind+1):(ind+k)] <- lapply(0L:(k-1L), fun, i)
+ ind <- ind+k
+ }
}
if((k%%2L)==0){
m <- k%/%2
@@ -522,8 +546,9 @@ allCircularSplits <- function(k, labels=NULL){
}
}
- if(is.null(labels)) labels=(as.character(1:k))
- attr(res, 'labels') =labels
+ if(is.null(labels)) labels <- as.character(1:k)
+ attr(res, 'labels') <- labels
+ attr(res, "cycle") <- 1:k
class(res)="splits"
res
}
@@ -544,7 +569,8 @@ splits2design <- function(obj, weight=NULL){
m = length(labels)
n=length(obj)
l = 1:m
- sl = sapply(obj, length)
+ sl = lengths(obj)
+# sl = sapply(obj, length)
p0 = sl * (m-sl)
p = c(0,cumsum(p0))
i = numeric(max(p))
@@ -641,9 +667,9 @@ addEdge <- function(network, desc, spl){
#neu zu alt verbinden
edge = rbind(edge, cbind(oldNodes, newNodes), deparse.level = 0)
index = c(index, rep(spl, length(oldNodes)) )
- network$splitIndex = index
network$edge = edge
network$Nnode = max(edge) - nTips
+ network$splitIndex = index
network
}
@@ -670,14 +696,18 @@ circNetwork <- function(x, ord=NULL){
index = match(spRes, x)
}
- l = sapply(oneWise(x, nTips), length)
- l2 = sapply(x, length)
+ l = lengths(oneWise(x, nTips))
+ l2 = lengths(x)
+# l = sapply(oneWise(x, nTips), length)
+# l2 = sapply(x, length)
+
# dm <- as.matrix(compatible2(x))
tmp <- countCycles(x, ord=ord)
ind = which(tmp == 2 & l2>1) # & l<nTips changed with ordering
- ind = ind[order(l[ind])]
+# ind = ind[order(l[ind])]
+ ind = ind[order(l2[ind], decreasing = TRUE)]
dm2 <- as.matrix(compatible2(x, x[ind]))
@@ -770,8 +800,8 @@ circNetwork <- function(x, ord=NULL){
class(res) = c("networx", "phylo")
attr(res, "order") = NULL
}
- res$Nnode = max(res$edge) - nTips
res$edge.length = weight[index] # ausserhalb
+ res$Nnode = max(res$edge) - nTips
res$splitIndex = index
# attr(res, "splits") = x
res$splits = x
@@ -830,7 +860,8 @@ as.networx.splits <- function(x, planar=FALSE, coord = c("none", "2D", "3D"), ..
attr(x, "weights") <- weight
x <- oneWise(x, nTips)
- l <- sapply(x, length)
+ l <- lengths(x)
+# l <- sapply(x, length)
if(any(l==nTips))x <- x[l!=nTips] # get rid of trivial splits
ext <- sum(l==1 | l==(nTips-1))
if(!is.null(attr(x, "cycle"))){
@@ -849,7 +880,8 @@ as.networx.splits <- function(x, planar=FALSE, coord = c("none", "2D", "3D"), ..
return(reorder(tmp))
}
- ll <- sapply(x, length)
+ ll <- lengths(x)
+# ll <- sapply(x, length)
ind <- tmp$splitIndex # match(sp, x)
ind2 = union(ind, which(ll==0)) # which(duplicated(x))
ind2 = union(ind2, which(ll==nTips))
@@ -863,8 +895,8 @@ as.networx.splits <- function(x, planar=FALSE, coord = c("none", "2D", "3D"), ..
class(tmp) = c("networx", "phylo")
}
}
- tmp$Nnode = max(tmp$edge) - nTips
tmp$edge.length = weight[tmp$splitIndex]
+ tmp$Nnode = max(tmp$edge) - nTips
attr(x, "cycle") <- c.ord
# attr(tmp, "splits") = x
tmp$splits <- x
@@ -891,7 +923,8 @@ as.networx.phylo <- function(x, ...){
spl <- as.splits(x)
spl <- spl[x$tree[,2]]
x$splitIndex <- 1:nrow(x$edge)
- attr(x, "splits") = spl
+# attr(x, "splits") = spl
+ x$splits <- spl
class(x) <- c("networx", "phylo")
x
}
@@ -1034,7 +1067,8 @@ coords <- function(obj, dim="3D"){
# g2 <- graph(t(obj$edge), directed=FALSE)
# g2 <- set.edge.attribute(g, "weight", value=rep(1, nrow(obj$edge))
if(dim=="3D"){
- coord <- layout_with_kk(g, dim=3)
+ coord <- layout_nicely(g, dim=3)
+# coord <- layout_with_kk(g, dim=3)
# coord <- layout.kamada.kawai(g, dim=3)
k = matrix(0, max(obj$splitIndex), 3)
for(i in ind1){
@@ -1052,7 +1086,8 @@ coords <- function(obj, dim="3D"){
}
}
else{
- coord <- layout_with_kk(g, dim=2)
+ coord <- layout_nicely(g, dim=2)
+# coord <- layout_with_kk(g, dim=2)
# coord <- layout.kamada.kawai(g, dim=2)
k = matrix(0, max(obj$splitIndex), 2)
for(i in ind1){
@@ -1320,7 +1355,8 @@ lento <- function (obj, xlim = NULL, ylim = NULL, main = "Lento plot",
labels = attr(obj, "labels")
l = length(labels)
if(!trivial){
- triv = sapply(obj, length)
+ triv = lengths(obj)
+# triv = sapply(obj, length)
ind = logical(length(obj))
ind[(triv >1) & (triv < (l-1))] = TRUE
if(length(col)==length(obj)) col=col[ind]
@@ -1412,12 +1448,14 @@ write.nexus.splits <- function (obj, file = "", weights=NULL, taxa=TRUE, append=
taxa.labels <- attr(obj, "labels")
ntaxa <- length(taxa.labels)
obj <- oneWise(obj, ntaxa)
- ind <- which(sapply(obj, length)==ntaxa)
+ ind <- which(lengths(obj)==ntaxa)
+# ind <- which(sapply(obj, length)==ntaxa)
if(length(ind))obj <- obj[-ind]
nsplits <- length(obj)
if(is.null(weights))weight <- attr(obj, "weights")
- if (is.null(weight))
- weight = numeric(nsplits) + 100
+ if (is.null(weight)) fwei <- FALSE
+ else fwei <- TRUE
+# if (is.null(weight)) weight = numeric(nsplits) + 100
if(!append){cat("#NEXUS\n\n", file = file)
cat("[Splits block for Spectronet or Splitstree]\n", file = file, append = TRUE)
cat("[generated by phangorn:\n", file = file, append = TRUE)
@@ -1433,7 +1471,10 @@ write.nexus.splits <- function (obj, file = "", weights=NULL, taxa=TRUE, append=
# SPLITS BLOCK
cat(paste("BEGIN SPLITS;\n\tDIMENSIONS ntax=", ntaxa, " nsplits=", nsplits,
";\n", sep = ""), file = file, append = TRUE)
- format = "\tFORMAT labels=left weights=yes"
+ format = "\tFORMAT labels=left"
+ # weights=yes
+ if(fwei) format = paste(format, "weights=yes")
+ else format = paste(format, "weights=no")
fcon = fint = flab = FALSE
if(!is.null(attr(obj, "confidences"))){
format = paste(format, "confidences=yes")
@@ -1461,10 +1502,13 @@ write.nexus.splits <- function (obj, file = "", weights=NULL, taxa=TRUE, append=
for (i in 1:nsplits){
slab <- ifelse(flab, attr(obj, "splitlabels")[i], i)
+ swei <- ifelse(fwei, paste(weight[i], "\t"), "")
scon <- ifelse(fcon, paste(attr(obj, "confidences")[i], "\t"), "")
sint <- ifelse(fint, paste(attr(obj, "intervals")[i], "\t"), "")
- cat("\t\t", slab, "\t", weight[i], "\t", scon, sint, paste(obj[[i]], collapse=" "),
- ",\n", file = file, append = TRUE, sep = "")
+# cat("\t\t", slab, "\t", weight[i], "\t", scon, sint, paste(obj[[i]], collapse=" "),
+# ",\n", file = file, append = TRUE, sep = "")
+ cat("\t\t", slab, "\t", swei, scon, sint, paste(obj[[i]], collapse=" "),
+ ",\n", file = file, append = TRUE, sep = "")
}
cat("\t;\nEND;\n", file = file, append = TRUE)
}
@@ -1693,13 +1737,13 @@ read.nexus.networx <- function(file, splits=TRUE){
vert <- swapRow(vert, oldLabel[i], i)
}
}
-
+ dimnames(edge) <- NULL
# y-axis differs between in R and SplitsTree
vert[,2] <- -vert[,2]
# translate=data.frame(as.numeric(TRANS[,1]), TRANS[,2], stringsAsFactors=FALSE)
plot <- list(vertices=vert)
- obj <- list(edge=edge, tip.label=TRANS$label, Nnode=max(edge)-ntaxa,
- edge.length=el, splitIndex=splitIndex, splits=spl, translate=TRANS)
+ obj <- list(edge=edge, tip.label=TRANS$label, edge.length=el, Nnode=max(edge)-ntaxa,
+ splitIndex=splitIndex, splits=spl, translate=TRANS)
obj$.plot <- list(vertices = vert, edge.color="black", edge.width=3, edge.lty = 1)
class(obj) <- c("networx", "phylo")
reorder(obj)
@@ -1778,12 +1822,12 @@ read.nexus.splits <- function(file)
cyc = as.integer(na.omit(as.numeric(strsplit(tmp, " ")[[1]])))
}
attr(res, "labels") = x
- attr(res, "weights") = weights
+ if(fwei)attr(res, "weights") = weights
if(fint)attr(res, "intervals") = intervals
if(fcon)attr(res, "confidences") = confidences
if(flab)attr(res, "splitlabels") = labels
attr(res, "cycle") = cyc
- class(res) = c("splits", "prop.part")
+ class(res) = "splits"
res
}
diff --git a/R/parsimony.R b/R/parsimony.R
index 92c820d..65e87b1 100644
--- a/R/parsimony.R
+++ b/R/parsimony.R
@@ -246,10 +246,16 @@ lowerBound <- function(x, cost=NULL){
y <- as.character(x)
states <- apply(y, 2, unique.default)
+
# duplicated function(x)x[duplicated(x)]="?" avoids looping
if(nr==1) nst <- length(states)
- else nst <- sapply(states, length)
-
+ else{
+ if(is.matrix(states)){
+ states = as.data.frame(states, stringsAsFactors=FALSE)
+ class(states) = "list"
+ }
+ nst <- sapply(states, length)
+ }
res = numeric(nr)
ust = sort(unique(nst))
@@ -647,7 +653,8 @@ optim.parsimony <- function(tree,data, method='fitch', cost=NULL, trace=1, rearr
}
-pratchet <- function (data, start=NULL, method="fitch", maxit=1000, k=10, trace=1, all=FALSE, rearrangements="SPR", ...)
+# perturbation="ratchet", "stochastic"
+pratchet <- function (data, start=NULL, method="fitch", maxit=1000, k=10, trace=1, all=FALSE, rearrangements="SPR", perturbation="ratchet", ...)
{
eps = 1e-08
# if(method=="fitch" && (is.null(attr(data, "compressed")) || attr(data, "compressed") == FALSE))
@@ -689,9 +696,18 @@ pratchet <- function (data, start=NULL, method="fitch", maxit=1000, k=10, trace=
result = list()
result[[1]] = tree
kmax = 1
+ nTips = length(tree$tip.label)
for (i in 1:maxit) {
- bstrees <- bootstrap.phyDat(data, FUN, tree = tree, bs = 1, trace = trace, method=method, rearrangements=rearrangements, ...)
- trees <- lapply(bstrees, optim.parsimony, data, trace = trace, method=method, rearrangements=rearrangements, ...)
+ if(perturbation=="ratchet"){
+ bstrees <- bootstrap.phyDat(data, FUN, tree = tree, bs = 1, trace = trace, method=method, rearrangements=rearrangements, ...)
+ trees <- lapply(bstrees, optim.parsimony, data, trace = trace, method=method, rearrangements=rearrangements, ...)
+ }
+ if(perturbation=="stochastic"){
+ treeNNI <- rNNI(tree, floor(nTips/2))
+ trees <- optim.parsimony(treeNNI, data, trace = trace,
+ method = method, rearrangements = rearrangements, ...)
+ trees <- list(trees)
+ }
if(inherits(result,"phylo"))m=1
else m = length(result)
if(m>0) trees[2 : (1+m)] = result[1:m]
@@ -766,9 +782,9 @@ ptree <- function (tree, data, type = "ACCTRAN", retData = FALSE)
stop("data must be of class phyDat")
if (is.null(attr(tree, "order")) || attr(tree, "order") ==
"cladewise")
- tree <- reorder(tree, "pruningwise")
- # if (!is.binary.tree(tree))
- # stop("Tree must be binary!")
+ tree <- reorder(tree, "postorder")
+# if (!is.binary.tree(tree))
+# stop("Tree must be binary!")
tmp = fitch(tree, data, site = "data")
nr = attr(data, "nr")
node <- tree$edge[, 1]
@@ -786,17 +802,17 @@ ptree <- function (tree, data, type = "ACCTRAN", retData = FALSE)
dat[, ind[2]], dat[, ind[3]], as.integer(nr))
dat[, root] = rSeq[[1]]
}
- result <- .C("ACCTRAN2", dat, as.integer(nr), numeric(nr),
- as.integer(node[l:1L]), as.integer(edge[l:1L]), l, as.double(weight),
- numeric(l), as.integer(nTips))
- el = result[[8]][l:1L]
+# result <- .C("ACCTRAN2", dat, as.integer(nr), numeric(nr),
+# as.integer(node[l:1L]), as.integer(edge[l:1L]), l, as.double(weight), numeric(l), as.integer(nTips))
+ result <- .C("ACCTRAN2", dat, as.integer(nr),
+ as.integer(node[l:1L]), as.integer(edge[l:1L]), l, as.integer(nTips))
+ el = result[[5]][l:1L]
if (!is.rooted(tree)) {
ind2 = which(node[] == root)
dat = matrix(result[[1]], nr, max(node))
result <- .C("ACCTRAN3", result[[1]], as.integer(nr),
numeric(nr), as.integer(node[(l - 3L):1L]), as.integer(edge[(l -
- 3L):1L]), l - 3L, as.double(weight), numeric(l),
- as.integer(nTips))
+ 3L):1L]), l - 3L, as.double(weight), numeric(l)) # , as.integer(nTips)
el = result[[8]][(l - 3L):1L]
pars = .C("fitchTripletACC4", dat[, root], dat[, ind[1]],
dat[, ind[2]], dat[, ind[3]], as.integer(nr), numeric(1),
@@ -809,7 +825,7 @@ ptree <- function (tree, data, type = "ACCTRAN", retData = FALSE)
else {
result <- .C("ACCTRAN3", result[[1]], as.integer(nr),
numeric(nr), as.integer(node[l:1L]), as.integer(edge[l:1L]),
- l, as.double(weight), numeric(l), as.integer(nTips))
+ l, as.double(weight), numeric(l)) # , as.integer(nTips)
el = result[[8]][l:1L]
}
tree$edge.length = el
diff --git a/R/phyDat.R b/R/phyDat.R
index 274e1b0..f428768 100644
--- a/R/phyDat.R
+++ b/R/phyDat.R
@@ -15,26 +15,27 @@ fast.table <- function (data)
bin <- bin[!is.na(bin)]
if (length(bin)) bin <- bin + 1
y <- tabulate(bin, pd)
- result=list(index = bin, weights = y, data = data[ind,])
+ result=list(index = bin, weights = y, data = data[ind,])
result
}
phyDat.default <- function (data, levels = NULL, return.index = TRUE, contrast = NULL,
- ambiguity = "?", compress=TRUE, ...)
+ ambiguity = "?", compress=TRUE, ...)
{
if (is.matrix(data))
nam = row.names(data)
else nam = names(data)
if(is.null(nam))stop("data object must contain taxa names")
+ if(inherits(data, "character")) data <- as.matrix(data)
if (inherits(data,"DNAbin"))
data = as.character(data)
if (is.matrix(data))
data = as.data.frame(t(data), stringsAsFactors = FALSE)
-# new 4.4.2016 bug fix (reported by Eli Levy Karin)
- if (is.vector(data) && !is.list(data))data = as.data.frame(t(data), stringsAsFactors = FALSE)
+ # new 4.4.2016 bug fix (reported by Eli Levy Karin)
+ # if (is.vector(data) && !is.list(data))data = as.data.frame(data, stringsAsFactors = FALSE)
else data = as.data.frame(data, stringsAsFactors = FALSE)
-# data = data.frame(as.matrix(data), stringsAsFactors = FALSE)
+ # data = data.frame(as.matrix(data), stringsAsFactors = FALSE)
if(length(data[[1]])==1) compress=FALSE
@@ -70,10 +71,26 @@ phyDat.default <- function (data, levels = NULL, return.index = TRUE, contrast =
contrast = rbind(contrast, matrix(1, k, l))
}
}
+
+
+# row.names(data) = as.character(1:p)
+# data = na.omit(data)
+# rn = as.numeric(rownames(data))
+
+ d = dim(data)
att = attributes(data)
- data = lapply(data, match, all.levels) # avoid unlist
+
+# print(system.time(match(unlist(data), all.levels)))
+
+ data = match(unlist(data), all.levels)
+ attr(data, "dim") = d
+ data = as.data.frame(data, stringsAsFactors=FALSE)
attributes(data) = att
+ # att = attributes(data)
+ # data = lapply(data, match, all.levels) # avoid unlist
+ # attributes(data) = att
+
row.names(data) = as.character(1:p)
data = na.omit(data)
@@ -100,6 +117,7 @@ phyDat.default <- function (data, levels = NULL, return.index = TRUE, contrast =
}
+
phyDat.DNA = function (data, return.index = TRUE)
{
if (is.matrix(data))
@@ -107,6 +125,7 @@ phyDat.DNA = function (data, return.index = TRUE)
else nam = names(data)
if (inherits(data,"DNAbin"))
data = as.character(data)
+ if(inherits(data, "character")) data <- as.matrix(data)
if (is.matrix(data))
data = as.data.frame(t(data), stringsAsFactors = FALSE)
else data = as.data.frame(data, stringsAsFactors = FALSE)
@@ -122,9 +141,22 @@ phyDat.DNA = function (data, return.index = TRUE)
0, 1, 1, 1, 1, 1, 1)), 18, 4, dimnames = list(NULL, c("a",
"c", "g", "t")))
- ddd = fast.table(data)
- data = ddd$data
- index = ddd$index
+# ddd = fast.table(data)
+ compress=TRUE
+ if(length(data[[1]])==1) compress=FALSE
+ if(compress){
+ ddd = fast.table(data)
+ data = ddd$data
+ weight = ddd$weight
+ index = ddd$index
+ }
+ else{
+ p = length(data[[1]])
+ weight = rep(1, p)
+ index = 1:p
+ }
+# data = ddd$data
+# index = ddd$index
q = length(data)
p = length(data[[1]])
d = dim(data)
@@ -144,7 +176,8 @@ phyDat.DNA = function (data, return.index = TRUE)
rn = as.numeric(rownames(data))
attr(data, "na.action") = NULL
- weight = ddd$weight[rn]
+ weight = weight[rn]
+# weight = ddd$weight[rn]
p = dim(data)[1]
names(data) = nam
attr(data, "row.names") = NULL
@@ -168,6 +201,7 @@ phyDat.AA <- function (data, return.index = TRUE)
else nam = names(data)
if (inherits(data,"DNAbin"))
data = as.character(data)
+ if(inherits(data, "character")) data <- as.matrix(data)
if (is.matrix(data))
data = as.data.frame(t(data), stringsAsFactors = FALSE)
else data = as.data.frame(data, stringsAsFactors = FALSE)
@@ -186,9 +220,22 @@ phyDat.AA <- function (data, return.index = TRUE)
AA[23:25, ] = 1
dimnames(AA) <- list(aa2, aa)
- ddd = fast.table(data)
- data = ddd$data
- index = ddd$index
+# ddd = fast.table(data)
+# data = ddd$data
+# index = ddd$index
+ compress=TRUE
+ if(length(data[[1]])==1) compress=FALSE
+ if(compress){
+ ddd = fast.table(data)
+ data = ddd$data
+ weight = ddd$weight
+ index = ddd$index
+ }
+ else{
+ p = length(data[[1]])
+ weight = rep(1, p)
+ index = 1:p
+ }
q = length(data)
p = length(data[[1]])
tmp <- vector("list", q)
@@ -209,8 +256,9 @@ phyDat.AA <- function (data, return.index = TRUE)
index = match(index, unique(index))
rn = as.numeric(rownames(data))
attr(data, "na.action") = NULL
-
- weight = ddd$weight[rn]
+
+# weight = ddd$weight[rn]
+ weight = weight[rn]
p = dim(data)[1]
names(data) = nam
attr(data, "row.names") = NULL
@@ -235,7 +283,7 @@ phyDat.codon <- function (data, return.index = TRUE)
else nam = names(data)
if (inherits(data,"DNAbin"))
data = as.character(data)
-
+ if(inherits(data, "character")) data <- as.matrix(data)
if (is.matrix(data))
data = as.data.frame(t(data), stringsAsFactors = FALSE)
else data = as.data.frame(data, stringsAsFactors = FALSE)
@@ -250,9 +298,21 @@ phyDat.codon <- function (data, return.index = TRUE)
sapply(starts, function(x) paste(seq[x:(x + 2L)], collapse=""))
}
- data = sapply(data, splseq)
-
- ddd = fast.table(data)
+ data = data.frame(lapply(data, splseq))
+# ddd = fast.table(data)
+ compress=TRUE
+ if(nrow(data)==1) compress=FALSE
+ if(compress){
+ ddd = fast.table(data)
+ data = ddd$data
+ weight = ddd$weight
+ index = ddd$index
+ }
+ else{
+ p = length(data[[1]])
+ weight = rep(1, p)
+ index = 1:p
+ }
codon = c("aaa", "aac", "aag", "aat", "aca", "acc", "acg", "act",
"aga", "agc", "agg", "agt", "ata", "atc", "atg", "att",
"caa", "cac", "cag", "cat", "cca", "ccc", "ccg", "cct", "cga",
@@ -266,8 +326,8 @@ phyDat.codon <- function (data, return.index = TRUE)
CODON <- diag(61)
dimnames(CODON) <- list(codon, codon)
- data = ddd$data
- index = ddd$index
+# data = ddd$data
+# index = ddd$index
q = length(data)
p = length(data[[1]])
tmp <- vector("list", q)
@@ -289,7 +349,8 @@ phyDat.codon <- function (data, return.index = TRUE)
rn = as.numeric(rownames(data))
attr(data, "na.action") = NULL
- weight = ddd$weight[rn]
+# weight = ddd$weight[rn]
+ weight = weight[rn]
p = dim(data)[1]
names(data) = nam
attr(data, "row.names") = NULL
@@ -307,15 +368,33 @@ phyDat.codon <- function (data, return.index = TRUE)
}
+dna2codon <- function(x){
+ if(!inherits(x, "phyDat"))stop("x needs to be of class phyDat!")
+ phyDat.codon(as.character(x))
+}
+
+
+codon2dna <- function(x){
+ if(!inherits(x, "phyDat"))stop("x needs to be of class phyDat!")
+ phyDat.DNA(as.character(x))
+}
+
+
as.phyDat <- function (x, ...){
if (inherits(x,"phyDat")) return(x)
UseMethod("as.phyDat")
}
-as.phyDat.DNAbin <- function(x,...) phyDat.DNA(x,...)
-
+as.phyDat.factor <- function(x, ...){
+ nam <- names(x)
+ lev <- levels(x)
+ x <- as.character(x)
+ names(x) <- nam
+ phyDat(x, type="USER", levels = lev, ...)
+}
+as.phyDat.DNAbin <- function(x,...) phyDat.DNA(x,...)
as.phyDat.alignment <- function (x, type="DNA",...)
@@ -367,7 +446,7 @@ as.MultipleAlignment <- function (x, ...){
}
-as.MultipleAlignment.phyDat <- function(x){
+as.MultipleAlignment.phyDat <- function(x, ...){
if (requireNamespace('Biostrings')){
z = as.character(x)
nam = rownames(z)
@@ -388,6 +467,9 @@ phyDat2MultipleAlignment <- as.MultipleAlignment.phyDat
as.phyDat.matrix <- function (x, ...) phyDat(data=x, ...)
+as.phyDat.character <- function (x, ...) phyDat(data=x, ...)
+
+
as.phyDat.data.frame <- function (x, ...) phyDat(data=x, ...)
@@ -413,24 +495,23 @@ acgt2ry <- function(obj){
}
+# replace as.character.phyDat weniger Zeilen, works also for codons
as.character.phyDat <- function (x, allLevels=TRUE, ...)
{
nr <- attr(x, "nr")
nc <- attr(x, "nc")
type <- attr(x, "type")
- if (type == "DNA") {
- labels <- c("a", "c", "g", "t", "u", "m", "r", "w", "s",
- "y", "k", "v", "h", "d", "b", "n", "?", "-")
- }
- if (type == "AA") {
- labels <- c("a", "r", "n", "d", "c", "q", "e", "g", "h",
- "i", "l", "k", "m", "f", "p", "s", "t", "w", "y",
- "v", "b", "z", "x", "-", "?")
- }
+ labels = attr(x, "allLevels")
+
+ if (!is.null(attr(x, "index"))) {
+ index = attr(x, "index")
+ if (is.data.frame(index))
+ index <- index[, 1]
+ }
+ else index = rep(1:nr, attr(x, "weight"))
if (type == "USER") {
- #levels
- if(allLevels)labels = attr(x, "allLevels")
- else{
+ #levels in acgt2ry
+ if(!allLevels){
tmp = attr(x, "levels")
contrast = attr(x, "contrast") # contrast=AC
contrast[contrast>0] = 1
@@ -438,70 +519,30 @@ as.character.phyDat <- function (x, allLevels=TRUE, ...)
contrast[rowSums(contrast)>1,] = 0
labels = rep(NA, length(attr(x, "allLevels")))
labels[ind] = tmp[contrast%*%c(1:length(tmp))]
- }
+ }
}
- result = matrix(NA, nrow = length(x), ncol = nr)
- for (i in 1:length(x)) result[i, ] <- labels[x[[i]]]
- if (is.null(attr(x, "index")))
- index = rep(1:nr, attr(x, "weight"))
- else {
- index = attr(x, "index")
- if (is.data.frame(index))
- index <- index[, 1]
+ if(type == "AA") labels <- toupper(labels)
+ if(type == "CODON"){
+ nr <- length(index)
+ result = matrix(NA, nrow = length(x), ncol = 3L*nr)
+ labels = strsplit(labels, "")
+ for (i in 1:length(x)) result[i, ] <- unlist(labels[ x[[i]][index] ])
}
- result = result[, index, drop = FALSE]
- rownames(result) = names(x)
- result
-}
-
-
-# replace as.character.phyDat 20 Zeilen weniger
-as.character.phyDat2 <- function (x, ...)
-{
- nr <- attr(x, "nr")
- nc <- attr(x, "nc")
- type <- attr(x, "type")
- labels = attr(x, "allLevels")
- result = matrix(NA, nrow = length(x), ncol = nr)
- for (i in 1:length(x)) result[i, ] <- labels[x[[i]]]
- if (is.null(attr(x, "index")))
- index = rep(1:nr, attr(x, "weight"))
else {
- index = attr(x, "index")
- if (is.data.frame(index))
- index <- index[, 1]
+ result = matrix(NA, nrow = length(x), ncol = nr)
+ for (i in 1:length(x)) result[i, ] <- labels[x[[i]]]
+ result = result[, index, drop = FALSE]
}
- result = result[, index, drop = FALSE]
rownames(result) = names(x)
result
}
-#as.data.frame.phyDat <- function(x, ...){
-# data.frame(t(as.character(x, ...)), stringsAsFactors=FALSE)
-#}
-
-# much faster
-# TODO as stringsAsFactors=FALSE
-# result[[i]] <- x[[i]] + factor levels setzen
-#
-as.data.frame.phyDatOld <- function(x, ...){
- nr <- attr(x, "nr")
- nc <- attr(x, "nc")
- labels <- attr(x, "allLevels")
- result <- vector("list", length(x))
- for (i in 1:length(x)) result[[i]] <- labels[x[[i]]]
- attr(result, "names") <- names(x)
- attr(result, "row.names") <- 1:nr
- attr(result, "class") <- "data.frame"
- result
-}
-
-
as.data.frame.phyDat <- function(x, ...){
nr <- attr(x, "nr")
nc <- attr(x, "nc")
labels <- attr(x, "allLevels")
+ if(attr(x, "type") == "AA") labels <- toupper(labels)
result <- vector("list", length(x))
if (is.null(attr(x, "index")))
index = rep(1:nr, attr(x, "weight"))
@@ -560,7 +601,7 @@ as.DNAbin.phyDat <- function (x, ...)
}
-if (getRversion() >= "2.15.1") utils::globalVariables("as.AAbin")
+#if (getRversion() >= "2.15.1") utils::globalVariables("as.AAbin")
as.AAbin.phyDat <- function(x,...) {
if(attr(x, "type")=="AA") return(as.AAbin(as.character(x, ...)))
else stop("x must be a amino acid sequence")
@@ -586,51 +627,42 @@ print.phyDat = function (x, ...)
}
-c.phyDat <- function(...){
- object <- as.list(substitute(list(...)))[-1]
- x <- list(...)
- n <- length(x)
- match.names <- function(a,b){
- if(any(!(a %in% b)))stop("names do not match previous names")
- }
- if (n == 1)
- return(x[[1]])
- type <- attr(x[[1]], "type")
- nr = numeric(n)
- nr[1] <- sum(attr(x[[1]], "weight"))
- levels <- attr(x[[1]], "levels")
- snames <- names(x[[1]])
- objNames<-as.character(object)
- if(any(duplicated(objNames))) objNames <- paste0(objNames, 1:n)
- tmp <- as.character(x[[1]])
- for(i in 2:n){
- match.names(snames,names(x[[i]]))
- x[[i]] <- getCols(x[[i]],snames)
- nr[i] <- sum(attr(x[[i]], "weight"))
- tmp <- cbind(tmp, as.character(x[[i]]))
- }
- if (type == "DNA")
- dat <- phyDat.DNA(tmp, return.index = TRUE)
- if (type == "AA")
- dat <- phyDat.AA(tmp, return.index = TRUE)
- if (type == "USER")
- dat <- phyDat.default(tmp, levels = levels, return.index = TRUE)
- if (type == "CODON")
- dat <- phyDat.codon(tmp, return.index = TRUE)
- attr(dat,"index") <- data.frame(index=attr(dat,"index"), genes=rep(objNames, nr))
- dat
+# in C++ to replace aggregate or use of distinct from dplyr / data.table
+aggr <- function(weight, ind){
+ res <- numeric(max(ind))
+ for(i in seq_along(weight))
+ res[ind[i]] = res[ind[i]] + weight[i]
+ res
}
+# data has to be a data.frame in cbind.phyDat
+fast.table2 <- function (data)
+{
+ if(!is.data.frame(data))
+ data <- as.data.frame(data, stringsAsFactors = FALSE)
+ da <- do.call("paste", data)
+ ind <- !duplicated(da)
+ levels <- da[ind]
+ cat <- factor(da,levels = levels)
+ nl <- length(levels)
+ bin <- (as.integer(cat) - 1L)
+ bin <- bin[!is.na(bin)]
+ if (length(bin)) bin <- bin + 1L
+ result=list(index = bin, pos = ind)
+ result
+}
+
+
# new cbind.phyDat
-cbindPD <- function(..., gaps="-"){
+cbind.phyDat <- function(..., gaps="-", compress=TRUE){
object <- as.list(substitute(list(...)))[-1]
x <- list(...)
n <- length(x)
if (n == 1)
return(x[[1]])
type <- attr(x[[1]], "type")
- nr = numeric(n)
+ nr <- numeric(n)
ATTR <- attributes(x[[1]])
@@ -639,96 +671,73 @@ cbindPD <- function(..., gaps="-"){
allLevels <- attr(x[[1]], "allLevels")
gapsInd <- match(gaps, allLevels)
snames <- vector("list", n) # names(x[[1]])
- vec = numeric(n+1)
- wvec = numeric(n+1)
- objNames<-as.character(object)
+ vec <- numeric(n+1)
+ wvec <- numeric(n+1)
+ objNames <- as.character(object)
if(any(duplicated(objNames))) objNames <- paste0(objNames, 1:n)
# tmp <- as.character(x[[1]])
for(i in 1:n){
- snames[[i]] = names(x[[i]])
- nr[i] <- attr(x[[i]], "nr")
- vec[i+1] = attr(x[[i]], "nr")
- wvec[i+1] = sum(attr(x[[i]], "weight"))
+ snames[[i]] <- names(x[[i]])
+ nr[i] <- sum(attr(x[[i]], "weight"))
+ vec[i+1] <- attr(x[[i]], "nr")
+ wvec[i+1] <- sum(attr(x[[i]], "weight"))
}
- vec = cumsum(vec)
- wvec = cumsum(wvec)
- snames = unique(unlist(snames))
+ vec <- cumsum(vec)
+ wvec <- cumsum(wvec)
+ snames <- unique(unlist(snames))
weight <- numeric(vec[n+1])
- index <- numeric(wvec[n+1])
+ index <- numeric(wvec[n+1])
+
ATTR$names <- snames
ATTR$nr <- vec[n+1]
- tmp = matrix(gapsInd, vec[n+1], length(snames), dimnames = list(NULL, snames))
+ tmp <- matrix(gapsInd, vec[n+1], length(snames), dimnames = list(NULL, snames))
tmp <- as.data.frame(tmp)
-
+ add.index <- TRUE
for(i in 1:n){
- nam = names(x[[i]])
+ nam <- names(x[[i]])
tmp[(vec[i]+1):vec[i+1], nam] <- x[[i]][nam]
weight[(vec[i]+1):vec[i+1]] <- attr(x[[i]], "weight")
- index[(wvec[i]+1):wvec[i+1]] <- attr(x[[i]], "index")
}
- ATTR$index <- index
+ if(compress){
+ ddd <- fast.table2(tmp)
+ tmp <- tmp[ddd$pos,]
+ weight <- aggregate(weight, by=list(ddd$index), FUN=sum)$x
+ }
+ for(i in 1:n){
+ tmp2 <- attr(x[[i]], "index")
+ if(!is.null(tmp2)){
+ if(is.data.frame(tmp2))index[(wvec[i]+1):wvec[i+1]] <- ddd$index[(vec[i]+1):vec[i+1]][tmp2[,1]]
+ else index[(wvec[i]+1):wvec[i+1]] <- ddd$index[(vec[i]+1):vec[i+1]][tmp2]
+ }
+ else add.index <- FALSE
+ }
+ if(add.index)ATTR$index <- data.frame(index = index, genes=rep(objNames, nr))
ATTR$weight <- weight
+ ATTR$nr <- length(weight)
attributes(tmp) <- ATTR
tmp
}
+c.phyDat <- cbind.phyDat
-cbind.phyDat <- function(..., gaps="-"){
- object <- as.list(substitute(list(...)))[-1]
- x <- list(...)
- n <- length(x)
- if (n == 1)
- return(x[[1]])
- type <- attr(x[[1]], "type")
- nr = numeric(n)
- nr[1] <- sum(attr(x[[1]], "weight"))
- levels <- attr(x[[1]], "levels")
- snames <- vector("list", n) # names(x[[1]])
- vec = numeric(n+1)
- objNames<-as.character(object)
- if(any(duplicated(objNames))) objNames <- paste0(objNames, 1:n)
- tmp <- as.character(x[[1]])
- for(i in 1:n){
- snames[[i]] = names(x[[i]]) #match.names(snames,names(x[[i]]))
- nr[i] <- sum(attr(x[[i]], "weight"))
- vec[i+1] = sum(attr(x[[i]], "weight"))
- }
- vec = cumsum(vec)
- snames = unique(unlist(snames))
- tmp = matrix(gaps, length(snames), vec[n+1], dimnames = list(snames, NULL))
-
- for(i in 1:n){
- nam = names(x[[i]])
- tmp[nam,(vec[i]+1):vec[i+1] ] <- as.character(x[[i]])
- }
- if (type == "DNA")
- dat <- phyDat.DNA(tmp, return.index = TRUE)
- if (type == "AA")
- dat <- phyDat.AA(tmp, return.index = TRUE)
- if (type == "USER")
- dat <- phyDat.default(tmp, levels = levels,
- return.index = TRUE)
- if (type == "CODON")
- dat <- phyDat.codon(tmp, return.index = TRUE)
- attr(dat,"index") <- data.frame(index=attr(dat,"index"), genes=rep(objNames, nr))
- dat
-}
-
-
-write.phyDat <- function(x, file, format="phylip",...){
- if(format=="fasta") write.dna(as.character(x), file, format="fasta", colsep = "", nbcol=-1, ...)
- if(format=="phylip") write.dna(as.character(x), file, format="sequential", ...)
+write.phyDat <- function(x, file, format="phylip", colsep = "", nbcol=-1, ...){
+ formats <- c("phylip", "nexus", "interleaved", "sequential", "fasta")
+ format <- match.arg(tolower(format), formats)
if(format=="nexus"){
- type = attr(x, "type")
- if(type=="DNA") write.nexus.data(as.list(as.data.frame(x)), file, format = "dna",...)
- else write.nexus.data(as.list(as.data.frame(x)), file, format = "protein", ...)
- }
+ type = attr(x, "type")
+ if(type=="DNA") write.nexus.data(as.list(as.data.frame(x)), file, format = "dna",...)
+ else write.nexus.data(as.list(as.data.frame(x)), file, format = "protein", ...)
}
+ else{
+ if(format=="phylip") format = "interleaved"
+ write.dna(as.character(x), file, format=format, colsep = colsep, nbcol=nbcol, ...)
+ }
+}
read.phyDat <- function(file, format="phylip", type="DNA", ...){
@@ -1091,7 +1100,6 @@ genlight2phyDat <- function(x, ambiguity=NA){
}
-if (getRversion() >= "2.15.1") utils::globalVariables("image.AAbin")
image.phyDat <- function(x, ...){
if(attr(x, "type")=="AA")image(as.AAbin(x), ...)
if(attr(x, "type")=="DNA")image(as.DNAbin(x), ...)
diff --git a/R/phylo.R b/R/phylo.R
index 69cdcff..1d30607 100644
--- a/R/phylo.R
+++ b/R/phylo.R
@@ -732,7 +732,7 @@ likelihoodRatchet <- function(obj, maxit=100, k=10,
obj
}
-
+
fs <- function (old.el, eig, parent.dat, child.dat, weight, g=g,
w=w, bf=bf, ll.0=ll.0, evi, getA=TRUE, getB=TRUE)
{
@@ -746,8 +746,8 @@ fs <- function (old.el, eig, parent.dat, child.dat, weight, g=g,
X <- .Call("getPrep", dad, child.dat, eig[[2]], evi, nr, nc)
.Call("FS4", eig, as.integer(length(bf)), as.double(old.el),
as.double(w), as.double(g), X, child.dat, dad, as.integer(length(w)),
- as.integer(length(weight)), as.double(bf), as.double(weight),
- as.double(ll.0), as.integer(getA), as.integer(getB))
+ as.integer(length(weight)), as.double(weight),
+ as.double(ll.0), as.integer(getA), as.integer(getB)) # as.double(bf),
}
@@ -774,8 +774,8 @@ fs3 <- function (old.el, eig, parent.dat, child, weight, g=g,
}
.Call("FS4", eig, as.integer(length(bf)), as.double(old.el),
as.double(w), as.double(g), X, child.dat, dad, as.integer(length(w)),
- as.integer(length(weight)), as.double(bf), as.double(weight),
- as.double(ll.0), as.integer(getA), as.integer(getB))
+ as.integer(length(weight)), as.double(weight),
+ as.double(ll.0), as.integer(getA), as.integer(getB)) # as.double(bf),
}
@@ -852,7 +852,8 @@ pml.move <- function(EDGE, el, data, g, w, eig, k, nTips, bf){
edge = as.integer(edge - 1L)
contrast = attr(data, "contrast")
nco = as.integer(dim(contrast)[1])
- tmp <- .Call("PML3", dlist=data, as.double(el), as.double(w), as.double(g), nr, nc, k, eig, as.double(bf), node, edge, nTips, nco, contrast, N=as.integer(length(edge)))
+ tmp <- .Call("PML3", dlist=data, as.double(el), as.double(g), nr, nc, k, eig, as.double(bf), node, edge, nTips, nco, contrast, N=as.integer(length(edge)))
+# as.double(w),
return(NULL)
}
@@ -2680,7 +2681,8 @@ pml.fit <- function (tree, data, bf = rep(1/length(levels), length(levels)),
nco = as.integer(dim(contrast)[1])
# dlist=data, nr, nc, weight, k ausserhalb definieren
# pmlPart einbeziehen
- resll <- .Call("PML0", dlist=data, el, as.double(w), as.double(g), nr, nc, k, eig, as.double(bf), node, edge, nTips, nco, contrast, N=as.integer(length(edge)))
+ # as.double(w),
+ resll <- .Call("PML0", dlist=data, el, as.double(g), nr, nc, k, eig, as.double(bf), node, edge, nTips, nco, contrast, N=as.integer(length(edge)))
# sort(INV at i)+1L
ind = which(ll.0>0) # automatic in INV gespeichert
@@ -2707,6 +2709,7 @@ pml <- function (tree, data, bf = NULL, Q = NULL, inv = 0, k = 1, shape = 1,
rate = 1, model=NULL, ...)
{
Mkv = FALSE
+ if(!is.null(model) && model=="Mkv")Mkv = TRUE
call <- match.call()
extras <- match.call(expand.dots = FALSE)$...
pmla <- c("wMix", "llMix")
@@ -2719,6 +2722,7 @@ pml <- function (tree, data, bf = NULL, Q = NULL, inv = 0, k = 1, shape = 1,
# warning("tree has no edge length, used nnls.phylo to assign them")
# tree <- nnls.phylo(tree, dist.ml(data))
# }
+ if(any(duplicated(tree$tip.label))) stop("tree must have unique labels!")
if (is.null(attr(tree, "order")) || attr(tree, "order") ==
"cladewise")
tree <- reorder(tree, "postorder")
@@ -3902,7 +3906,8 @@ index2tree2 <- function(x, tree, root=length(tree$tip.label)+1L){
# weight, nr, nc, contrast, nco (Reihenfolge beibehalten)
# INV raus
-optimQuartet <- function (tree, data, eig, w, g, bf, rate, ll.0=ll.0, nTips,
+# evi, eve, contrast2 ausserhalb definieren
+optimQuartet <- function (tree, data, eig, w, g, bf, rate, ll.0, nTips,
weight, nr, nc, contrast, nco, llcomp =-Inf,
control = pml.control(epsilon = 1e-08, maxit = 5, trace=0), ...)
{
@@ -3925,7 +3930,7 @@ optimQuartet <- function (tree, data, eig, w, g, bf, rate, ll.0=ll.0, nTips,
eps = 1
iter = 0
- treeP = tree
+# treeP = tree
child = tree$edge[, 2]
parent = tree$edge[, 1]
@@ -3953,24 +3958,25 @@ optimQuartet <- function (tree, data, eig, w, g, bf, rate, ll.0=ll.0, nTips,
as.double(contrast2), nco, data, as.double(weight), as.double(ll.0))
iter = iter + 1
# tree$edge.length = EL[tree$edge[,2]]
- treeP$edge.length = EL # [treeP$edge[,2]]
- newll <- pml.quartet(treeP, data, bf=bf, g=g, w=w, eig=eig, ll.0=ll.0, k=k, nTips=nTips,
- weight=weight, nr=nr, nc=nc, contrast=contrast, nco=nco)
+ tree$edge.length = EL # [treeP$edge[,2]] # vormals treeP
+ newll <- pml.quartet(tree, data, bf=bf, g=g, w=w, eig=eig, ll.0=ll.0, k=k, nTips=nTips,
+ weight=weight, nr=nr, nc=nc, contrast=contrast, nco=nco) # vormals treeP
eps = ( old.ll - newll ) / newll
if( (eps<0) || (newll < llcomp) ) return(list(tree=oldtree, logLik=old.ll, c(eps, iter)))
- oldtree = treeP
+ oldtree = tree # vormals treeP
if(control$trace>1) cat(old.ll, " -> ", newll, "\n")
old.ll = newll
# loli = parent[1]
}
if(control$trace>0) cat(start.ll, " -> ", newll, "\n")
- list(tree=treeP, logLik=newll, c(eps, iter))
+ list(tree=tree, logLik=newll, c(eps, iter)) # vormals treeP
}
# weight, nr, nc, contrast, nco rein
-# inv, INV raus
+# inv, INV, site, ... raus
+# wMix, rate last
pml.quartet <- function (tree, data, bf = rep(.25, 4), k = 1, rate = 1, g, w,
eig, ll.0 = NULL, ind.ll0=NULL, llMix = NULL, wMix = 0, nTips,
weight, nr, nc, contrast, nco, ..., site=FALSE)
diff --git a/R/sankoff.R b/R/sankoff.R
index 53baec7..768debf 100644
--- a/R/sankoff.R
+++ b/R/sankoff.R
@@ -27,8 +27,6 @@ sankoffNew <- function (tree, data, cost = NULL, site = 'pscore')
# on.exit(.C("sankoff_free"))
# .C("sankoff_init", as.integer())
-
-
# for (i in 1:length(data)) storage.mode(data[[i]]) = "double"
if(inherits(tree,"phylo")) return(fit.sankoffNew(tree, data, cost, returnData =site))
diff --git a/R/simSeq.R b/R/simSeq.R
index 3eff2e1..8c9cb71 100644
--- a/R/simSeq.R
+++ b/R/simSeq.R
@@ -9,12 +9,23 @@ simSeq <- function (x, ...)
simSeq.phylo = function(x, l=1000, Q=NULL, bf=NULL, rootseq=NULL, type = "DNA", model="USER",
levels = NULL, rate=1, ancestral=FALSE, ...){
- pt <- match.arg(type, c("DNA", "AA", "USER"))
+ pt <- match.arg(type, c("DNA", "AA", "USER", "CODON"))
if (pt == "DNA")
levels <- c("a", "c", "g", "t")
if (pt == "AA")
levels <- c("a", "r", "n", "d", "c", "q", "e", "g", "h", "i",
"l", "k", "m", "f", "p", "s", "t", "w", "y", "v")
+ if (pt == "CODON"){
+ levels <- c("aaa", "aac", "aag", "aat", "aca", "acc", "acg", "act",
+ "aga", "agc", "agg", "agt", "ata", "atc", "atg", "att",
+ "caa", "cac", "cag", "cat", "cca", "ccc", "ccg", "cct", "cga",
+ "cgc", "cgg", "cgt", "cta", "ctc", "ctg", "ctt", "gaa", "gac",
+ "gag", "gat", "gca", "gcc", "gcg", "gct", "gga", "ggc", "ggg",
+ "ggt", "gta", "gtc", "gtg", "gtt", "tac", "tat",
+ "tca", "tcc", "tcg", "tct", "tgc", "tgg", "tgt", "tta",
+ "ttc", "ttg", "ttt")
+ Q <- as.numeric(.syn > 0)
+ }
if (pt == "USER")
if(is.null(levels))stop("levels have to be supplied if type is USER")
@@ -59,6 +70,10 @@ simSeq.phylo = function(x, l=1000, Q=NULL, bf=NULL, rootseq=NULL, type = "DNA",
if(pt=="DNA") return(phyDat.DNA(as.data.frame(res, stringsAsFactors = FALSE), return.index=TRUE))
if(pt=="AA") return(phyDat.AA(as.data.frame(res, stringsAsFactors = FALSE), return.index=TRUE))
if(pt=="USER") return(phyDat.default(as.data.frame(res, stringsAsFactors = FALSE), levels = levels, return.index=TRUE))
+ if(pt=="CODON"){
+ res <- apply(res, 2, function(x)unlist(strsplit(x, "")))
+ return(phyDat.codon(as.data.frame(res, stringsAsFactors = FALSE)))
+ }
}
diff --git a/R/superTree.R b/R/superTree.R
new file mode 100644
index 0000000..90b97f6
--- /dev/null
+++ b/R/superTree.R
@@ -0,0 +1,135 @@
+# now more memoryefficient
+# from phytools code by Liam Revell with a few changes
+my.supertree<-function(trees, trace=0, ...){
+ # some minor error checking
+ if(!inherits(trees,"multiPhylo")) stop("trees must be object of class 'multiPhylo.'")
+
+ labels <- lapply(trees, function(x)sort(x$tip.label))
+ ulabels <- unique(labels)
+ lul <- length(ulabels)
+ # compute matrix representation phylogenies
+ X<-vector("list", lul) # list of bipartitions
+ characters<-0 # number of characters
+ weights <- NULL
+ species<-trees[[1]]$tip.label
+ for(i in 1:lul){
+ pos <- match(labels, ulabels[i])
+ ind <- which(!is.na(pos))
+ temp<-prop.part(trees[ind]) # find all bipartitions
+ # create matrix representation of trees[[i]] in X[[i]]
+ X[[i]]<-matrix(0,nrow=length(trees[[ind[1]]]$tip.label),ncol=length(temp)-1)
+ for(j in 1:ncol(X[[i]])) X[[i]][c(temp[[j+1]]),j]<-1
+ rownames(X[[i]])<-attr(temp,"labels") # label rows
+ # if(i==1) species<-trees[[ind[1]]]$tip.label
+ # else
+ species<-union(species,trees[[ind[1]]]$tip.label) # accumulate labels
+ characters<-characters+ncol(X[[i]]) # count characters
+ weights <- c(weights, attr(temp, "number")[-1])
+ }
+ XX<-matrix(data="?",nrow=length(species),ncol=characters,dimnames=list(species))
+ j<-1
+ for(i in 1:length(X)){
+ # copy each of X into supermatrix XX
+ XX[rownames(X[[i]]),c(j:((j-1)+ncol(X[[i]])))]<-X[[i]][1:nrow(X[[i]]),1:ncol(X[[i]])]
+ j<-j+ncol(X[[i]])
+ }
+ # compute contrast matrix for phangorn
+ contrast<-matrix(data=c(1,0,0,1,1,1),3,2,dimnames=list(c("0","1","?"),c("0","1")),byrow=TRUE)
+ # convert XX to phyDat object
+ XX<-phyDat(XX,type="USER",contrast=contrast, compress=FALSE)
+ attr(XX, "weight") <- weights
+ # estimate supertree
+ supertree<-pratchet(XX,all=TRUE, trace=trace, ...)
+ return(supertree)
+}
+
+
+# Robinson-Foulds supertree
+fun.rf <- function(x, tree) sum(RF.dist(x, tree))
+fun.spr <- function(x, tree) sum(SPR.dist(x, tree))
+
+dist.superTree <- function(tree, trace=0, fun, start=NULL, multicore=FALSE, mc.cores = NULL){
+ if(multicore && is.null(mc.cores)){
+ mc.cores <- detectCores()
+ }
+ if(is.null(start))start <- superTree(tree, rooted=FALSE)
+ if(inherits(start, "multiPhylo")) start <- start[[1]]
+ best_tree <- unroot(start)
+ best <- fun(best_tree, tree)
+ if(trace>0)cat("best score so far:", best, "\n")
+ eps <- TRUE
+ while(eps){
+ nni_trees <- nni(best_tree)
+ if (multicore) {
+ tmp <- mclapply(nni_trees, fun, tree, mc.cores = mc.cores)
+ tmp <-unlist(tmp)
+ }
+ else tmp <- sapply(nni_trees, fun, tree)
+ if(min(tmp) < best){
+ ind <- which.min(tmp)
+ best_tree <- nni_trees[[ind]]
+ best <- tmp[ind]
+ if(trace>0)cat("best score so far:", best, "\n")
+ }
+ else eps <- FALSE
+ }
+ attr(best_tree, "score") <- best
+ best_tree
+}
+
+
+# we want a rooted supertree
+superTree = function(tree, method="MRP", rooted=FALSE, trace=0, ...){
+ fun = function(x){
+ x=reorder(x, "postorder")
+ nTips = length(x$tip.label)
+ x$edge[x$edge>nTips] = x$edge[x$edge>nTips] + 2L
+ l=nrow(x$edge)
+ oldroot = x$edge[l,1L]
+ x$edge=rbind(x$edge,matrix(c(rep(nTips+2,2),oldroot,nTips+1),2L,2L))
+ x$edge.length=c(x$edge.length, 100, 100)
+ x$tip.label=c(x$tip.label, "ZZZ")
+ x$Nnode=x$Nnode+1L
+ x
+ }
+ if(!is.null(attr(tree, "TipLabel")))tree = .uncompressTipLabel(tree)
+ tree = unclass(tree)
+ if(rooted) tree = lapply(tree, fun)
+ class(tree)="multiPhylo"
+ res = my.supertree(tree, trace=trace, ...)
+ if(rooted){
+ if(inherits(res,"multiPhylo")){
+ res = lapply(res, root, "ZZZ")
+ res = lapply(res, drop.tip, "ZZZ")
+ class(res) = "multiPhylo"
+ }
+ else{
+ res = root(res, "ZZZ")
+ res = drop.tip(res, "ZZZ")
+ }
+ }
+ if(inherits(res,"multiPhylo")){
+ fun = function(x){
+ x$edge.length <- rep(.1, nrow(x$edge))
+ x
+ }
+ res <- lapply(res, fun)
+ res <- lapply(res, reorder, "postorder")
+ class(res) = "multiPhylo"
+ }
+ else{
+ res$edge.length = rep(.1, nrow(res$edge))
+ res <- reorder(res, "postorder")
+ }
+ if(method=="MRP") return(res)
+
+ tree <- unroot(tree)
+ tree <- reorder(tree, "postorder")
+# tree <- lapply(tree, unroot)
+# tree <- lapply(tree, reorder, "postorder")
+# class(tree) = "multiPhylo"
+
+ if(method=="RF") res <- dist.superTree(tree, trace=trace, fun.rf, start=res)
+ if(method=="SPR") res <- dist.superTree(tree, trace=trace, fun.spr, start=res)
+ res
+}
diff --git a/R/treeManipulation.R b/R/treeManipulation.R
index 8f4e9f2..18cab99 100644
--- a/R/treeManipulation.R
+++ b/R/treeManipulation.R
@@ -551,7 +551,8 @@ kSPR = function(tree, k=NULL){
distN = dn(tree)[-c(1:l), -c(1:l)]
distN[upper.tri(distN)]=Inf
dN = distN[lower.tri(distN)]
- tab = table(dN)
+ tab = tabulate(dN) # should be much faster
+# tab = table(dN)
tab[1] = tab[1] * 2
tab[-1] = tab[-1] * 8
if(is.null(k)) k = 1:length(tab)
@@ -774,7 +775,7 @@ Ancestors <- function (x, node, type = c("all", "parent"))
{
parents <- x$edge[, 1]
child <- x$edge[, 2]
- pvector <- numeric(max(x$edge)) # parents
+ pvector <- integer(max(x$edge)) # parents
pvector[child] <- parents
type <- match.arg(type)
if (type == "parent")
@@ -794,23 +795,49 @@ Ancestors <- function (x, node, type = c("all", "parent"))
}
+#allChildren <- function(x){
+# l = length(x$tip.label)
+# if(l<20){
+# parent = x$edge[,1]
+# children = x$edge[,2]
+# res = vector("list", max(x$edge))
+# for(i in 1:length(parent)) res[[parent[i]]] = c(res[[parent[i]]], children[i])
+# return(res)
+# }
+# else{
+# if (is.null(attr(x, "order")) || attr(x, "order") == "cladewise")
+# x <- reorder(x, "postorder")
+# parent = x$edge[,1]
+# children = x$edge[,2]
+# res <- .Call("AllChildren", as.integer(children), as.integer(parent), as.integer(max(x$edge))) # , PACKAGE="phangorn"
+# return(res)
+# }
+#}
+
+
+# alternative version using Rcpp
allChildren <- function(x){
- l = length(x$tip.label)
- if(l<20){
- parent = x$edge[,1]
- children = x$edge[,2]
- res = vector("list", max(x$edge))
- for(i in 1:length(parent)) res[[parent[i]]] = c(res[[parent[i]]], children[i])
- return(res)
- }
- else{
- if (is.null(attr(x, "order")) || attr(x, "order") == "cladewise")
- x <- reorder(x, "postorder")
- parent = x$edge[,1]
- children = x$edge[,2]
- res <- .Call("AllChildren", as.integer(children), as.integer(parent), as.integer(max(x$edge))) # , PACKAGE="phangorn"
- return(res)
- }
+ l = length(x$tip.label)
+ if(l<20){
+ parent = x$edge[,1]
+ children = x$edge[,2]
+ res = vector("list", max(x$edge))
+ for(i in 1:length(parent)) res[[parent[i]]] = c(res[[parent[i]]], children[i])
+ return(res)
+ }
+ else{
+ if (is.null(attr(x, "order")) || attr(x, "order") == "cladewise")
+ x <- reorder(x, "postorder")
+ allChildrenCPP(x$edge)
+ }
+}
+
+
+allDescendants <- function(x){
+ if (is.null(attr(x, "order")) || attr(x, "order") == "cladewise")
+ x <- reorder(x, "postorder")
+ nTips <- as.integer(length(x$tip.label))
+ allDescCPP(x$edge, nTips)
}
@@ -824,6 +851,8 @@ Descendants = function(x, node, type=c("tips","children","all")){
type <- match.arg(type)
if(type=="children") return(Children(x, node))
if(type=="tips") return(bip(x)[node])
+ # new version using Rcpp
+ if(length(node)>10) return(allDescendants(x)[node])
ch = allChildren(x) # out of the loop
isInternal = logical(max(x$edge))
isInternal[ unique(x$edge[,1]) ] =TRUE
diff --git a/R/treedist.R b/R/treedist.R
index 20bd0dd..772630b 100644
--- a/R/treedist.R
+++ b/R/treedist.R
@@ -56,7 +56,8 @@ cophenetic.networx <- function(x){
SHORTwise <- function (x, nTips, delete=FALSE)
{
v <- 1:nTips
- l <- sapply(x, length)
+# l <- sapply(x, length)
+ l <- lengths(x)
lv = floor(nTips/2)
for (i in 1:length(x)) {
if(l[i]>lv){
@@ -83,12 +84,121 @@ oneWise <- function (x, nTips=NULL)
for (i in 2:length(x)) {
y <- x[[i]]
if (y[1] != 1)
- x[[i]] <- v[-y]
+ y <- v[-y]
+ if (y[1] != 1)
+ y <- v[-y]
+ x[[i]] <- y
}
x
}
+# leomrtns addition
+sprdist <- function (tree1, tree2)
+{
+ tree1 = unroot(tree1)
+ tree2 = unroot(tree2)
+ lt1 = length(tree1$tip.label)
+ lt2 = length(tree2$tip.label)
+ # checking labels is obligatory for spr (user should prune one of them beforehand?)
+ ind <- match(tree1$tip.label, tree2$tip.label)
+ if (any(is.na(ind)) | lt1 != lt2 ) stop("trees have different labels")
+ tree2$tip.label <- tree2$tip.label[ind]
+ ind2 <- match(1:length(ind), tree2$edge[, 2])
+ tree2$edge[ind2, 2] <- order(ind)
+ # same as in original treedist (will create list of strings with shorter side of splits)
+ tree1 = reorder(tree1, "postorder")
+ tree2 = reorder(tree2, "postorder")
+ if(!is.binary.tree(tree1) | !is.binary.tree(tree2))message("Trees are not binary!")
+ # possibly replace bip with bipart
+ bp1 = bip(tree1); bp1 <- SHORTwise(bp1, lt1)
+ bp2 = bip(tree2); bp2 <- SHORTwise(bp2, lt2)
+
+ bp1 <- bp1[ lengths(bp1)>1 ] # only internal nodes
+ bp2 <- bp2[ lengths(bp2)>1 ]
+ if (length(bp1) != length(bp2)) stop ("number of bipartitions given to C_sprdist are not the same")
+ # OBS: SPR distance works w/ incompatible splits only, but it needs common cherries (to replace by single leaf)
+ spr <- .Call("C_sprdist", bp1, bp2, lt1);
+ names(spr) <- c("spr", "spr_extra", "rf", "hdist")
+# list("spr" = spr[1], "spr_extra" = spr[2], "rf" = spr[3], "hdist"= spr[4]);
+ spr
+}
+
+
+SPR1 <- function(trees){
+ trees <- .compressTipLabel(trees)
+ trees <- .uncompressTipLabel(trees)
+ trees <- lapply(trees, unroot)
+ trees <- lapply(trees, reorder, "postorder")
+
+ nTips <- length(trees[[1]]$tip.label)
+
+ fun <- function(x, nTips){
+ bp <- bipart(x)
+ bp <- SHORTwise(bp, nTips)
+ bp <- bp[ lengths(bp)>1 ]
+ bp
+ }
+
+ BP <- lapply(trees, fun, nTips)
+ k <- 1
+ l <- length(trees)
+ SPR <- numeric((l * (l - 1))/2)
+ for (i in 1:(l - 1)){
+ bp <- BP[[i]]
+ for (j in (i + 1):l){
+ SPR[k] <- .Call("C_sprdist", bp, BP[[j]], nTips)[1]
+ k=k+1
+ }
+ }
+ attr(SPR, "Size") <- l
+ if(!is.null(names(trees)))attr(SPR, "Labels") <- names(trees)
+ attr(SPR, "Diag") <- FALSE
+ attr(SPR, "Upper") <- FALSE
+ class(SPR) <- "dist"
+ return(SPR)
+}
+
+
+SPR2 <- function(tree, trees){
+ trees <- .compressTipLabel(trees)
+ tree <- checkLabels(tree, attr(trees, "TipLabel"))
+ trees <- .uncompressTipLabel(trees)
+ if (any(sapply(trees, is.rooted))) {
+ trees <- lapply(trees, unroot)
+ }
+ trees <- lapply(trees, reorder, "postorder")
+ tree <- unroot(tree)
+ nTips <- length(tree$tip.label)
+
+ fun <- function(x, nTips){
+ bp <- bipart(x)
+ bp <- SHORTwise(bp, nTips)
+ bp <- bp[ lengths(bp)>1 ]
+ bp
+ }
+
+ bp <- fun(tree, nTips)
+ k <- 1
+ l <- length(trees)
+ SPR <- numeric(l)
+ for (i in 1:l){
+ SPR[i] <- .Call("C_sprdist", bp, fun(trees[[i]], nTips), nTips)[1]
+ }
+ if(!is.null(names(trees)))names(SPR) <- names(trees)
+ return(SPR)
+}
+
+
+SPR.dist <- function(tree1, tree2=NULL){
+ if(inherits(tree1, "multiPhylo") && is.null(tree2))return(SPR1(tree1))
+ if(inherits(tree1, "phylo") && inherits(tree2, "phylo"))return(sprdist(tree1, tree2)[1])
+ if(inherits(tree1, "phylo") && inherits(tree2, "multiPhylo"))return(SPR2(tree1, tree2))
+ if(inherits(tree2, "phylo") && inherits(tree1, "multiPhylo"))return(SPR2(tree2, tree1))
+ return(NULL)
+}
+
+
treedist <- function (tree1, tree2, check.labels=TRUE)
{
tree1 = unroot(tree1)
@@ -324,7 +434,7 @@ wRF1 <- function(trees, normalize=FALSE, check.labels = TRUE){
}
}
attr(wRF, "Size") <- l
- if(!is.null(names(trees)))attr(KF, "Labels") <- names(trees)
+ if(!is.null(names(trees)))attr(wRF, "Labels") <- names(trees)
attr(wRF, "Diag") <- FALSE
attr(wRF, "Upper") <- FALSE
class(wRF) <- "dist"
@@ -359,6 +469,7 @@ mRF2 <- function(tree, trees, normalize=FALSE, check.labels = TRUE){
message("Some trees are rooted. Unrooting all trees.\n")
trees <- lapply(trees, unroot)
}
+ if(is.rooted(tree)) tree <- unroot(tree)
if (any(sapply(trees, function(x) !is.binary.tree(x)))) {
message("Some trees are not binary. Result may not what you expect!")
}
@@ -433,14 +544,24 @@ mRF<-function(trees, normalize=FALSE){
RF0 <- function(tree1, tree2=NULL, normalize=FALSE, check.labels = TRUE, rooted=FALSE){
r1 = is.rooted(tree1)
r2 = is.rooted(tree2)
- if(r1 != r2){
- message("one tree is unrooted, unrooted both")
- }
if(!rooted){
- if(r1) tree1<-unroot(tree1)
- if(r2) tree2<-unroot(tree2)
+ if(r1) {
+ tree1<-unroot(tree1)
+ r1 <- FALSE
+ }
+ if(r2) {
+ tree2<-unroot(tree2)
+ r2 <- FALSE
+ }
}
-
+ else{
+ if(r1 != r2) {
+ message("one tree is unrooted, unrooted both")
+ tree1<-unroot(tree1)
+ tree2<-unroot(tree2)
+ r1 <- r2 <- FALSE
+ }
+ }
if (check.labels) {
ind <- match(tree1$tip.label, tree2$tip.label)
if (any(is.na(ind)) | length(tree1$tip.label) !=
@@ -451,17 +572,13 @@ RF0 <- function(tree1, tree2=NULL, normalize=FALSE, check.labels = TRUE, rooted=
ind2 <- match(1:length(ind), tree2$edge[, 2])
tree2$edge[ind2, 2] <- order(ind)
}
-
- if(!r1 | !r2){
- if(r1) tree1 = unroot(tree1)
- if(r2) tree2 = unroot(tree2)
- }
if(!is.binary.tree(tree1) | !is.binary.tree(tree2))message("Trees are not binary!")
bp1 = bipart(tree1)
bp2 = bipart(tree2)
+ nTips <- length(tree1$tip.label)
if(!rooted){
- bp1 <- SHORTwise(bp1, length(tree1$tip.label))
- bp2 <- SHORTwise(bp2, length(tree2$tip.label))
+ bp1 <- SHORTwise(bp1, nTips)
+ bp2 <- SHORTwise(bp2, nTips)
}
RF = sum(match(bp1, bp2, nomatch=0L)==0L) + sum(match(bp2, bp1, nomatch=0L)==0L)
if(normalize) RF <- RF / (Nnode(tree1) + Nnode(tree2) - 2)
@@ -679,7 +796,10 @@ pd1 <- function(tree, trees, check.labels=TRUE, path=TRUE){
}
trees <- .uncompressTipLabel(trees)
unclass(trees)
- if(path)trees <- lapply(trees, unroot)
+ if(path){
+ trees <- lapply(trees, unroot)
+ tree <- unroot(tree)
+ }
trees <- lapply(trees, reorder, "postorder")
l <- length(trees)
dt <- coph(tree, path)
diff --git a/R/zzz.R b/R/zzz.R
index 12b83b2..0f6750b 100644
--- a/R/zzz.R
+++ b/R/zzz.R
@@ -4,6 +4,11 @@
.aamodels <- c("WAG", "JTT", "LG", "Dayhoff", "cpREV", "mtmam", "mtArt", "MtZoa", "mtREV24", "VT","RtREV", "HIVw", "HIVb", "FLU","Blosum62","Dayhoff_DCMut","JTT_DCMut")
+.dnamodels <- c("JC", "F81", "K80", "HKY", "TrNe", "TrN",
+ "TPM1", "K81", "TPM1u", "TPM2", "TPM2u", "TPM3", "TPM3u",
+ "TIM1e", "TIM1", "TIM2e", "TIM2", "TIM3e", "TIM3", "TVMe",
+ "TVM", "SYM", "GTR")
+
# if g[i] is smaller .gEps inv is increased w[i]
.gEps <- 1e-30
diff --git a/README.md b/README.md
index 9cf418f..f532211 100644
--- a/README.md
+++ b/README.md
@@ -1,7 +1,7 @@
[](https://travis-ci.org/KlausVigo/phangorn)
-[](http://cran.r-project.org/package=phangorn)
+[](https://cran.r-project.org/package=phangorn)
[](http://www.gnu.org/licenses/gpl-2.0.html)
-[](http://cran.r-project.org/package=phangorn)
+[](https://cran.r-project.org/package=phangorn)
[](http://depsy.org/package/r/phangorn)
[](https://codecov.io/github/KlausVigo/phangorn?branch=master)
@@ -17,7 +17,7 @@ You can install
You may need to install first the Biostrings package from bioconductor
```
-source("http://bioconductor.org/biocLite.R")
+source("https://bioconductor.org/biocLite.R")
biocLite("Biostrings")
```
diff --git a/build/vignette.rds b/build/vignette.rds
index 7f43017..5baba21 100644
Binary files a/build/vignette.rds and b/build/vignette.rds differ
diff --git a/inst/doc/Ancestral.pdf b/inst/doc/Ancestral.pdf
index f6050c7..ad558fc 100644
Binary files a/inst/doc/Ancestral.pdf and b/inst/doc/Ancestral.pdf differ
diff --git a/inst/doc/IntertwiningTreesAndNetworks.html b/inst/doc/IntertwiningTreesAndNetworks.html
index 4ab62d8..f273268 100644
--- a/inst/doc/IntertwiningTreesAndNetworks.html
+++ b/inst/doc/IntertwiningTreesAndNetworks.html
@@ -12,7 +12,7 @@
<meta name="author" content="Klaus Schliep, Alastair Potts, David Morrison and Guido Grimm" />
-<meta name="date" content="2016-06-20" />
+<meta name="date" content="2016-12-03" />
<title>Intertwining phylogenetic trees and networks: R Example Script</title>
@@ -70,7 +70,7 @@ code > span.in { color: #60a0b0; font-weight: bold; font-style: italic; } /* Inf
<h1 class="title toc-ignore">Intertwining phylogenetic trees and networks: R Example Script</h1>
<h4 class="author"><em>Klaus Schliep, Alastair Potts, David Morrison and Guido Grimm</em></h4>
-<h4 class="date"><em>June 20, 2016</em></h4>
+<h4 class="date"><em>December 03, 2016</em></h4>
@@ -222,7 +222,7 @@ y <-<span class="st"> </span><span class="kw">plot</span>(y,<span class="st">
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">cnet <-<span class="st"> </span><span class="kw">consensusNet</span>(raxml.bootstrap,<span class="dt">prob=</span><span class="fl">0.10</span>)
edge.col <-<span class="st"> </span><span class="kw">createLabel</span>(cnet, Nnet, <span class="dt">label=</span><span class="st">"black"</span>, <span class="st">"edge"</span>, <span class="dt">nomatch=</span><span class="st">"grey"</span>)
cnet <-<span class="st"> </span><span class="kw">plot.networx</span>(cnet, <span class="st">"2D"</span>, <span class="dt">show.edge.label =</span> <span class="ot">TRUE</span>, <span class="dt">edge.color=</span>edge.col)</code></pre></div>
-<p><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAqAAAAKgCAMAAABz4j/3AAADAFBMVEUAAAABAQECAgIDAwMEBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUWFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJycoKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6Ojo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tMTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1eXl5fX19gYGBhYWFiYmJjY2NkZ [...]
+<p><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAqAAAAKgCAMAAABz4j/3AAADAFBMVEUAAAABAQECAgIDAwMEBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUWFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJycoKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6Ojo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tMTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1eXl5fX19gYGBhYWFiYmJjY2NkZ [...]
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">edge.col <-<span class="st"> </span><span class="kw">createLabel</span>(Nnet, cnet, <span class="dt">label=</span><span class="st">"black"</span>, <span class="st">"edge"</span>, <span class="dt">nomatch=</span><span class="st">"grey"</span>)
z <-<span class="st"> </span><span class="kw">plot.networx</span>(Nnet, <span class="st">"2D"</span>, <span class="dt">show.edge.label =</span> <span class="ot">TRUE</span>, <span class="dt">edge.color=</span>edge.col)</code></pre></div>
<p><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAqAAAAKgCAMAAABz4j/3AAADAFBMVEUAAAABAQECAgIDAwMEBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUWFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJycoKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6Ojo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tMTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1eXl5fX19gYGBhYWFiYmJjY2NkZ [...]
diff --git a/inst/doc/Networx.Rmd b/inst/doc/Networx.Rmd
index 129ba78..05c69ab 100644
--- a/inst/doc/Networx.Rmd
+++ b/inst/doc/Networx.Rmd
@@ -11,7 +11,7 @@ vignette: >
---
-This tutorial gives a basic introduction for constructing phylogenetic networks and adding parameters to trees or networx objects using [phangorn](http://cran.r-project.org/package=phangorn) [@Schliep2011] in R.
+This tutorial gives a basic introduction for constructing phylogenetic networks and adding parameters to trees or networx objects using [phangorn](https://cran.r-project.org/package=phangorn) [@Schliep2011] in R.
Splits graphs or phylogenetic networks are a useful way to display conflicting data or to summarize different trees. Here, we present two popular networks, consensus networks [@Holland2004]
and Neighbor-Net [@Bryant2004].
Trees or networks are often missing either edge weights or edge support values. We show how to improve a tree/networx object by adding support values or estimating edge weights using non-negative Least-Squares (nnls).
diff --git a/inst/doc/Networx.html b/inst/doc/Networx.html
index 21e1b4a..c020946 100644
--- a/inst/doc/Networx.html
+++ b/inst/doc/Networx.html
@@ -12,7 +12,7 @@
<meta name="author" content="Klaus Schliep" />
-<meta name="date" content="2016-06-20" />
+<meta name="date" content="2016-12-03" />
<title>Splits and Networx</title>
@@ -70,11 +70,11 @@ code > span.in { color: #60a0b0; font-weight: bold; font-style: italic; } /* Inf
<h1 class="title toc-ignore">Splits and Networx</h1>
<h4 class="author"><em>Klaus Schliep</em></h4>
-<h4 class="date"><em>June 20, 2016</em></h4>
+<h4 class="date"><em>December 03, 2016</em></h4>
-<p>This tutorial gives a basic introduction for constructing phylogenetic networks and adding parameters to trees or networx objects using <a href="http://cran.r-project.org/package=phangorn">phangorn</a> <span class="citation">(Schliep 2011)</span> in R. Splits graphs or phylogenetic networks are a useful way to display conflicting data or to summarize different trees. Here, we present two popular networks, consensus networks <span class="citation">(Holland et al. 2004)</span> and Neigh [...]
+<p>This tutorial gives a basic introduction for constructing phylogenetic networks and adding parameters to trees or networx objects using <a href="https://cran.r-project.org/package=phangorn">phangorn</a> <span class="citation">(Schliep 2011)</span> in R. Splits graphs or phylogenetic networks are a useful way to display conflicting data or to summarize different trees. Here, we present two popular networks, consensus networks <span class="citation">(Holland et al. 2004)</span> and Neig [...]
Trees or networks are often missing either edge weights or edge support values. We show how to improve a tree/networx object by adding support values or estimating edge weights using non-negative Least-Squares (nnls).</p>
<p>We first load the phangorn package and a few data sets we use in this vignette.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">library</span>(phangorn)
@@ -93,7 +93,7 @@ tree <-<span class="st"> </span><span class="kw">plotBS</span>(tree, bs, <spa
<p><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAASAAAAEgCAMAAAAjXV6yAAAC91BMVEUAAAABAQECAgIDAwMEBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQWFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJycoKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6Ojo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tMTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1eXl5fX19gYGBhYWFiYmJjY2NkZGRlZ [...]
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">cnet <-<span class="st"> </span><span class="kw">consensusNet</span>(bs, .<span class="dv">3</span>)
<span class="kw">plot</span>(cnet, <span class="st">"2D"</span>, <span class="dt">show.edge.label=</span><span class="ot">TRUE</span>)</code></pre></div>
-<p><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAASAAAAEgCAMAAAAjXV6yAAADAFBMVEUAAAABAQECAgIDAwMEBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUWFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJycoKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6Ojo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tMTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1eXl5fX19gYGBhYWFiYmJjY2NkZ [...]
+<p><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAASAAAAEgCAMAAAAjXV6yAAAC/VBMVEUAAAABAQECAgIDAwMEBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUWFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJycoKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6Ojo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tMTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1eXl5fX19gYGBhYWFiYmJjY2NkZ [...]
<p>In many cases, <code>consensusNet</code> will return more than two incompatible (competing) splits. This cannot be plotted as a planar (2-dimensional) graph. Such as situation requires a n-dimensional graph, where the maximum number of dimensions equals the maximum number of incompatible splits. For example, if we have three alternative incompatible splits: (A,B)|(C,D) vs. (A,C)|(B,D) vs. (A,D)|(B,C), we need a 3-dimensional graph to show all three alternatives. A nice way to get stil [...]
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">plot</span>(cnet)
<span class="co"># rotate 3d plot</span>
@@ -113,7 +113,7 @@ tree <-<span class="st"> </span><span class="kw">plotBS</span>(tree, bs, <spa
nnet <-<span class="st"> </span><span class="kw">neighborNet</span>(dm)
<span class="kw">par</span>(<span class="st">"mar"</span> =<span class="st"> </span><span class="kw">rep</span>(<span class="dv">1</span>, <span class="dv">4</span>))
<span class="kw">plot</span>(nnet, <span class="st">"2D"</span>)</code></pre></div>
-<p><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAASAAAAEgCAMAAAAjXV6yAAAC+lBMVEUAAAABAQECAgIDAwMEBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUWFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJycoKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6Ojo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tMTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1eXl5fX19gYGBhYWFiYmJjY2NkZ [...]
+<p><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAASAAAAEgCAMAAAAjXV6yAAAC/VBMVEUAAAABAQECAgIDAwMEBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUWFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJycoKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6Ojo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tMTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1eXl5fX19gYGBhYWFiYmJjY2NkZ [...]
<p>The advantage of Neighbor-Net is that it returns a circular split system which can be always displayed in a planar (2D) graph. The plots displayed in <code>phangorn</code> may (however) not be planar, but re-plotting can give you a planar graph. This unwanted behavior will be improved in future version. The rendering of the <code>networx</code> is done using the the fantastic igraph package <span class="citation">(Csardi and Nepusz 2006)</span>.</p>
</div>
<div id="adding-support-values" class="section level2">
@@ -122,14 +122,14 @@ nnet <-<span class="st"> </span><span class="kw">neighborNet</span>(dm)
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">nnet <-<span class="st"> </span><span class="kw">addConfidences</span>(nnet, tree)
<span class="kw">par</span>(<span class="st">"mar"</span> =<span class="st"> </span><span class="kw">rep</span>(<span class="dv">1</span>, <span class="dv">4</span>))
<span class="kw">plot</span>(nnet, <span class="st">"2D"</span>, <span class="dt">show.edge.label=</span><span class="ot">TRUE</span>)</code></pre></div>
-<p><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAASAAAAEgCAMAAAAjXV6yAAADAFBMVEUAAAABAQECAgIDAwMEBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUWFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJycoKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6Ojo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tMTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1eXl5fX19gYGBhYWFiYmJjY2NkZ [...]
+<p><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAASAAAAEgCAMAAAAjXV6yAAAC/VBMVEUAAAABAQECAgIDAwMEBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUWFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJycoKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6Ojo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tMTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1eXl5fX19gYGBhYWFiYmJjY2NkZ [...]
<p>Analogously, we can also add support values to a tree:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">tree2 <-<span class="st"> </span><span class="kw">rNNI</span>(tree, <span class="dv">2</span>)
tree2 <-<span class="st"> </span><span class="kw">addConfidences</span>(tree2, tree)
<span class="co"># several support values are missing</span>
<span class="kw">par</span>(<span class="st">"mar"</span> =<span class="st"> </span><span class="kw">rep</span>(<span class="dv">1</span>, <span class="dv">4</span>))
<span class="kw">plot</span>(tree2, <span class="dt">show.node.label=</span><span class="ot">TRUE</span>)</code></pre></div>
-<p><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAASAAAAEgCAMAAAAjXV6yAAAC61BMVEUAAAABAQECAgIDAwMEBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUWFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJycoKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6Ojo7Ozs8PDw+Pj4/Pz9AQEBBQUFCQkJERERFRUVGRkZHR0dISEhJSUlKSkpLS0tMTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1eXl5fX19gYGBhYWFiYmJjY2NkZGRlZWVnZ [...]
+<p><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAASAAAAEgCAMAAAAjXV6yAAAC7lBMVEUAAAABAQECAgIDAwMEBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUWFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJycoKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6Ojo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJERERFRUVGRkZISEhJSUlKSkpLS0tMTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1eXl5fX19gYGBhYWFiYmJjY2NkZGRlZWVnZ [...]
</div>
<div id="estimating-edge-weights-nnls" class="section level2">
<h2>Estimating edge weights (nnls)</h2>
@@ -137,7 +137,7 @@ tree2 <-<span class="st"> </span><span class="kw">addConfidences</span>(tree2
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">cnet <-<span class="st"> </span><span class="kw">nnls.networx</span>(cnet, dm)
<span class="kw">par</span>(<span class="st">"mar"</span> =<span class="st"> </span><span class="kw">rep</span>(<span class="dv">1</span>, <span class="dv">4</span>))
<span class="kw">plot</span>(cnet, <span class="st">"2D"</span>, <span class="dt">show.edge.label=</span><span class="ot">TRUE</span>)</code></pre></div>
-<p><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAASAAAAEgCAMAAAAjXV6yAAADAFBMVEUAAAABAQECAgIDAwMEBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUWFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJycoKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6Ojo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tMTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1eXl5fX19gYGBhYWFiYmJjY2NkZ [...]
+<p><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAASAAAAEgCAMAAAAjXV6yAAADAFBMVEUAAAABAQECAgIDAwMEBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUWFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJycoKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6Ojo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tMTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1eXl5fX19gYGBhYWFiYmJjY2NkZ [...]
<div id="import-and-export-networks-advanced-functions-for-networx-objects" class="section level3">
<h3>Import and export networks, advanced functions for networx objects</h3>
<p>The functions <code>read.nexus.networx</code> and <code>write.nexus.networx</code> can read and write nexus files for or from SplitsTree <span class="citation">(Huson and Bryant 2006)</span>. Check-out the new vignette IntertwiningTreesAndNetworks for additional functions, examples, and advanced applications.</p>
diff --git a/inst/doc/Trees.pdf b/inst/doc/Trees.pdf
index dc9a6bf..40c819c 100644
Binary files a/inst/doc/Trees.pdf and b/inst/doc/Trees.pdf differ
diff --git a/inst/doc/phangorn-specials.pdf b/inst/doc/phangorn-specials.pdf
index d5653c8..7ec2b88 100644
Binary files a/inst/doc/phangorn-specials.pdf and b/inst/doc/phangorn-specials.pdf differ
diff --git a/man/SOWH.test.Rd b/man/SOWH.test.Rd
index 19076c2..e061d3e 100644
--- a/man/SOWH.test.Rd
+++ b/man/SOWH.test.Rd
@@ -21,7 +21,7 @@ It makes extensive use \code{simSeq} and \code{optim.pml} and can take quite lon
}
\value{
an object of class SOWH. That is a list with three elements, one is a matrix
- containing for each bootstrap replicate the (log-) likelihood of the restricted and unrestricted estimate and two pml objetcs of the restricted and unrestricted model.
+ containing for each bootstrap replicate the (log-) likelihood of the restricted and unrestricted estimate and two pml objects of the restricted and unrestricted model.
}
\references{
Goldman, N., Anderson, J. P., and Rodrigo, A. G. (2000) Likelihood
diff --git a/man/ancestral.pml.Rd b/man/ancestral.pml.Rd
index cf475c1..4ba4a2a 100644
--- a/man/ancestral.pml.Rd
+++ b/man/ancestral.pml.Rd
@@ -11,7 +11,7 @@ Ancestral character reconstruction.
Marginal reconstruction of the ancestral character states.
}
\usage{
-ancestral.pml(object, type = c("ml", "bayes"))
+ancestral.pml(object, type = c("marginal", "ml", "bayes"))
ancestral.pars(tree, data, type = c("MPR", "ACCTRAN"), cost = NULL)
pace(tree, data, type = c("MPR", "ACCTRAN"), cost = NULL)
plotAnc(tree, data, i, col=NULL, cex.pie=par("cex"), pos="bottomright", ...)
diff --git a/man/as.splits.Rd b/man/as.splits.Rd
index 60e75f5..93041a0 100644
--- a/man/as.splits.Rd
+++ b/man/as.splits.Rd
@@ -7,9 +7,12 @@
\alias{as.matrix.splits}
\alias{as.Matrix}
\alias{as.Matrix.splits}
+\alias{c.splits}
+\alias{distinct.splits}
\alias{print.splits}
\alias{write.splits}
\alias{allSplits}
+\alias{allCircularSplits}
\alias{compatible}
\alias{write.nexus.splits}
\alias{read.nexus.splits}
@@ -33,15 +36,17 @@ as.splits(x, ...)
\method{as.prop.part}{splits}(x, ...)
compatible(obj)
allSplits(k, labels = NULL)
+allCircularSplits(k, labels = NULL)
write.nexus.splits(obj, file = "", weights=NULL, taxa=TRUE, append=FALSE)
read.nexus.splits(file)
addTrivialSplits(obj)
+distinct.splits(...)
}
%- maybe also 'usage' for other objects documented here.
\arguments{
\item{x}{An object of class phylo or multiPhylo.}
\item{maxp}{integer, default from \code{options(max.print)}, influences how many entries of large matrices are printed at all.}
- \item{zero.print}{character which should be printed for zeroes.}
+ \item{zero.print}{character which should be printed for zeros.}
\item{one.print}{character which should be printed for ones.}
\item{\dots}{Further arguments passed to or from other methods.}
\item{obj}{an object of class splits.}
@@ -75,6 +80,8 @@ The internal representation is likely to change.
\examples{
(sp <- as.splits(rtree(5)))
write.nexus.splits(sp)
+spl <- allCircularSplits(5)
+plot(as.networx(spl), "2D")
}
% Add one or more standard keywords, see file 'KEYWORDS' in the
% R documentation directory.
diff --git a/man/densiTree.Rd b/man/densiTree.Rd
index c094ed5..d9023d9 100644
--- a/man/densiTree.Rd
+++ b/man/densiTree.Rd
@@ -1,6 +1,5 @@
\name{densiTree}
\alias{densiTree}
-%- Also NEED an '\alias' for EACH other topic documented here.
\title{
Plots a densiTree.
}
@@ -11,7 +10,6 @@ An R function to plot trees similar to those produced by DensiTree.
densiTree(x, type = "cladogram", alpha = 1/length(x), consensus = NULL, optim = FALSE,
scaleX = FALSE, col = 1, width = 1, cex = 0.8, ...)
}
-%- maybe also 'usage' for other objects documented here.
\arguments{
\item{x}{
an object of class \code{multiPhylo}.
@@ -58,10 +56,9 @@ Remco R. Bouckaert (2010) DensiTree: making sense of sets of phylogenetic trees
Klaus Schliep \email{klaus.schliep at gmail.com}
}
-%% ~Make other sections like Warning with \section{Warning }{....} ~
\seealso{
-\code{\link{plot.phylo}}, \code{\link{plot.networx}}
+\code{\link{plot.phylo}}, \code{\link{plot.networx}}
}
\examples{
data(Laurasiatherian)
@@ -72,7 +69,7 @@ bs <- bootstrap.phyDat(Laurasiatherian, FUN =
densiTree(bs, optim=TRUE, type="cladogram", col="blue")
densiTree(bs, optim=TRUE, type="phylogram", col="green")
\dontrun{
-# phylogram are nice to show different age estimates
+# phylograms are nice to show different age estimates
require(PhyloOrchard)
data(BinindaEmondsEtAl2007)
BinindaEmondsEtAl2007 <- .compressTipLabel(BinindaEmondsEtAl2007)
diff --git a/man/dist.p.Rd b/man/dist.p.Rd
index a84748e..2b2938c 100644
--- a/man/dist.p.Rd
+++ b/man/dist.p.Rd
@@ -14,7 +14,7 @@ dist.p(x, cost="polymorphism", ignore.indels=TRUE)
\arguments{
\item{x}{a matrix containing DNA sequences; this must be of class
"phyDat" (use as.phyDat to convert from DNAbin objects).}
- \item{cost}{A cost matrix or "polymorphism" for a pre defined one.}
+ \item{cost}{A cost matrix or "polymorphism" for a predefined one.}
\item{ignore.indels}{a logical indicating whether gaps are treated
as fifth state or not. Warning, each gap site is treated as a
characters, so an an indel that spans a number of base positions
diff --git a/man/hadamard.Rd b/man/hadamard.Rd
index 298a861..2805370 100644
--- a/man/hadamard.Rd
+++ b/man/hadamard.Rd
@@ -50,24 +50,22 @@ LogDet transforms, and maximum likelihood. \emph{PhD thesis}.
\seealso{\code{\link{distanceHadamard}}, \code{\link{lento}}, \code{\link{plot.networx}}}
\examples{
-H = hadamard(3)
-v = 1:8
+H <- hadamard(3)
+v <- 1:8
H%*%v
fhm(v)
data(yeast)
-dat = as.character(yeast)
+
# RY-coding
-dat2 = dat
-dat2[dat=="a" | dat=="g"] = "r"
-dat2[dat=="c" | dat=="t"] = "y"
-dat2 = phyDat(dat2, type="USER", levels=c("r","y"), ambiguity=NULL)
-fit2 = h2st(dat2)
+dat_ry <- acgt2ry(yeast)
+fit2 <- h2st(dat_ry)
lento(fit2)
# write.nexus.splits(fit2, file = "test.nxs")
# read this file into Spectronet or Splitstree to show the network
\dontrun{
+dat = as.character(yeast)
dat4 = phyDat(dat, type="USER", levels=c("a","c", "g", "t"), ambiguity=NULL)
fit4 = h4st(dat4)
diff --git a/man/maxCladeCred.Rd b/man/maxCladeCred.Rd
index 5bfce57..4aa8ddc 100644
--- a/man/maxCladeCred.Rd
+++ b/man/maxCladeCred.Rd
@@ -27,7 +27,7 @@ logical, if FALSE the tree with highest maximum bipartition credibility is retur
\details{
So far just the best tree is returned. No annotations or transformations of edge length are performed.
-If a list of partition is provided then the clade credibity is computed for the trees in x.
+If a list of partition is provided then the clade credibility is computed for the trees in x.
}
\value{
a tree (an object of class \code{phylo}) with the highest clade credibility or a numeric vector of clade credibilities for each tree.
diff --git a/man/midpoint.Rd b/man/midpoint.Rd
index 9da9268..059e8b9 100644
--- a/man/midpoint.Rd
+++ b/man/midpoint.Rd
@@ -22,7 +22,7 @@ getRoot(tree)
}
\details{
\code{pruneTree} prunes back a tree and produces a consensus tree, for trees already containing nodelabels.
-It assumes that nodelabels are numerical or character genereated from numericals, it uses as.numeric(as.character(tree$node.labels))
+It assumes that nodelabels are numerical or character that allows conversion to numerical, it uses as.numeric(as.character(tree$node.labels))
to convert them.
\code{midpoint} so far does not transform node.labels properly.
}
diff --git a/man/modelTest.Rd b/man/modelTest.Rd
index ddb924f..62abb2f 100644
--- a/man/modelTest.Rd
+++ b/man/modelTest.Rd
@@ -5,7 +5,7 @@
ModelTest
}
\description{
-Comparison of different nucleotide or amino acid substition models
+Comparison of different nucleotide or amino acid substitution models
}
\usage{
modelTest(object, tree=NULL, model = c("JC", "F81", "K80", "HKY", "SYM", "GTR"),
@@ -18,7 +18,7 @@ modelTest(object, tree=NULL, model = c("JC", "F81", "K80", "HKY", "SYM", "GTR"),
\item{tree}{a phylogenetic tree.}
\item{model}{a vector containing the substitution models to compare with each other
or "all" to test all available models}
- \item{G}{logical, TRUE (default) if (discrete) Gamma modelshould be tested}
+ \item{G}{logical, TRUE (default) if (discrete) Gamma model should be tested}
\item{I}{logical, TRUE (default) if invariant sites should be tested}
\item{FREQ}{logical, FALSE (default) if TRUE amino acid frequencies will be estimated.}
\item{k}{number of rate classes}
diff --git a/man/nni.Rd b/man/nni.Rd
index a2ceb48..e1ae201 100644
--- a/man/nni.Rd
+++ b/man/nni.Rd
@@ -25,7 +25,7 @@ rNNI(tree, moves=1, n=length(moves))
}
\author{ Klaus Schliep \email{klaus.schliep at gmail.com} }
-\seealso{\code{\link{allTrees}}}
+\seealso{\code{\link{allTrees}}, \code{\link{SPR.dist}}}
\examples{
tree = unroot(rtree(20))
trees1 <- nni(tree)
diff --git a/man/parsimony.Rd b/man/parsimony.Rd
index f1b7b9c..5022126 100644
--- a/man/parsimony.Rd
+++ b/man/parsimony.Rd
@@ -14,7 +14,7 @@
\code{parsimony} returns the parsimony score of a tree using either the sankoff or the fitch algorithm.
\code{optim.parsimony} tries to find the maximum parsimony tree using either Nearest Neighbor Interchange (NNI)
-rearrangements or sub tree pruning and regrafting (SPR). \code{pratchet} implements the parsimony ratchet (Nixon, 1999) and is the prefered way to search for the best tree.
+rearrangements or sub tree pruning and regrafting (SPR). \code{pratchet} implements the parsimony ratchet (Nixon, 1999) and is the preferred way to search for the best tree.
\code{random.addition} can be used to produce starting trees.
\code{CI} and \code{RI} computes the consistency and retention index.
}
@@ -22,7 +22,7 @@ rearrangements or sub tree pruning and regrafting (SPR). \code{pratchet} impleme
parsimony(tree, data, method="fitch", ...)
optim.parsimony(tree, data, method="fitch", cost=NULL, trace=1, rearrangements="SPR", ...)
pratchet(data, start=NULL, method="fitch", maxit=1000, k=10, trace=1, all=FALSE,
- rearrangements="SPR", ...)
+ rearrangements="SPR", perturbation="ratchet", ...)
fitch(tree, data, site = "pscore")
sankoff(tree, data, cost = NULL, site = "pscore")
random.addition(data, method="fitch")
@@ -43,6 +43,7 @@ acctran(tree, data)
\item{maxit}{maximum number of iterations in the ratchet.}
\item{k}{number of rounds ratchet is stopped, when there is no improvement.}
\item{all}{return all equally good trees or just one of them.}
+ \item{perturbation}{whether using a ratchet or stochastic (nni) for shuffling the tree.}
\item{...}{Further arguments passed to or from other methods (e.g. model="sankoff" and cost matrix).}
\item{sitewise}{return CI/RI for alignment or sitewise}
}
diff --git a/man/phyDat.Rd b/man/phyDat.Rd
index 75aeb0a..e75a439 100644
--- a/man/phyDat.Rd
+++ b/man/phyDat.Rd
@@ -3,7 +3,12 @@
\alias{as.AAbin.phyDat}
\alias{as.phyDat.DNAbin}
\alias{as.phyDat.alignment}
+\alias{as.phyDat.character}
+\alias{as.phyDat.data.frame}
+\alias{as.phyDat.factor}
+\alias{as.phyDat.matrix}
\alias{as.phyDat.MultipleAlignment}
+\alias{as.MultipleAlignment}
\alias{as.MultipleAlignment.phyDat}
\alias{as.data.frame.phyDat}
\alias{as.character.phyDat}
@@ -22,6 +27,8 @@
\alias{phyDat2alignment}
\alias{phyDat2MultipleAlignment}
\alias{image.phyDat}
+\alias{dna2codon}
+\alias{codon2dna}
\title{Conversion among Sequence Formats}
\description{
These functions transform several DNA formats into the \code{phyDat} format.
@@ -30,9 +37,13 @@ These functions transform several DNA formats into the \code{phyDat} format.
\usage{
phyDat(data, type = "DNA", levels = NULL, return.index=TRUE, ...)
read.phyDat(file, format="phylip", type="DNA", ...)
-write.phyDat(x, file, format="phylip",...)
+write.phyDat(x, file, format="phylip", colsep = "", nbcol=-1, ...)
\method{as.phyDat}{DNAbin}(x, ...)
\method{as.phyDat}{alignment}(x, type="DNA", ...)
+\method{as.phyDat}{character}(x, ...)
+\method{as.phyDat}{data.frame}(x, ...)
+\method{as.phyDat}{factor}(x, ...)
+\method{as.phyDat}{matrix}(x, ...)
\method{as.phyDat}{MultipleAlignment}(x, ...)
%\method{as.AAbin}{phyDat}(x, ...)
\method{as.character}{phyDat}(x, allLevels = TRUE, ...)
@@ -44,6 +55,8 @@ phyDat2alignment(x)
allSitePattern(n, levels=c("a","c","g","t"), names=NULL)
acgt2ry(obj)
baseFreq(obj, freq=FALSE, all=FALSE, drop.unused.levels=FALSE)
+dna2codon(x)
+codon2dna(x)
}
%- maybe also 'usage' for other objects documented here.
\arguments{
@@ -56,6 +69,8 @@ baseFreq(obj, freq=FALSE, all=FALSE, drop.unused.levels=FALSE)
\item{format}{File format of the sequence alignment (see details).
Several popular formats are supported: "phylip", "interleaved", "sequential", "clustal", "fasta" or "nexus", or any unambiguous abbreviation of these.
}
+ \item{colsep}{a character used to separate the columns (a single space by default).}
+ \item{nbcol}{a numeric specifying the number of columns per row (-1 by default); may be negative implying that the nucleotides are printed on a single line.}
\item{n}{Number of sequences.}
\item{names}{Names of sequences.}
\item{subset}{a subset of taxa.}
@@ -66,7 +81,7 @@ baseFreq(obj, freq=FALSE, all=FALSE, drop.unused.levels=FALSE)
\item{freq}{logical, if 'TRUE', frequencies or counts are returned otherwise proportions}
\item{all}{all a logical; if all = TRUE, all counts of bases, ambiguous codes, missing data, and alignment gaps are returned as defined in the contrast.}
\item{drop.unused.levels}{logical, drop unused levels}
- \item{incomparables}{for compatability with unique.}
+ \item{incomparables}{for compatibility with unique.}
\item{identical}{if TRUE (default) sequences have to be identical, if FALSE sequences are considered
duplicates if distance between sequences is zero (happens frequently with ambiguous sites).}
\item{...}{further arguments passed to or from other methods.}
diff --git a/man/plot.networx.Rd b/man/plot.networx.Rd
index 6ec1337..beb9c4d 100644
--- a/man/plot.networx.Rd
+++ b/man/plot.networx.Rd
@@ -120,7 +120,7 @@ cyclic ordering. These objects can be nicely plotted in "2D".
So far not all parameters behave the same on the the rgl "3D"
and basic graphic "2D" device.
-Often it is easier (and safer) to supply vectors of graphial parameters for splits (e.g. splits.color).
+Often it is easier (and safer) to supply vectors of graphical parameters for splits (e.g. splits.color).
These overwrite values edge.color.
}
\note{
diff --git a/man/pml.Rd b/man/pml.Rd
index dbf5196..e98fb32 100644
--- a/man/pml.Rd
+++ b/man/pml.Rd
@@ -131,10 +131,12 @@ L.-T. Nguyen, H.A. Schmidt, A. von Haeseler, and B.Q. Minh (2015) IQ-TREE: A fas
Vos, R. A. (2003) Accelerated Likelihood Surface Exploration: The Likelihood Ratchet. \emph{Systematic Biology}, \bold{52(3)}, 368--373
Yang, Z., and R. Nielsen (1998) Synonymous and nonsynonymous rate variation in nuclear genes of mammals. \emph{Journal of Molecular Evolution}, \bold{46}, 409-418.
+
+Lewis, P.O. (2001) A likelihood approach to estimating phylogeny from discrete morphological character data. \emph{Systematic Biology} \bold{50}, 913--925.
}
\author{Klaus Schliep \email{klaus.schliep at gmail.com}}
\seealso{
-\code{\link{bootstrap.pml}}, \code{\link{modelTest}}, \code{\link{pmlPart}}, \code{\link{pmlMix}}, \code{\link{plot.phylo}}, \code{\link{SH.test}}
+\code{\link{bootstrap.pml}}, \code{\link{modelTest}}, \code{\link{pmlPart}}, \code{\link{pmlMix}}, \code{\link{plot.phylo}}, \code{\link{SH.test}}, \code{\link{ancestral.pml}}
}
% \note{For small trees the likelihood seems to be very similar to Paup* or PhyML.}
\examples{
diff --git a/man/pml.fit.Rd b/man/pml.fit.Rd
index e19952f..ddcc3f4 100644
--- a/man/pml.fit.Rd
+++ b/man/pml.fit.Rd
@@ -32,34 +32,19 @@ discrete.gamma(alpha, k)
\item{alpha}{Shape parameter of the gamma distribution.}
\item{k}{Number of intervals of the discrete gamma distribution.}
\item{Q}{A vector containing the lower triangular part of the rate matrix.}
- \item{levels}{
-%% ~~Describe \code{levels} here~~
-}
+ \item{levels}{The alphabet used e.g. c("a", "c", "g", "t") for DNA}
\item{inv}{Proportion of invariable sites.}
\item{rate}{Rate.}
- \item{g}{
-%% ~~Describe \code{g} here~~
-}
- \item{w}{
-%% ~~Describe \code{w} here~~
-}
+ \item{g}{vector of quantiles (default is NULL)}
+ \item{w}{vector of probabilities (default is NULL)}
\item{eig}{Eigenvalue decomposition of Q}
\item{INV}{Sparse representation of invariant sites}
- \item{ll.0}{
-%% ~~Describe \code{ll.0} here~~
- }
- \item{llMix}{
-%% ~~Describe \code{llMix} here~~
- }
- \item{wMix}{
-%% ~~Describe \code{wMix} here~~
-
- }
+ \item{ll.0}{default is NULL}
+ \item{llMix}{default is NULL}
+ \item{wMix}{default is NULL}
\item{\dots}{Further arguments passed to or from other methods.}
- \item{site}{
-%% ~~Describe \code{site} here~~
-}
-}
+ \item{site}{return the log-likelihood or vector of sitewise likelihood values}
+}
\details{
These functions are exported to be used in different packages so far only in the package coalescentMCMC, but are not intended for end user. Most of the functions call C code and are far less forgiving if the import is not what they expect than \code{pml}.
}
diff --git a/man/simSeq.Rd b/man/simSeq.Rd
index 25e990b..cf8de37 100644
--- a/man/simSeq.Rd
+++ b/man/simSeq.Rd
@@ -24,7 +24,8 @@ other root sequence is randomly generated.}
\item{model}{Amino acid models: one of "WAG", "JTT", "Dayhoff" or "LG"}
\item{levels}{ \code{levels} takes a character vector of the different bases,
default is for nucleotide sequences, only used when type = "USER".}
-\item{rate}{rate. }
+\item{rate}{mutation rate or scaler for the edge length,
+a numerical value greater than zero. }
\item{ancestral}{Return ancestral sequences?}
\item{\dots}{Further arguments passed to or from other methods.}
}
@@ -33,7 +34,7 @@ default is for nucleotide sequences, only used when type = "USER".}
It is quite flexible and allows to generate DNA, RNA, amino acids or binary sequences.
It is possible to give a \code{pml} object as input simSeq return a \code{phyDat}
from these model.
-There is also a more low level version, which lacks rate variation, but one can combine different alignments having their own rate (see example).
+There is also a more low level version, which lacks rate variation, but one can combine different alignments having their own rate (see example). The rate parameter acts like a scaler for the edge lengths.
}
\value{
\code{simSeq} returns an object of class phyDat.
diff --git a/man/superTree.Rd b/man/superTree.Rd
index f920dab..293b279 100644
--- a/man/superTree.Rd
+++ b/man/superTree.Rd
@@ -6,9 +6,9 @@
Super Tree and Species Tree methods
}
\description{
-These function \code{superTree} allows the estimation of a rooted supertree from a set of rooted trees using Matrix representation parsimony. \code{coalSpeciesTree} estimates species trees and can multiple individuals per species.}
+These function \code{superTree} allows the estimation of a supertree from a set of trees using either Matrix representation parsimony, Robinson-Foulds or SPR as criterion. \code{coalSpeciesTree} estimates species trees and can multiple individuals per species.}
\usage{
-superTree(tree, method = "pratchet", rooted=TRUE, ...)
+superTree(tree, method = "MRP", rooted=FALSE, trace=0, ...)
coalSpeciesTree(tree, X, sTree = NULL)
}
%- maybe also 'usage' for other objects documented here.
@@ -17,7 +17,8 @@ coalSpeciesTree(tree, X, sTree = NULL)
an object of class \code{multiPhylo}
}
\item{method}{
-An argument defining which algorithm is used to optimize the tree.
+An argument defining which algorithm is used to optimize the tree.
+Possible are "MRP", "NNI", and "SPR".
}
\item{rooted}{
should the resulting supertrees be rooted.
@@ -25,6 +26,9 @@ should the resulting supertrees be rooted.
\item{X}{
A \code{phyDat} object to define which individual belongs to which species.
}
+\item{trace}{
+defines how much information is printed during optimization.
+}
\item{sTree}{
A species tree which forces the topology.
}
@@ -35,7 +39,7 @@ further arguments passed to or from other methods.
\details{
The function \code{superTree} extends the function mrp.supertree from Liam Revells,
with artificial adding an outgroup on the root of the trees.
-This allows to root the supertree afterwards. The functions is internally used in DensiTree.
+This allows to root the supertree afterwards. The functions is internally used in DensiTree. The implementation for the RF- and SPR-supertree are very basic so far and assume that all trees share the same set of taxa.
\code{coalSpeciesTree} estimates a single linkage tree as suggested by Liu et al. (2010) from the element wise minima of the cophenetic matrices of the gene trees. It extends \code{speciesTree} in ape as it allows that have several individuals per gene tree.
}
@@ -55,14 +59,20 @@ Emmanuel Paradies
%% ~Make other sections like Warning with \section{Warning }{....} ~
\seealso{
-\code{mrp.supertree}, \code{\link{speciesTree}}, \code{\link{densiTree}}, \code{\link{hclust}}
+\code{mrp.supertree}, \code{\link{speciesTree}}, \code{\link{densiTree}}, \code{\link{RF.dist}}, \code{\link{SPR.dist}}
}
\examples{
data(Laurasiatherian)
set.seed(1)
bs <- bootstrap.phyDat(Laurasiatherian, FUN = function(x)upgma(dist.hamming(x)), bs=50)
class(bs) <- 'multiPhylo'
-plot(superTree(bs))
+
+mrp_st <- superTree(bs, rooted=TRUE)
+plot(superTree(mrp_st))
+\dontrun{
+rf_st <- superTree(bs, method = "RF")
+spr_st <- superTree(bs, method = "SPR")
+}
}
% Add one or more standard keywords, see file 'KEYWORDS' in the
% R documentation directory.
diff --git a/man/treedist.Rd b/man/treedist.Rd
index 109f462..b36fe09 100644
--- a/man/treedist.Rd
+++ b/man/treedist.Rd
@@ -4,6 +4,8 @@
\alias{wRF.dist}
\alias{KF.dist}
\alias{path.dist}
+\alias{sprdist}
+\alias{SPR.dist}
%\alias{print.treedist}
\title{ Distances between trees }
\description{
@@ -11,9 +13,11 @@
}
\usage{
treedist(tree1, tree2, check.labels = TRUE)
+sprdist(tree1, tree2)
RF.dist(tree1, tree2=NULL, normalize=FALSE, check.labels=TRUE, rooted=FALSE)
wRF.dist(tree1, tree2=NULL, normalize=FALSE, check.labels=TRUE, rooted=FALSE)
KF.dist(tree1, tree2=NULL, check.labels=TRUE, rooted=FALSE)
+SPR.dist(tree1, tree2=NULL)
path.dist(tree1, tree2=NULL, check.labels=TRUE, use.weight=FALSE)
}
\arguments{
@@ -52,6 +56,8 @@ where \eqn{i(T_1)} denotes the number of internal edges and \eqn{v_s(T_1, T_2)}
If edge weights should be considered \code{wRF.dist} calculates the weighted RF distance (Robinson & Foulds 1981). and \code{KF.dist} calculates the branch score distance (Kuhner & Felsenstein 1994).
\code{path.dist} computes the path difference metric as described in Steel and Penny 1993).
+\code{sprdist} computes the approximate SPR distance (Oliveira Martins et al. 2008, de Oliveira Martins 2016).
+
For large number of trees the distance functions can use a lot of memory!
% The function used internally is 2 * (nt - m) where nt is the number of tips and
@@ -60,6 +66,10 @@ For large number of trees the distance functions can use a lot of memory!
}
\references{
+de Oliveira Martins L., Leal E., Kishino H. (2008) \emph{Phylogenetic Detection of Recombination with a Bayesian Prior on the Distance between Trees}. PLoS ONE \bold{3(7)}. e2651. doi: 10.1371/journal.pone.0002651
+
+de Oliveira Martins L., Mallo D., Posada D. (2016) \emph{A Bayesian Supertree Model for Genome-Wide Species Tree Reconstruction}. Syst. Biol. \bold{65(3)}: 397-416, doi:10.1093/sysbio/syu082
+
Steel M. A. and Penny P. (1993) \emph{Distributions of tree comparison metrics - some new results}, Syst. Biol., \bold{42(2)}, 126--141
Kuhner, M. K. and Felsenstein, J. (1994) \emph{A simulation comparison of phylogeny algorithms under equal and unequal evolutionary rates}, Molecular Biology and Evolution, \bold{11(3)}, 459--468
@@ -68,12 +78,18 @@ D.F. Robinson and L.R. Foulds (1981) \emph{Comparison of phylogenetic trees}, Ma
D.F. Robinson and L.R. Foulds (1979) Comparison of weighted labelled trees. In Horadam, A. F. and Wallis, W. D. (Eds.), \emph{Combinatorial Mathematics VI: Proceedings of the Sixth Australian Conference on Combinatorial Mathematics, Armidale, Australia}, 119--126
}
-\author{ Klaus P. Schliep \email{klaus.schliep at gmail.com}}
-\seealso{\code{\link[ape]{dist.topo}}}
+\author{ Klaus P. Schliep \email{klaus.schliep at gmail.com},
+
+Leonardo de Oliveira Martins
+}
+\seealso{\code{\link[ape]{dist.topo}}, \code{\link{nni}}, \code{\link{superTree}}}
\examples{
tree1 <- rtree(100, rooted=FALSE)
tree2 <- rSPR(tree1, 3)
RF.dist(tree1, tree2)
treedist(tree1, tree2)
+sprdist(tree1, tree2)
+trees <- rSPR(tree1, 1:5)
+SPR.dist(tree1, trees)
}
\keyword{ classif }% at least one, from doc/KEYWORDS
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
new file mode 100644
index 0000000..c4efb1d
--- /dev/null
+++ b/src/RcppExports.cpp
@@ -0,0 +1,30 @@
+// Generated by using Rcpp::compileAttributes() -> do not edit by hand
+// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+#include <Rcpp.h>
+
+using namespace Rcpp;
+
+// allDescCPP
+List allDescCPP(IntegerMatrix orig, int nTips);
+RcppExport SEXP phangorn_allDescCPP(SEXP origSEXP, SEXP nTipsSEXP) {
+BEGIN_RCPP
+ Rcpp::RObject rcpp_result_gen;
+ Rcpp::RNGScope rcpp_rngScope_gen;
+ Rcpp::traits::input_parameter< IntegerMatrix >::type orig(origSEXP);
+ Rcpp::traits::input_parameter< int >::type nTips(nTipsSEXP);
+ rcpp_result_gen = Rcpp::wrap(allDescCPP(orig, nTips));
+ return rcpp_result_gen;
+END_RCPP
+}
+// allChildrenCPP
+List allChildrenCPP(IntegerMatrix orig);
+RcppExport SEXP phangorn_allChildrenCPP(SEXP origSEXP) {
+BEGIN_RCPP
+ Rcpp::RObject rcpp_result_gen;
+ Rcpp::RNGScope rcpp_rngScope_gen;
+ Rcpp::traits::input_parameter< IntegerMatrix >::type orig(origSEXP);
+ rcpp_result_gen = Rcpp::wrap(allChildrenCPP(orig));
+ return rcpp_result_gen;
+END_RCPP
+}
diff --git a/src/dist.c b/src/dist.c
index 8577ec2..0217938 100644
--- a/src/dist.c
+++ b/src/dist.c
@@ -102,10 +102,11 @@ SEXP PWI(SEXP LEFT, SEXP RIGHT, SEXP L, SEXP N, SEXP W, SEXP LI){
void C_fhm(double *v, int *n){
unsigned int level, i, j;
unsigned int start, step, num_splits;
+ unsigned int max_n = (unsigned int)*n;
double vi, vj;
num_splits = (1 << (*n));
step = 1;
- for(level = 0; level < (*n); level++){
+ for(level = 0; level < max_n; level++){
start = 0;
while(start < (num_splits-1)){
for(i = start; i < (start + step); i++){
@@ -157,8 +158,8 @@ void distance_hadamard(double *d, int n) {
d[0] = 0.0;
}
-
-void pairwise_distances(double *dm, int n, int num_splits, double *d) {
+// int num_splits raus
+void pairwise_distances(double *dm, int n, double *d) {
int k=0;
unsigned int offset;
for (int i = 0; i < (n-1); ++i) {
@@ -182,7 +183,7 @@ SEXP dist2spectra(SEXP dm, SEXP nx, SEXP ns) {
SEXP result;
PROTECT(result = allocVector(REALSXP, nsp));
res = REAL(result);
- pairwise_distances(REAL(dm), n, nsp, res);
+ pairwise_distances(REAL(dm), n, res); //nsp,
distance_hadamard(res, n);
UNPROTECT(1);
return(result);
diff --git a/src/fitch.c b/src/fitch.c
index ca4e0d6..5dac621 100644
--- a/src/fitch.c
+++ b/src/fitch.c
@@ -256,7 +256,9 @@ SEXP FITCH(SEXP dat, SEXP nrx, SEXP node, SEXP edge, SEXP l, SEXP weight, SEXP m
/*
ACCTRAN
*/
-void fitchT(int *dat1, int *dat2, int *nr, double *pars, double *weight, double *w){
+
+// , double *pars, double *weight, double *w raus
+void fitchT(int *dat1, int *dat2, int *nr){
int k;
int tmp;
for(k = 0; k < (*nr); k++){
@@ -290,6 +292,7 @@ void fitchT3(int *dat1, int *dat2, int *nr, double *pars, double *weight, double
// return lower and upper bound for the number of changes
// upper bound very conservative
+/*
void countMPR(double *res, int *dat1, int *dat2, int *nr, double *weight, int *external){
int k;
int tmp;
@@ -310,18 +313,21 @@ void countMPR(double *res, int *dat1, int *dat2, int *nr, double *weight, int *e
}
}
}
+*/
-
-void ACCTRAN2(int *dat, int *nr, double *pars, int *node, int *edge, int *nl, double *weight, double *pvec, int *nTips)
+//void ACCTRAN2(int *dat, int *nr, double *pars, int *node, int *edge, int *nl, double *weight, double *pvec, int *nTips)
+void ACCTRAN2(int *dat, int *nr, int *node, int *edge, int *nl, int *nTips)
{
int i;
- for (i=0; i< *nl; i++) {
- if(edge[i]>nTips[0]) fitchT(&dat[(edge[i]-1L) * (*nr)], &dat[(node[i]-1) * (*nr)], nr, pars, weight, &pvec[i]);
+ for (i=0; i< *nl; i++) {
+ // , pars, weight, &pvec[i]
+ if(edge[i]>nTips[0]) fitchT(&dat[(edge[i]-1L) * (*nr)], &dat[(node[i]-1) * (*nr)], nr);
}
}
-void ACCTRAN3(int *dat, int *nr, double *pars, int *node, int *edge, int *nl, double *weight, double *pvec, int *nTips)
+// , int *nTips
+void ACCTRAN3(int *dat, int *nr, double *pars, int *node, int *edge, int *nl, double *weight, double *pvec)
{
int i;
for (i=0; i< *nr; i++)pars[i]=0.0;
diff --git a/src/ml.c b/src/ml.c
index d46bf5c..9769923 100644
--- a/src/ml.c
+++ b/src/ml.c
@@ -32,7 +32,7 @@ const double LOG_SCALE_EPS = -22.18070977791824915926;
// 2^64 = 18446744073709551616
-static double *LL;
+static double *LL, *WEIGHTS;
static int *SCM, *XXX;
void ll_free(){
@@ -71,10 +71,12 @@ void ll_init2(int *data, int *weights, int *nr, int *nTips, int *nc, int *k)
{
int i;
LL = (double *) calloc(*nr * *nc * *k * *nTips, sizeof(double));
+ WEIGHTS = (double *) calloc(*nr, sizeof(double));
XXX = (int *) calloc(*nr * *nTips, sizeof(int));
SCM = (int *) calloc(*nr * *k * *nTips, sizeof(int)); // * 2L
for(i =0; i < (*nr * *k * *nTips); i++) SCM[i] = 0L;
for(i =0; i < (*nr * *nTips); i++) XXX[i] = data[i];
+ for(i =0; i < *nr; i++) WEIGHTS[i] = weights[i];
}
@@ -150,7 +152,9 @@ void rowMinScale(int *dat, int n, int k, int *res){
static R_INLINE void getP(double *eva, double *ev, double *evi, int m, double el, double w, double *result){
int i, j, h;
- double tmp[m], res;
+ double res; //tmp[m],
+ double *tmp;
+ tmp = malloc(m * sizeof(double));
for(i = 0; i < m; i++) tmp[i] = exp(eva[i] * w * el);
for(i = 0; i < m; i++){
for(j = 0; j < m; j++){
@@ -159,6 +163,7 @@ static R_INLINE void getP(double *eva, double *ev, double *evi, int m, double el
result[i+j*m] = res;
}
}
+ free(tmp);
}
@@ -225,7 +230,7 @@ void lll(SEXP dlist, double *eva, double *eve, double *evei, double *el, double
F77_CALL(dgemv)(transa, nr, nc, &one, &ans[ni * rc], nr, bf, &ONE, &zero, TMP, &ONE);
}
-
+/*
// neue Version: keine SEXP (dlist)
void lll0(int *X, double *eva, double *eve, double *evei, double *el, double g, int *nr, int *nc, int *node, int *edge, int nTips, double *contrast, int nco, int n, int *scaleTmp, double *bf, double *TMP, double *ans){
int ni, ei, j, i, rc; // R_len_t i, n = length(node);
@@ -257,7 +262,8 @@ void lll0(int *X, double *eva, double *eve, double *evei, double *el, double g,
scaleMatrix(&ans[ni * rc], nr, nc, scaleTmp);
F77_CALL(dgemv)(transa, nr, nc, &one, &ans[ni * rc], nr, bf, &ONE, &zero, TMP, &ONE);
}
-
+*/
+
// this seems to work perfectly
void lll3(SEXP dlist, double *eva, double *eve, double *evei, double *el, double g, int *nr, int *nc, int *node, int *edge,
int nTips, double *contrast, int nco, int n, int *scaleTmp, double *bf, double *TMP, double *ans, int *SC){
@@ -298,7 +304,7 @@ void lll3(SEXP dlist, double *eva, double *eve, double *evei, double *el, double
F77_CALL(dgemv)(transa, nr, nc, &one, &ans[ni * rc], nr, bf, &ONE, &zero, TMP, &ONE);
}
-
+/*
SEXP PML_NEW2(SEXP EL, SEXP W, SEXP G, SEXP NR, SEXP NC, SEXP K, SEXP eig, SEXP bf, SEXP node, SEXP edge, SEXP NTips, SEXP nco, SEXP contrast, SEXP N){
int nr=INTEGER(NR)[0], nc=INTEGER(NC)[0], k=INTEGER(K)[0], i, indLL;
int nTips = INTEGER(NTips)[0], *SC;
@@ -327,6 +333,7 @@ SEXP PML_NEW2(SEXP EL, SEXP W, SEXP G, SEXP NR, SEXP NC, SEXP K, SEXP eig, SEXP
return TMP;
}
+
// TODO: parallelize
SEXP PML_NEW(SEXP EL, SEXP W, SEXP G, SEXP NR, SEXP NC, SEXP K, SEXP eig, SEXP bf, SEXP node, SEXP edge, SEXP NTips, SEXP nco, SEXP contrast, SEXP N){
int nr=INTEGER(NR)[0], nc=INTEGER(NC)[0], k=INTEGER(K)[0], i, indLL, n=INTEGER(N)[0], ncontr=INTEGER(nco)[0];
@@ -353,9 +360,11 @@ SEXP PML_NEW(SEXP EL, SEXP W, SEXP G, SEXP NR, SEXP NC, SEXP K, SEXP eig, SEXP b
UNPROTECT(1);
return TMP;
}
-
-// in pml.move inside optimRooted
-SEXP PML3(SEXP dlist, SEXP EL, SEXP W, SEXP G, SEXP NR, SEXP NC, SEXP K, SEXP eig, SEXP bf, SEXP node, SEXP edge, SEXP NTips, SEXP nco, SEXP contrast, SEXP N){
+ */
+
+// in pml.move inside optimRooted mit SCM
+// , SEXP W
+SEXP PML3(SEXP dlist, SEXP EL, SEXP G, SEXP NR, SEXP NC, SEXP K, SEXP eig, SEXP bf, SEXP node, SEXP edge, SEXP NTips, SEXP nco, SEXP contrast, SEXP N){
int nr=INTEGER(NR)[0], nc=INTEGER(NC)[0], k=INTEGER(K)[0], i, indLL;
int nTips = INTEGER(NTips)[0], *SC;
double *g=REAL(G), *tmp, logScaleEPS;
@@ -378,8 +387,8 @@ SEXP PML3(SEXP dlist, SEXP EL, SEXP W, SEXP G, SEXP NR, SEXP NC, SEXP K, SEXP ei
return TMP;
}
-
-SEXP PML0(SEXP dlist, SEXP EL, SEXP W, SEXP G, SEXP NR, SEXP NC, SEXP K, SEXP eig, SEXP bf, SEXP node, SEXP edge, SEXP NTips, SEXP nco, SEXP contrast, SEXP N){
+// , SEXP W
+SEXP PML0(SEXP dlist, SEXP EL, SEXP G, SEXP NR, SEXP NC, SEXP K, SEXP eig, SEXP bf, SEXP node, SEXP edge, SEXP NTips, SEXP nco, SEXP contrast, SEXP N){
int nr=INTEGER(NR)[0], nc=INTEGER(NC)[0], k=INTEGER(K)[0], i, indLL;
int nTips = INTEGER(NTips)[0], *SC;
double *g=REAL(G), *tmp, logScaleEPS;
@@ -536,9 +545,9 @@ void goUp(double *dad, int *child, double *contrast, double *P, int nr, int nc,
for(int j=0; j<(nc * nr); j++) dad[j]*=res[j];
}
-
+// , double *w
void updateLLQ(SEXP dlist, int pa, int ch, double *eva, double *eve, double*evei,
- double el, double *w, double *g, int nr,
+ double el, double *g, int nr,
int nc, int ntips, double *contrast, int nco, int k,
double *tmp, double *P){
int i;
@@ -557,9 +566,9 @@ void updateLLQ(SEXP dlist, int pa, int ch, double *eva, double *eve, double*evei
}
-
+// double *w,
void updateLL2(SEXP dlist, int pa, int ch, double *eva, double *eve, double*evei,
- double el, double *w, double *g, int nr,
+ double el, double *g, int nr,
int nc, int ntips, double *contrast, int nco, int k,
double *tmp, double *P){
int i;
@@ -775,9 +784,9 @@ SEXP getM3(SEXP dad, SEXP child, SEXP P, SEXP nr, SEXP nc){
return(RESULT);
}
-
+// SEXP basefreq,
SEXP FS4(SEXP eig, SEXP nc, SEXP el, SEXP w, SEXP g, SEXP X, SEXP dad, SEXP child, SEXP ld, SEXP nr,
- SEXP basefreq, SEXP weight, SEXP f0, SEXP retA, SEXP retB)
+ SEXP weight, SEXP f0, SEXP retA, SEXP retB)
{
SEXP RESULT, EL, P;
double *tmp, *f, *wgt=REAL(weight), edle, ledle, newedle, eps=10, *eva=REAL(VECTOR_ELT(eig,0));
@@ -843,8 +852,8 @@ SEXP FS4(SEXP eig, SEXP nc, SEXP el, SEXP w, SEXP g, SEXP X, SEXP dad, SEXP chil
return (RESULT);
}
-
-SEXP FS5(SEXP eig, SEXP nc, SEXP el, SEXP w, SEXP g, SEXP X, SEXP ld, SEXP nr, SEXP basefreq, SEXP weight, SEXP f0)
+//, SEXP basefreq
+SEXP FS5(SEXP eig, SEXP nc, SEXP el, SEXP w, SEXP g, SEXP X, SEXP ld, SEXP nr, SEXP weight, SEXP f0)
{
SEXP RESULT; // EL, P;
double *tmp, *f, *wgt=REAL(weight), edle, ledle, newedle, eps=10, *eva=REAL(VECTOR_ELT(eig,0));
@@ -1046,7 +1055,7 @@ SEXP optE(SEXP PARENT, SEXP CHILD, SEXP ANC, SEXP eig, SEXP EVI, SEXP EL,
}
}
fs3(eva, nc, oldel, w, g, X, k, nr, weight, f0, res);
- updateLL2(dlist, pa, ch, eva, eve, evei, res[0], w, g, nr,
+ updateLL2(dlist, pa, ch, eva, eve, evei, res[0], g, nr,
nc, ntips, contrast, nco, k, tmp, P);
el[ch-1L] = res[0]; //edgeLengthIndex
if (ch > ntips) loli = ch;
@@ -1119,9 +1128,9 @@ SEXP optQrtt(SEXP PARENT, SEXP CHILD, SEXP eig, SEXP EVI, SEXP EL,
// go up
// if i=2 go down
- if(m==2)updateLLQ(dlist, ch, pa, eva, eve, evei, res[0], w, g, nr,
+ if(m==2)updateLLQ(dlist, ch, pa, eva, eve, evei, res[0], g, nr,
nc, ntips, contrast, nco, k, tmp, P);
- else updateLLQ(dlist, pa, ch, eva, eve, evei, res[0], w, g, nr,
+ else updateLLQ(dlist, pa, ch, eva, eve, evei, res[0], g, nr,
nc, ntips, contrast, nco, k, tmp, P);
el[m] = res[0];
}
diff --git a/src/phangorn.c b/src/phangorn.c
index b081fa0..34035eb 100644
--- a/src/phangorn.c
+++ b/src/phangorn.c
@@ -108,7 +108,9 @@ static R_INLINE void getP00(double *eva, double *ev, double *evi, int m, double
static R_INLINE void getPP(double *eva, double *ev, double *evi, int m, double el, double w, double *result){
int i, j, h;
- double tmp[m];
+// double tmp[m];
+ double *tmp;
+ tmp = malloc(m * sizeof(double));
for(i = 0; i < m; i++) tmp[i] = exp(eva[i] * w * el);
for(i = 0; i < m; i++){
for(j = 0; j < m; j++){
@@ -116,12 +118,15 @@ static R_INLINE void getPP(double *eva, double *ev, double *evi, int m, double e
for(h = 0; h < m; h++) result[i+j*m] += ev[i + h*m] * tmp[h] * evi[h + j*m];
}
}
+ free(tmp);
}
void getdP(double *eva, double *ev, double *evi, int m, double el, double w, double *result){
int i, j, h;
- double tmp[m], res;
+ double res; //tmp[m],
+ double *tmp;
+ tmp = malloc(m * sizeof(double));
for(i = 0; i < m; i++) tmp[i] = (eva[i] * w * el) * exp(eva[i] * w * el);
for(i = 0; i < m; i++){
for(j = 0; j < m; j++){
@@ -130,12 +135,15 @@ void getdP(double *eva, double *ev, double *evi, int m, double el, double w, dou
result[i+j*m] = res;
}
}
+ free(tmp);
}
void getdP2(double *eva, double *ev, double *evi, int m, double el, double w, double *result){
int i, j, h;
- double tmp[m], res;
+ double res; //tmp[m],
+ double *tmp;
+ tmp = malloc(m * sizeof(double));
for(i = 0; i < m; i++) tmp[i] = (eva[i] * w) * exp(eva[i] * w * el);
for(i = 0; i < m; i++){
for(j = 0; j < m; j++){
@@ -144,12 +152,15 @@ void getdP2(double *eva, double *ev, double *evi, int m, double el, double w, do
result[i+j*m] = res;
}
}
+ free(tmp);
}
void getd2P(double *eva, double *ev, double *evi, int m, double el, double w, double *result){
int i, j, h;
- double tmp[m], res;
+ double res; //tmp[m],
+ double *tmp;
+ tmp = malloc(m * sizeof(double));
for(i = 0; i < m; i++) tmp[i] = (eva[i] * w * el) * (eva[i] * w * el) * exp(eva[i] * w * el);
for(i = 0; i < m; i++){
for(j = 0; j < m; j++){
@@ -158,12 +169,15 @@ void getd2P(double *eva, double *ev, double *evi, int m, double el, double w, do
result[i+j*m] = res;
}
}
+ free(tmp);
}
void getd2P2(double *eva, double *ev, double *evi, int m, double el, double w, double *result){
int i, j, h;
- double tmp[m], res;
+ double res; //tmp[m],
+ double *tmp;
+ tmp = malloc(m * sizeof(double));
for(i = 0; i < m; i++) tmp[i] = (eva[i] * w) * (eva[i] * w) * exp(eva[i] * w * el);
for(i = 0; i < m; i++){
for(j = 0; j < m; j++){
@@ -172,6 +186,7 @@ void getd2P2(double *eva, double *ev, double *evi, int m, double el, double w, d
result[i+j*m] = res;
}
}
+ free(tmp);
}
diff --git a/src/phangorn_help.cpp b/src/phangorn_help.cpp
new file mode 100644
index 0000000..c4a299d
--- /dev/null
+++ b/src/phangorn_help.cpp
@@ -0,0 +1,48 @@
+#include <Rcpp.h>
+using namespace Rcpp;
+
+// import: edge matrix, number of tips
+// export: Descendants(x, 1:max(x$edge), "all")
+// [[Rcpp::export]]
+List allDescCPP(IntegerMatrix orig, int nTips) {
+ IntegerVector parent = orig( _, 0);
+ IntegerVector children = orig( _, 1);
+ int m = max(parent);
+ // create list for results
+ std::vector< std::vector<int> > out(m) ;
+ for(int i = 0; i<nTips; i++){
+ out[i].push_back(i + 1L);
+ }
+ std::vector<int> x;
+ std::vector<int> y;
+ for(int i = 0; i<parent.size(); i++){
+ out[parent[i]-1L].push_back(children[i]);
+ if(children[i] > nTips){
+ y = out[children[i] -1L];
+ out[parent[i]-1L].insert( out[parent[i]-1L].end(), y.begin(), y.end() );
+ }
+ }
+ // return the list
+ return wrap(out);
+}
+
+
+// shorter and easier to understand replacement of C function
+// import: edge matrix
+// export: list of children
+// [[Rcpp::export]]
+List allChildrenCPP(IntegerMatrix orig) {
+ IntegerVector parent = orig( _, 0);
+ IntegerVector children = orig( _, 1);
+ int m = max(parent);
+ // create list for results
+ std::vector< std::vector<int> > out(m) ;
+ for(int i = 0; i<parent.size(); i++){
+ out[parent[i]-1L].push_back(children[i]);
+ }
+ return wrap(out);
+}
+
+
+
+
diff --git a/src/sankoff.c b/src/sankoff.c
index 2f5f953..68b0076 100644
--- a/src/sankoff.c
+++ b/src/sankoff.c
@@ -48,7 +48,9 @@ void rowMin2(double *dat, int n, int k, double *res){
}
}
-
+/*
+ * never used
+
void rowMinInt(int *dat, int n, int k, int *res){
int i, h;
int x;
@@ -58,11 +60,13 @@ void rowMinInt(int *dat, int n, int k, int *res){
res[i] = x;
}
}
-
+ */
void sankoff4(double *dat, int n, double *cost, int k, double *result){
int i, j, h;
- double tmp[k], x;
+ double x; //tmp[k],
+ double *tmp;
+ tmp = malloc(k * sizeof(double));
for(i = 0; i < n; i++){
for(j = 0; j < k; j++){
for(h = 0; h< k; h++){tmp[h] = dat[i + h*n] + cost[h + j*k];}
@@ -71,6 +75,7 @@ void sankoff4(double *dat, int n, double *cost, int k, double *result){
result[i+j*n] += x;
}
}
+ free(tmp);
}
diff --git a/src/sprdist.c b/src/sprdist.c
new file mode 100644
index 0000000..082f413
--- /dev/null
+++ b/src/sprdist.c
@@ -0,0 +1,1047 @@
+/*
+ * sprdist.c
+ *
+ * (c) 2016 Leonardo de Oliveira Martins (leomrtns at gmail.com)
+ *
+ *
+ * This code may be distributed under the GNU GPL
+ *
+ */
+
+# define USE_RINTERNALS
+
+#include <Rmath.h>
+#include <math.h>
+#include <R.h>
+#include <Rinternals.h>
+//#include <stdint.h> /* standard integer types (int32_t typedef etc.) [C99]*/
+
+
+#define true 1U /*!< Boolean TRUE */
+#define false 0U /*!< Boolean FALSE */
+
+typedef unsigned char bool;
+typedef struct splitset_struct* splitset;
+typedef struct hungarian_struct* hungarian;
+typedef struct bipartition_struct* bipartition;
+typedef struct bipsize_struct* bipsize;
+
+struct splitset_struct
+{
+ int size, spsize, spr, spr_extra, rf, hdist; /*! \brief spr, extra prunes for spr, rf distances and hdist=assignment cost */
+ int n_g, n_s, n_agree, n_disagree;
+ bipartition *g_split, *s_split, *agree, *disagree;
+ bipartition prune;
+ hungarian h; /* hungarian method for solving the assignment between edges */
+ bool match; /*! \brief do we want to calculate the minimum cost assignment */
+};
+
+struct hungarian_struct
+{
+ int **cost, *col_mate; /*! \brief cost matrix, and col_mate[row] with column match for row */
+ int size, /*! \brief assignment size. Cost is a square matrix, so size should be an overestimate where "missing" nodes are added w/ cost zero */
+ initial_cost, /*! \brief sum of lowest input cost values for each column. The hungarian method rescales them so that minimum per column is zero */
+ final_cost; /*! \brief our final cost is on rescaled cost matrix, therefore to restore the "classical" optimal cost one should sum it with initial_cost */
+ int *unchosen_row, *row_dec, *slack_row, *row_mate, *parent_row, *col_inc, *slack; /* aux vectors */
+};
+
+/*! \brief Bit-string representation of splits. */
+struct bipartition_struct
+{
+ unsigned long long *bs; /*! \brief Representation of a bipartition by a vector of integers (bitstrings). */
+ int n_ones; /*! \brief Counter (number of "one"s) */
+ bipsize n; /*! \brief number of bits (leaves), vector size and mask */
+ int ref_counter; /*! \brief How many times this struct is being referenced */
+};
+
+struct bipsize_struct
+{
+ unsigned long long mask;/*! \brief mask to make sure we consider only active positions (of last bitstring) */
+ int ints, bits, original_size; /*! \brief Vector size and total number of elements (n_ints = n_bits/(8*sizeof(long long)) +1). */
+ int ref_counter; /*! \brief How many times this struct is being referenced */
+};
+
+
+/*! \brief Allocate space for splitset structure (two vectors of bipartitions), for simple comparisons */
+splitset new_splitset (int nleaves, int nsplits);
+/*! \brief free memory allocated for splitset structure */
+void del_splitset (splitset split);
+/*! \brief low level function that does the actual SPR and hdist calculation based on a filled splitset struct */
+int dSPR_topology_lowlevel (splitset split);
+/*! \brief function used by qsort for a vector of bipartitions (from smaller to larger) */
+int compare_splitset_bipartition_increasing (const void *a1, const void *a2);
+
+/* BELOW: low level functions that work with bipartitions */
+void split_create_agreement_list (splitset split);
+void split_remove_agree_edges (splitset split, bipartition *b, int *nb);
+void split_remove_duplicates (bipartition *b, int *nb);
+void split_compress_agreement (splitset split);
+void split_create_disagreement_list (splitset split);
+void split_disagreement_assign_match (splitset split);
+void split_find_small_disagreement (splitset split);
+void split_remove_small_disagreement (splitset split);
+void split_minimize_subtrees (splitset split);
+void split_remove_redundant_bit (splitset split, int id);
+void split_replace_bit (splitset split, int to, int from);
+void split_new_size (splitset split, int size, bool update_bipartitions);
+void split_swap_position (bipartition *b, int i1, int i2);
+
+/* BELOW: Hungarian method for bipartite matching (assignment) */
+hungarian new_hungarian (int size);
+void hungarian_reset (hungarian p);
+void hungarian_update_cost (hungarian p, int row, int col, int cost);
+void del_hungarian (hungarian p);
+void hungarian_solve (hungarian p, int this_size);
+
+/* BELOW: Memory-efficient, fast bipartition comparisions based on 64bit representation of splits */
+bipartition new_bipartition (int size);
+bipsize new_bipsize (int size);
+bipartition new_bipartition_copy_from (const bipartition from);
+bipartition new_bipartition_from_bipsize (bipsize n);
+void del_bipartition (bipartition bip);
+void del_bipsize (bipsize n);
+void bipsize_resize (bipsize n, int nbits);
+void bipartition_initialize (bipartition bip, int position);
+void bipartition_zero (bipartition bip);
+void bipartition_set (bipartition bip, int position);
+void bipartition_set_lowlevel (bipartition bip, int i, int j);
+void bipartition_unset (bipartition bip, int position);
+void bipartition_unset_lowlevel (bipartition bip, int i, int j);
+void bipartition_copy (bipartition to, const bipartition from);
+void bipartition_OR (bipartition result, const bipartition b1, const bipartition b2, bool update_count);
+void bipartition_AND (bipartition result, const bipartition b1, const bipartition b2, bool update_count);
+void bipartition_ANDNOT (bipartition result, const bipartition b1, const bipartition b2, bool update_count);
+void bipartition_XOR (bipartition result, const bipartition b1, const bipartition b2, bool update_count);
+void bipartition_XORNOT (bipartition result, const bipartition b1, const bipartition b2, bool update_count);
+void bipartition_NOT (bipartition result, const bipartition bip);
+void bipartition_count_n_ones (const bipartition bip);
+void bipartition_to_int_vector (const bipartition b, int *id, int vecsize);
+bool bipartition_is_equal (const bipartition b1, const bipartition b2);
+bool bipartition_is_equal_bothsides (const bipartition b1, const bipartition b2);
+bool bipartition_is_larger (const bipartition b1, const bipartition b2);
+void bipartition_flip_to_smaller_set (bipartition bip);
+bool bipartition_is_bit_set (const bipartition bip, int position);
+bool bipartition_contains_bits (const bipartition b1, const bipartition b2);
+// void bipartition_print_to_stdout (const bipartition b1);
+void bipartition_replace_bit_in_vector (bipartition *bvec, int n_b, int to, int from, bool reduce);
+void bipartition_resize_vector (bipartition *bvec, int n_b);
+
+/*! \brief Main SPR calculation function, to be used within R */
+SEXP C_sprdist (SEXP bp1, SEXP bp2, SEXP lt) {
+ int i, j, n_leaves = INTEGER(lt)[0];
+ SEXP result;
+ double *res;
+ splitset split;
+
+ PROTECT(result = allocVector(REALSXP, 4));
+ res = REAL(result);
+
+// if (length(bp1) != length(bp2)) error ("number of bipartitions given to C_sprdist are not the same");
+
+ split = new_splitset (n_leaves, length(bp1));
+ for (i=0; i < split->size; i++) {
+// for (j=0; j < length(VECTOR_ELT (bp1, i)); j++) printf (";;%d ", INTEGER (VECTOR_ELT (bp1, i))[j]);
+ for (j=0; j < length(VECTOR_ELT (bp1, i)); j++) bipartition_set (split->g_split[i], INTEGER (VECTOR_ELT (bp1, i))[j] - 1);
+ for (j=0; j < length(VECTOR_ELT (bp2, i)); j++) bipartition_set (split->s_split[i], INTEGER (VECTOR_ELT (bp2, i))[j] - 1);
+ }
+ dSPR_topology_lowlevel (split);
+ res[0] = split->spr;
+ res[1] = split->spr_extra;
+ res[2] = split->rf;
+ res[3] = split->hdist;
+ del_splitset (split);
+
+ UNPROTECT(1); // result
+ return(result);
+}
+
+/* functions below should not be called outside this scope */
+
+splitset
+new_splitset (int nleaves, int nsplits)
+{
+ splitset split;
+ int i;
+
+ split = (splitset) malloc (sizeof (struct splitset_struct));
+ split->n_g = split->n_s = split->size = nsplits;
+ split->n_agree = split->n_disagree = 0;
+ split->prune = NULL;
+ split->match = true; /* do we want to calculate the assignment matching cost (using hungarian() )? */
+ split->spr = split->spr_extra = split->rf = split->hdist = 0;
+
+ split->g_split = (bipartition*) malloc (split->size * sizeof (bipartition));
+ split->s_split = (bipartition*) malloc (split->size * sizeof (bipartition));
+ split->g_split[0] = new_bipartition (nleaves);
+ split->s_split[0] = new_bipartition (nleaves);
+ for (i = 1; i < split->size; i++) {
+ split->g_split[i] = new_bipartition_from_bipsize (split->g_split[0]->n); /* use same bipsize */
+ split->s_split[i] = new_bipartition_from_bipsize (split->s_split[0]->n);
+ }
+
+
+ split->agree = (bipartition*) malloc (split->size * sizeof (bipartition));
+ split->disagree = (bipartition*) malloc (split->size * split->size * sizeof (bipartition));
+ split->agree[0] = new_bipartition (nleaves); // this bipsize will be recycled below
+ split->disagree[0] = new_bipartition (nleaves);
+ for (i = 1; i < split->size; i++) split->agree[i] = new_bipartition_from_bipsize (split->agree[0]->n);
+ for (i = 1; i < split->size * split->size; i++) split->disagree[i] = new_bipartition_from_bipsize (split->disagree[0]->n);
+ split->prune = new_bipartition_from_bipsize (split->disagree[0]->n);
+
+ split->h = new_hungarian (split->size);
+
+ return split;
+}
+
+void
+del_splitset (splitset split)
+{
+ int i;
+ if (!split) return;
+
+ del_bipartition (split->prune);
+ if (split->disagree) {
+ for (i = split->size * split->size - 1; i >= 0; i--) del_bipartition (split->disagree[i]);
+ free (split->disagree);
+ }
+ if (split->agree) {
+ for (i = split->size - 1; i >= 0; i--) del_bipartition (split->agree[i]);
+ free (split->agree);
+ }
+ if (split->g_split) {
+ for (i = split->size - 1; i >= 0; i--) del_bipartition (split->g_split[i]);
+ free (split->g_split);
+ }
+ del_hungarian (split->h);
+ free (split);
+}
+
+int
+compare_splitset_bipartition_increasing (const void *a1, const void *a2)
+{ /* similar to bipartition_is_larger() */
+ bipartition *b1 = (bipartition *) a1;
+ bipartition *b2 = (bipartition *) a2;
+ int i;
+
+ if ((*b1)->n_ones > (*b2)->n_ones) return 1;
+ if ((*b1)->n_ones < (*b2)->n_ones) return -1;
+
+ for (i = (*b1)->n->ints - 1; (i >= 0) && ((*b1)->bs[i] == (*b2)->bs[i]); i--); /* find position of distinct bipartition elem*/
+ if (i < 0) return 0; /* identical bipartitions */
+ if ((*b1)->bs[i] > (*b2)->bs[i]) return 1;
+ else return -1;
+}
+
+
+int
+dSPR_topology_lowlevel (splitset split)
+{
+ int i = 0, mismatch = -1;
+
+ for (i=0; i < split->size; i++) {
+ bipartition_flip_to_smaller_set (split->g_split[i]);
+ bipartition_flip_to_smaller_set (split->s_split[i]);
+ }
+ qsort (split->g_split, split->size, sizeof (bipartition), compare_splitset_bipartition_increasing);
+ qsort (split->s_split, split->size, sizeof (bipartition), compare_splitset_bipartition_increasing);
+ //for (i = 0; i < split->n_g; i++) bipartition_print_to_stdout (split->g_split[i]); printf ("G ::DEBUG 0 ::\n");
+ //for (i = 0; i < split->n_s; i++) bipartition_print_to_stdout (split->s_split[i]); printf ("S\n");
+
+ i++; /* to trick -Werror, since we don't use it unless for debug */
+ while (mismatch) {
+ split_create_agreement_list (split); // vector of identical bipartitions
+ split_compress_agreement (split); // iterative replacement of cherry by new leaf
+
+// for (i = 0; i < split->n_g; i++) bipartition_print_to_stdout (split->g_split[i]); printf ("G ::DEBUG::\n");
+// for (i = 0; i < split->n_s; i++) bipartition_print_to_stdout (split->s_split[i]); printf ("S\n");
+// for (i = 0; i < split->n_agree; i++) bipartition_print_to_stdout (split->agree[i]); printf ("A\n");
+
+ if (mismatch == -1) split->rf = split->n_g + split->n_s;
+ mismatch = (split->n_g > 0) && (split->n_s > 0); // all edges were in agreement
+ if (!mismatch) return split->spr;
+
+ split_create_disagreement_list (split); // vector of smallest disagreements
+ split_disagreement_assign_match (split); /* assignment matching between edges using hungarian method (split->hdist after first time) */
+
+ split_remove_duplicates (split->disagree, &(split->n_disagree)); // some elements are equal; this function also qsorts
+ split_find_small_disagreement (split); // could also be one leaf only
+
+ //for (i = 0; i < split->n_disagree; i++) { bipartition_print_to_stdout (split->disagree[i]); printf ("\n"); }
+ //printf ("{%d} prune: ", split->n_disagree); bipartition_print_to_stdout (split->prune); printf ("\n");
+
+ split->spr++;
+ split_remove_small_disagreement (split);
+
+ split_minimize_subtrees (split);
+ mismatch = (split->n_g > 0) && (split->n_s > 0); // all edges were in agreement
+ }
+ return split->spr;
+}
+
+void
+split_create_agreement_list (splitset split)
+{
+ int s, g;
+ for (g = 0; g < split->n_g; g++) for (s = 0; s < split->n_s; s++)
+ if (bipartition_is_equal (split->g_split[g], split->s_split[s])) {
+ bipartition_copy (split->agree[split->n_agree++], split->g_split[g]);
+ split->n_g--; split_swap_position (split->g_split, g, split->n_g); /* if we don't swap them, we lose ref to "old" value on g_split[] */
+ split->n_s--; split_swap_position (split->s_split, s, split->n_s);
+ g--; s = split->n_s; /* pretend loop finished, examine again with new values */
+ }
+ split_remove_agree_edges (split, split->g_split, &(split->n_g));
+ split_remove_agree_edges (split, split->s_split, &(split->n_s));
+}
+
+void
+split_remove_agree_edges (splitset split, bipartition *b, int *nb)
+{
+ int i, a;
+ for (i = 0; i < (*nb); i++) for (a = 0; a < split->n_agree; a++)
+ if (bipartition_is_equal (b[i], split->agree[a])) {
+ (*nb)--;
+ split_swap_position (b, i, (*nb));
+ i--;
+ a = split->n_agree; /* loop again over new value */
+ }
+}
+
+void
+split_remove_duplicates (bipartition *b, int *nb)
+{
+ int i, j;
+ bipartition pivot;
+ if ((*nb) < 2) return; /* only if we have a vector with > 1 element */
+ qsort (b, (*nb), sizeof (bipartition), compare_splitset_bipartition_increasing);
+
+ for (i = (*nb) - 1; i >= 1; i--)
+ if (bipartition_is_equal (b[i], b[i-1])) {
+ pivot = b[i]; /* do not lose a pointer to this element */
+ for (j = i; j < (*nb)-1; j++) b[j] = b[j+1];
+ b[j] = pivot; /* j = (*nb) - 1, which will become obsolete through next line --> (*nb)-- */
+ (*nb)--;
+ }
+}
+
+void
+split_compress_agreement (splitset split)
+{
+ int i, j, pair[2];
+
+ for (i = 0; i < split->n_agree; i++) if (split->agree[i]->n_ones == 2) { /* cherry in common, can be represented by just one leaf */
+ bipartition_to_int_vector (split->agree[i], pair, 2);
+ split_remove_redundant_bit (split, pair[1]);
+ split_new_size (split,split->agree[0]->n->bits - 1, false); /* false = do not recalculate every bipartition's last elem */
+ bipartition_resize_vector (split->agree, split->n_agree);
+ for (j = 0; j < split->n_agree; j++) { /* minimize subtree size and remove single leaves */
+ bipartition_flip_to_smaller_set (split->agree[j]); /* agree only */
+ if (split->agree[j]->n_ones < 2) split_swap_position (split->agree, j--, --split->n_agree);
+ }
+ i = -1; /* redo all iterations, with new info (agree[] will be smaller) */
+ }
+ bipartition_resize_vector (split->g_split, split->n_g);
+ bipartition_resize_vector (split->s_split, split->n_s);
+}
+
+void
+split_create_disagreement_list (splitset split)
+{
+ int g, s;
+
+ for (g = 0; g < split->n_g; g++) for (s = 0; s < split->n_s; s++) {
+ bipartition_XOR (split->disagree[g * split->n_s + s], split->g_split[g], split->s_split[s], true); /* true means to calculate n_ones */
+ bipartition_flip_to_smaller_set (split->disagree[g * split->n_s + s]);
+ }
+ split->n_disagree = split->n_g * split->n_s;
+}
+
+void
+split_disagreement_assign_match (splitset split)
+{ /* also calculates split->hdist */
+ int g, s, max_n, sum = 0;
+
+ if (split->n_g > split->n_s) max_n = split->n_g;
+ else max_n = split->n_s;
+ if (max_n < 2) return;
+
+ hungarian_reset (split->h);
+ for (g = 0; g < split->n_g; g++) for (s = 0; s < split->n_s; s++)
+ hungarian_update_cost (split->h, g, s, split->disagree[g * split->n_s + s]->n_ones);
+ hungarian_solve (split->h, max_n);
+ /* now split->h->col_mate will have the pairs */
+ /* if we do the matching below it becomes much faster, but we may miss the best prune subtrees in a few cases (do not compromise the algo) */
+ split->n_disagree = 0;
+ for (g = 0; g < max_n; g++) if ((g < split->n_g) && ( split->h->col_mate[g] < split->n_s)) { /* some matchings might be to dummy edges */
+ bipartition_XOR (split->disagree[split->n_disagree], split->g_split[g], split->s_split[split->h->col_mate[g]], true); /* true means to calculate n_ones */
+ bipartition_flip_to_smaller_set (split->disagree[split->n_disagree++]);
+ sum += split->disagree[split->n_disagree-1]->n_ones;
+ }
+ if (split->match) { split->hdist = split->h->final_cost+split->h->initial_cost; split->match = false; }
+}
+
+void
+split_find_small_disagreement (splitset split)
+{
+ bipartition dis;
+ int a, d;
+
+ bipartition_copy (split->prune, split->disagree[0]); /* smallest, in case we don't find a better one in loop below */
+ if (split->prune->n_ones < 2) return;
+
+ dis = new_bipartition_from_bipsize (split->disagree[0]->n);
+ for (d = 0; d < split->n_disagree; d++) for (a = 0; a < split->n_agree; a++) {
+ if ((split->disagree[d]->n_ones == split->agree[a]->n_ones) ||
+ (split->disagree[d]->n_ones == (split->agree[a]->n->bits - split->agree[a]->n_ones))) {
+ bipartition_XOR (dis, split->disagree[d], split->agree[a], true);
+ if (!dis->n_ones) { bipartition_copy (split->prune, split->disagree[d]); d = split->n_disagree; a = split->n_agree; }
+ else if (dis->n_ones == dis->n->bits) { bipartition_NOT (split->prune, split->disagree[d]); d = split->n_disagree; a = split->n_agree; }
+ }
+ }
+ /* check if prune nodes are all on same side of a tree or if they are actually two SPRs (one from each tree) */
+ for (d = 0; d < split->n_g; d++) {
+ if (!bipartition_contains_bits (split->g_split[d], split->prune)) {
+ bipartition_NOT (dis, split->g_split[d]);
+ if (!bipartition_contains_bits (dis, split->prune)) { split->spr_extra++; d = split->n_g; }
+ }
+ }
+ del_bipartition (dis);
+}
+
+void
+split_remove_small_disagreement (splitset split)
+{
+ int *index, i, j = split->prune->n_ones - 1, k = 0, size = split->agree[0]->n->bits;
+
+ index = (int*) malloc (split->prune->n_ones * sizeof (int));
+ bipartition_to_int_vector (split->prune, index, split->prune->n_ones);
+
+ for (i = size - 1; i >= (size - split->prune->n_ones); i--) {
+ if (index[k] >= (size - split->prune->n_ones)) i = -1;
+ else {
+ if (i == index[j]) j--;
+ else split_replace_bit (split, index[k++], i);
+ }
+ }
+
+ split_new_size (split,size - split->prune->n_ones, true);
+ if (index) free (index);
+}
+
+void
+split_minimize_subtrees (splitset split)
+{
+ int i;
+
+ for (i = 0; i < split->n_s; i++) {
+ bipartition_flip_to_smaller_set (split->s_split[i]);
+ if (split->s_split[i]->n_ones < 2) { split->n_s--; split_swap_position (split->s_split, i, split->n_s); i--; }
+ }
+ for (i = 0; i < split->n_g; i++) {
+ bipartition_flip_to_smaller_set (split->g_split[i]);
+ if (split->g_split[i]->n_ones < 2) { split->n_g--; split_swap_position (split->g_split, i, split->n_g); i--; }
+ }
+ for (i = 0; i < split->n_agree; i++) {
+ bipartition_flip_to_smaller_set (split->agree[i]);
+ if (split->agree[i]->n_ones < 2) { split->n_agree--; split_swap_position (split->agree, i, split->n_agree); i--; }
+ }
+}
+
+void
+split_remove_redundant_bit (splitset split, int id)
+{
+ int last = split->agree[0]->n->bits-1;
+ if (id < last) split_replace_bit (split, id, last);
+}
+
+void
+split_replace_bit (splitset split, int to, int from)
+{
+ if (from <= to) return;
+ /*not needed for disagree[] */
+ bipartition_replace_bit_in_vector (split->agree, split->n_agree, to, from, true);
+ bipartition_replace_bit_in_vector (split->g_split, split->n_g, to, from, true);
+ bipartition_replace_bit_in_vector (split->s_split, split->n_s, to, from, true);
+}
+
+void
+split_new_size (splitset split, int size, bool update_bipartitions)
+{
+ bipsize_resize (split->g_split[0]->n, size);
+ bipsize_resize (split->s_split[0]->n, size);
+ bipsize_resize (split->agree[0]->n, size);
+ bipsize_resize (split->disagree[0]->n, size);
+ if (update_bipartitions) {
+ bipartition_resize_vector (split->g_split, split->n_g);
+ bipartition_resize_vector (split->s_split, split->n_s);
+ bipartition_resize_vector (split->agree, split->n_agree);
+ }
+}
+
+void
+split_swap_position (bipartition *b, int i1, int i2)
+{
+ bipartition pivot = b[i1];
+ b[i1] = b[i2];
+ b[i2] = pivot;
+}
+
+
+/* The hungarian method below is copied from http://www.informatik.uni-freiburg.de/~stachnis/misc.html
+ * The (edited) original message follows:
+ *
+ ** libhungarian by Cyrill Stachniss, 2004 Solving the Minimum Assignment Problem using the
+ ** Hungarian Method. ** This file may be freely copied and distributed! **
+ **
+ ** Parts of the used code was originally provided by the "Stanford GraphGase", but I made changes to this code.
+ ** As asked by the copyright node of the "Stanford GraphGase", I hereby proclaim that this file are *NOT* part of the
+ ** "Stanford GraphGase" distrubition! */
+
+void
+hungarian_reset (hungarian p)
+{
+ int i, j;
+
+ for (i = 0; i < p->size; i++) {
+ p->col_mate[i] = p->unchosen_row[i] = p->row_dec[i] = p->slack_row[i] = p->row_mate[i] = p->parent_row[i] = p->col_inc[i] = p->slack[i] = 0;
+ for (j = 0; j < p->size; j++) p->cost[i][j] = 0;
+ }
+ p->final_cost = 0;
+}
+
+hungarian
+new_hungarian (int size)
+{
+ int i;
+ hungarian p;
+
+ p = (hungarian) malloc (sizeof (struct hungarian_struct));
+ p->size = size; /* n_rows = n_columns; if it's not, fill with zeroes (no cost) */
+ p->cost = (int**) malloc (size * sizeof (int*));
+ for (i = 0; i < p->size; i++)
+ p->cost[i] = (int*) malloc (size * sizeof (int));
+ /* edges would be assignment_matrix[ i * ncols + col_mate[i] ] = true; and other elems "false" (but we don't use the matrix notation) */
+ p->col_mate = (int*) malloc (size * sizeof (int)); /* for a given row node, col_mate[row] is the assigned col node */
+ p->unchosen_row = (int*) malloc (size * sizeof (int));
+ p->row_dec = (int*) malloc (size * sizeof (int));
+ p->slack_row = (int*) malloc (size * sizeof (int));
+ p->row_mate = (int*) malloc (size * sizeof (int));
+ p->parent_row = (int*) malloc (size * sizeof (int));
+ p->col_inc = (int*) malloc (size * sizeof (int));
+ p->slack = (int*) malloc (size * sizeof (int));
+
+ hungarian_reset (p);
+ return p;
+}
+
+void
+hungarian_update_cost (hungarian p, int row, int col, int cost)
+{
+ if (row >= p->size) return;
+ if (col >= p->size) return;
+ p->cost[row][col] = cost;
+}
+
+void
+del_hungarian (hungarian p)
+{
+ int i;
+ if (!p) return;
+ if (p->cost) {
+ for (i = p->size - 1; i >= 0; i--) if (p->cost[i]) free (p->cost[i]);
+ free (p->cost);
+ }
+ free (p->col_mate); /* this is the important one, with i assigned to col_mate[i] */
+ free (p->slack);
+ free (p->col_inc);
+ free (p->parent_row);
+ free (p->row_mate);
+ free (p->slack_row);
+ free (p->row_dec);
+ free (p->unchosen_row);
+ free (p);
+}
+
+void
+hungarian_solve (hungarian p, int this_size)
+{
+ int i, j, nrows = this_size, ncols = this_size, k, l, s, t, q, unmatched;
+ p->final_cost = p->initial_cost = 0;
+
+ if (this_size > p->size) { p->final_cost = -1; return; } /* we don't call biomcmc_error(), but it *is* an error! */
+
+ for (l = 0; l < ncols; l++) { // Begin subtract column minima in order to start with lots of zeroes 12
+ s = p->cost[0][l];
+ for (k = 1; k < nrows; k++) if (p->cost[k][l] < s) s = p->cost[k][l];
+ p->initial_cost += s; /* this should be added to final_cost to have classical assignment cost; here we distinguish them */
+ if (s!=0) for (k = 0; k < nrows; k++) p->cost[k][l] -= s;
+ } // End subtract column minima in order to start with lots of zeroes 12
+
+ // Begin initial state 16
+ t=0;
+ for (l = 0; l < ncols; l++) { // n => num_cols
+ p->row_mate[l]= -1;
+ p->parent_row[l]= -1;
+ p->col_inc[l]=0;
+ p->slack[l]= 0x7FFFFFFF;
+ }
+ for (k = 0; k < nrows; k++) { // m => num_rows
+ s = p->cost[k][0];
+ for (l = 1; l < ncols; l++) if (p->cost[k][l] < s) s = p->cost[k][l];
+ p->row_dec[k]=s;
+ for (l = 0; l < ncols; l++) if ((s==p->cost[k][l]) && (p->row_mate[l] < 0)) {
+ p->col_mate[k] = l;
+ p->row_mate[l] = k; // fprintf(stderr, "matching col %d==row %d\n",l,k);
+ goto row_done;
+ }
+ p->col_mate[k] = -1; // fprintf(stderr, "node %d: unmatched row %d\n",t,k);
+ p->unchosen_row[t++] = k;
+row_done:
+ ;
+ }
+ // End initial state 16
+
+ // Begin Hungarian algorithm 18
+ if (t==0) goto done;
+ unmatched=t;
+ while (1) {
+ q=0; // fprintf(stderr, "Matched %d rows.\n",m-t);
+ while (1) {
+ while (q<t) {
+ { // Begin explore node q of the forest 19
+ k = p->unchosen_row[q];
+ s=p->row_dec[k];
+ for (l=0;l<ncols;l++) if (p->slack[l]) {
+ int del;
+ del = p->cost[k][l] - s + p->col_inc[l];
+ if (del < p->slack[l]) {
+ if (del==0) {
+ if (p->row_mate[l]<0) goto breakthru;
+ p->slack[l]=0;
+ p->parent_row[l]=k; // fprintf(stderr, "node %d: row %d==col %d--row %d\n", t,row_mate[l],l,k);
+ p->unchosen_row[t++]=p->row_mate[l];
+ }
+ else { p->slack[l]=del; p->slack_row[l]=k; }
+ }
+ }
+ } // End explore node q of the forest 19
+ q++;
+ }
+
+ // Begin introduce a new zero into the matrix 21
+ s = 0x7FFFFFFF;
+ for (l = 0;l < ncols; l++) if (p->slack[l] && p->slack[l] < s) s = p->slack[l];
+ for (q = 0; q < t; q++) p->row_dec[ p->unchosen_row[q] ] += s;
+ for (l = 0; l < ncols; l++) if (p->slack[l]) {
+ p->slack[l]-=s;
+ if (p->slack[l]==0) { // Begin look at a new zero 22
+ k = p->slack_row[l]; // fprintf(stderr, "Decreasing uncovered elements by %d produces zero at [%d,%d]\n", s,k,l);
+ if (p->row_mate[l]<0) {
+ for (j=l+1;j<ncols;j++) if (p->slack[j]==0) p->col_inc[j]+=s;
+ goto breakthru;
+ }
+ else {
+ p->parent_row[l]=k; // fprintf(stderr, "node %d: row %d==col %d--row %d\n",t,row_mate[l],l,k);
+ p->unchosen_row[t++]=p->row_mate[l];
+ }
+ } // End look at a new zero 22
+
+ }
+ else p->col_inc[l]+=s;
+ // End introduce a new zero into the matrix 21
+ }
+breakthru:
+ // fprintf(stderr, "Breakthrough at node %d of %d!\n",q,t);
+ while (1) { // Begin update the matching 20
+ j=p->col_mate[k];
+ p->col_mate[k]=l;
+ p->row_mate[l]=k; // fprintf(stderr, "rematching col %d==row %d\n",l,k);
+ if (j<0) break;
+ k=p->parent_row[j];
+ l=j;
+ } // End update the matching 20
+ if (--unmatched==0) goto done;
+ // Begin get ready for another stage 17
+ t=0;
+ for (l=0;l<ncols;l++) {
+ p->parent_row[l]= -1;
+ p->slack[l]=0x7FFFFFFF;
+ }
+ for (k=0;k<nrows;k++) if (p->col_mate[k]<0) p->unchosen_row[t++]=k; // fprintf(stderr, "node %d: unmatched row %d\n",t,k);
+ // End get ready for another stage 17
+ }
+done:
+ // Begin doublecheck the solution 23
+ for (k = 0; k < nrows; k++) for (l=0;l<ncols;l++) if (p->cost[k][l] < p->row_dec[k] - p->col_inc[l]) { p->final_cost = -1; return;} //printf ("\n**\n");
+ for (k = 0; k < nrows; k++) {
+ l=p->col_mate[k];
+ if ((l < 0) || (p->cost[k][l] != p->row_dec[k] - p->col_inc[l])) { p->final_cost = -1; return; }
+ }
+ k=0;
+ for (l=0;l<ncols;l++) if (p->col_inc[l]) k++;
+ if (k>nrows) { p->final_cost = -1; return; }
+ // End doublecheck the solution 23
+ // End Hungarian algorithm 18
+
+ for (k = 0; k < nrows; ++k) for (l = 0; l < ncols; ++l) p->cost[k][l] = p->cost[k][l] - p->row_dec[k] + p->col_inc[l];
+ for (i = 0; i < nrows; i++) p->final_cost += p->row_dec[i];
+ for (i = 0; i < ncols; i++) p->final_cost -= p->col_inc[i]; // fprintf(stderr, "Cost is %d\n",cost);
+}
+ /* BELOW are the bipartition functions (memory-efficient storage of bipartitions on 64 bits */
+
+int BitStringSize = 8 * sizeof (unsigned long long);
+
+bipartition
+new_bipartition (int size)
+{
+ bipartition bip;
+ int i;
+
+ bip = (bipartition) malloc (sizeof (struct bipartition_struct));
+ bip->n = new_bipsize (size);
+ bip->n_ones = 0;
+ bip->ref_counter = 1;
+
+ bip->bs = (unsigned long long*) malloc (bip->n->ints * sizeof (unsigned long long));
+ for (i=0; i < bip->n->ints; i++) bip->bs[i] = 0LL;
+
+ return bip;
+}
+
+bipsize
+new_bipsize (int size)
+{
+ bipsize n;
+ int i;
+
+ n = (bipsize) malloc (sizeof (struct bipsize_struct));
+ n->bits = n->original_size = size;
+ n->ref_counter = 1;
+ n->ints = size/BitStringSize + 1;
+ n->mask = 0LL;
+ for (i=0; i < n->bits%BitStringSize; i++) n->mask |= (1LL << i); /* disregard other bits */
+
+ return n;
+}
+
+bipartition
+new_bipartition_copy_from (const bipartition from)
+{
+ bipartition bip;
+ int i;
+
+ bip = (bipartition) malloc (sizeof (struct bipartition_struct));
+ bip->n = new_bipsize (from->n->bits);
+ bip->n_ones = from->n_ones;
+ bip->ref_counter = 1;
+
+ bip->bs = (unsigned long long*) malloc (bip->n->ints * sizeof (unsigned long long));
+ for (i=0; i < bip->n->ints; i++) bip->bs[i] = from->bs[i];
+
+ return bip;
+}
+
+bipartition
+new_bipartition_from_bipsize (bipsize n)
+{
+ bipartition bip;
+ int i;
+
+ bip = (bipartition) malloc (sizeof (struct bipartition_struct));
+ bip->n = n;
+ bip->n->ref_counter++;
+ bip->n_ones = 0;
+ bip->ref_counter = 1;
+
+ bip->bs = (unsigned long long*) malloc (bip->n->ints * sizeof (unsigned long long));
+ for (i=0; i < bip->n->ints; i++) bip->bs[i] = 0LL;
+
+ return bip;
+}
+
+void
+del_bipartition (bipartition bip)
+{
+ if (bip) {
+ if (--bip->ref_counter) return;
+ if (bip->bs) free (bip->bs);
+ del_bipsize (bip->n);
+ free (bip);
+ }
+}
+
+void
+del_bipsize (bipsize n)
+{
+ if (n) {
+ if (--n->ref_counter) return;
+ free (n);
+ }
+}
+
+void
+bipsize_resize (bipsize n, int nbits)
+{
+ int i;
+ n->bits = nbits;
+ n->ints = nbits/BitStringSize + 1; // might be smaller than original bs size
+ n->mask = 0LL;
+ for (i=0; i < nbits%BitStringSize; i++) n->mask |= (1LL << i); /* disregard other bits */
+}
+
+void
+bipartition_initialize (bipartition bip, int position)
+{
+ int i, j;
+ for (i=0; i < bip->n->ints; i++) bip->bs[i] = 0LL;
+ j = position%BitStringSize;
+ i = position/BitStringSize;
+
+ bip->bs[i] = (1LL << j);
+ bip->n_ones = 1;
+}
+
+void
+bipartition_zero (bipartition bip)
+{
+ int i;
+ for (i=0; i < bip->n->ints; i++) bip->bs[i] = 0LL;
+ bip->n_ones = 0;
+}
+
+void
+bipartition_set (bipartition bip, int position)
+{
+ bipartition_set_lowlevel (bip, position/BitStringSize, position%BitStringSize);
+}
+
+void
+bipartition_set_lowlevel (bipartition bip, int i, int j)
+{
+ if (bip->bs[i] & (1LL << j)) return; // bit already set
+ bip->bs[i] |= (1LL << j);
+ bip->n_ones++; /* doesn't work if we reduce space later (check replace_int_in_vector() ) */
+}
+
+void
+bipartition_unset (bipartition bip, int position)
+{
+ bipartition_unset_lowlevel (bip, position/BitStringSize, position%BitStringSize);
+}
+
+void
+bipartition_unset_lowlevel (bipartition bip, int i, int j)
+{
+ if (!(bip->bs[i] & (1LL << j))) return; // bit already unset
+ bip->bs[i] &= ~(1LL << j);
+ bip->n_ones--;
+}
+
+void
+bipartition_copy (bipartition to, const bipartition from)
+{
+ int i;
+ for (i=0; i < to->n->ints; i++) to->bs[i] = from->bs[i];
+ to->n_ones = from->n_ones;
+}
+
+void
+bipartition_OR (bipartition result, const bipartition b1, const bipartition b2, bool update_count)
+{
+ int i;
+ for (i=0; i < result->n->ints; i++) result->bs[i] = b1->bs[i] | b2->bs[i];
+ result->bs[i-1] &= b1->n->mask; /* do not change last bits (do not belong to bipartition) */
+ if (update_count) bipartition_count_n_ones (result);
+ else result->n_ones = b1->n_ones + b2->n_ones; // works on topologies where b1 and b2 are disjoint
+}
+
+void
+bipartition_AND (bipartition result, const bipartition b1, const bipartition b2, bool update_count)
+{
+ int i;
+ for (i=0; i < result->n->ints; i++) result->bs[i] = b1->bs[i] & b2->bs[i];
+ result->bs[i-1] &= b1->n->mask; /* do not change last bits (do not belong to bipartition) */
+ if (update_count) bipartition_count_n_ones (result);
+ else result->n_ones = 0;// update_count = false should be used only when you don't care about this value (temp var)
+}
+
+void
+bipartition_ANDNOT (bipartition result, const bipartition b1, const bipartition b2, bool update_count)
+{
+ int i;
+ for (i=0; i < result->n->ints; i++) result->bs[i] = b1->bs[i] & (~b2->bs[i]);
+ result->bs[i-1] &= b1->n->mask; /* do not change last bits (do not belong to bipartition) */
+ if (update_count) bipartition_count_n_ones (result);
+ else result->n_ones = 0;// update_count = false should be used only when you don't care about this value (temp var)
+}
+
+void
+bipartition_XOR (bipartition result, const bipartition b1, const bipartition b2, bool update_count)
+{
+ int i;
+ for (i=0; i < result->n->ints; i++) result->bs[i] = b1->bs[i] ^ b2->bs[i];
+ result->bs[i-1] &= b1->n->mask; /* do not change last bits (do not belong to bipartition) */
+ if (update_count) bipartition_count_n_ones (result);
+ else result->n_ones = 0;// update_count = false should be used only when you don't care about this value (temp var)
+}
+
+void
+bipartition_XORNOT (bipartition result, const bipartition b1, const bipartition b2, bool update_count)
+{ /* equivalent to XOR followed by NOT */
+ int i;
+ for (i=0; i < result->n->ints; i++) result->bs[i] = b1->bs[i] ^ (~b2->bs[i]);
+ result->bs[i-1] &= b1->n->mask; /* do not change last bits (do not belong to bipartition) */
+ if (update_count) bipartition_count_n_ones (result);
+ else result->n_ones = 0;// update_count = false should be used only when you don't care about this value (temp var)
+}
+
+void
+bipartition_NOT (bipartition result, const bipartition bip)
+{
+ int i;
+ for (i=0; i < result->n->ints; i++) result->bs[i] = ~bip->bs[i];
+ result->bs[i-1] &= bip->n->mask; /* do not invert last bits (do not belong to bipartition) */
+ result->n_ones = bip->n->bits - bip->n_ones;
+}
+
+void
+bipartition_count_n_ones (const bipartition bip)
+{
+ int i;
+ unsigned long long j;
+ bip->n_ones = 0;
+/* // Naive approach
+ for (i=0; i < bip->n_ints - 1; i++) for (j=0; j < BitStringSize; j++) bip->n_ones += ((bip->bs[i] >> j) & 1LL);
+ for (j=0; j < bip->n_bits%BitStringSize; j++) bip->n_ones += ((bip->bs[i] >> j) & 1LL);
+ */
+ // clear the least significant bit set per iteration (Peter Wegner in CACM 3 (1960), 322, mentioned in K&R)
+ for (i=0; i < bip->n->ints; i++) for (j = bip->bs[i]; j; bip->n_ones++) j &= j - 1LL;
+}
+
+bool
+bipartition_is_equal (const bipartition b1, const bipartition b2)
+{
+ int i;
+ if (b1->n_ones != b2->n_ones) return false;
+ if (b1->n->ints != b2->n->ints) return false;
+ for (i=0; i < b1->n->ints - 1; i++) if (b1->bs[i] != b2->bs[i]) return false;
+ b1->bs[i] &= b1->n->mask; b2->bs[i] &= b2->n->mask; /* apply mask before comparing last elems */
+ if (b1->bs[i] != b2->bs[i]) return false;
+ return true;
+}
+
+bool
+bipartition_is_equal_bothsides (const bipartition b1, const bipartition b2)
+{
+ int i;
+ bool equal = true;
+ for (i=0; (i < b1->n->ints - 1) && (equal); i++) if (b1->bs[i] != b2->bs[i]) equal = false;
+ if ((equal) && ((b1->bs[i] & b1->n->mask) != (b2->bs[i] & b2->n->mask))) equal = false;
+ if (equal) return true; /* the biparitions are already the same, without flipping the bits */
+ /* now we compare one bipartition with the complement of the other */
+ for (i=0; (i < b1->n->ints - 1); i++) if (b1->bs[i] != ~b2->bs[i]) return false;
+ if ((b1->bs[i] & b1->n->mask) != ((~b2->bs[i]) & b2->n->mask)) return false;
+ return true; /* they are the exact complement of one another */
+}
+
+bool
+bipartition_is_larger (const bipartition b1, const bipartition b2)
+{
+ int i;
+ if (b1->n_ones > b2->n_ones) return true;
+ if (b1->n_ones < b2->n_ones) return false;
+
+ for (i = b1->n->ints - 1; (i >= 0) && (b1->bs[i] == b2->bs[i]); i--); /* find position of distinct bipartition elem*/
+
+ if (i < 0) return false; /* identical bipartitions */
+ if (b1->bs[i] > b2->bs[i]) return true;
+ else return false;
+}
+
+void
+bipartition_flip_to_smaller_set (bipartition bip)
+{
+ int i = bip->n->ints - 1; /* most significant position -- consistent with is_larger() above, using OLD algo below */
+ if ((2 * bip->n_ones) < bip->n->bits) return; /* it is already the smaller set */
+ /* OLD always x is different from ~x, so we just look at last element ("largest digits of number") */
+ // if (((2 * bip->n_ones) == bip->n->bits) && (bip->bs[i] < (bip->n->mask & ~bip->bs[i]))) return;
+ /* NEW: resolve ties by always showing the same "side" of bipartition, that is, the one having an arbitrary leaf (first one, in our case) */
+ if (((2 * bip->n_ones) == bip->n->bits) && (bip->bs[0] & 1LL)) return;
+
+ for (i=0; i < bip->n->ints; i++) bip->bs[i] = ~bip->bs[i]; /* like bipartition_NOT() */
+ bip->bs[i-1] &= bip->n->mask; /* do not invert last bits (do not belong to bipartition) */
+ bip->n_ones = bip->n->bits - bip->n_ones;
+ return;
+}
+
+bool
+bipartition_is_bit_set (const bipartition bip, int position)
+{
+ if (bip->bs[(int)(position/BitStringSize)] & (1LL << (int)(position%BitStringSize))) return true;
+ return false;
+}
+
+bool
+bipartition_contains_bits (const bipartition b1, const bipartition b2)
+{ /* generalization of bipartition_is_bit_set(); b1 contains or not b2 */
+ int i;
+ if (b1->n_ones < b2->n_ones) return false;
+ for (i=0; i < b1->n->ints; i++) if ((b2->bs[i]) && (b2->bs[i] != (b1->bs[i] & b2->bs[i]))) return false;
+ return true;
+}
+
+void
+bipartition_to_int_vector (const bipartition b, int *id, int vecsize)
+{
+ int i, j, k = 0;
+ for (i=0; i < b->n->ints; i++) for (j=0; (j < BitStringSize) && (k < vecsize); j++) if ( ((b->bs[i] >> j) & 1LL) ) id[k++] = i * BitStringSize + j;
+}
+
+/*
+void
+bipartition_print_to_stdout (const bipartition b1)
+{
+ int i, j;
+ for (i = 0; i < b1->n->ints - 1; i++) {
+ for (j = 0; j < BitStringSize; j++) printf ("%d", (int)((b1->bs[i] >> j) & 1LL));
+ printf (".");
+ }
+ for (j = 0; j < b1->n->bits%BitStringSize; j++) printf ("%d", (int)((b1->bs[i] >> j) & 1LL));
+ printf ("[%d] ", b1->n_ones);
+}
+*/
+
+void
+bipartition_replace_bit_in_vector (bipartition *bvec, int n_b, int to, int from, bool reduce)
+{ /* copy info from position "from" to position "to" */
+ int k, j = from%BitStringSize, i = from/BitStringSize, j2 = to%BitStringSize, i2 = to/BitStringSize;
+
+ /* boolean "reduce" means that bitstring space will be reduced (last bits will be removed), therefore the update of
+ * n_ones is different from default bipartition_set() behaviour: it's not an extra "1" (that is, one that did not
+ * contribute to n_ones), but an existing "1" that change places. Schematically:
+ * from -> to | normal n_ones count | when bitstring is reduced afterwards
+ * 0 -> 0 | 0 | 0
+ * 0 -> 1 | -1 | -1
+ * 1 -> 0 | +1 | 0 (since it's a leaf that belonged to position "from" and now is on position "to")
+ * 1 -> 1 | 0 | -1 (in fact one of the two "1"s dissapeared after reducing the bitstring
+ * (the above description is outdated since I rewrote by hand the bit functions -- observe how we must erase 1 values from "from") */
+ if (reduce) for (k = 0; k < n_b; k++) { // copy 0 or 1 values, erasing "from" values to avoid problems after reducing space (hanging 1s out of range)
+ if ( ((bvec[k]->bs[i] >> j) & 1LL) && ((bvec[k]->bs[i2] >> j2) & 1LL) ) { bvec[k]->n_ones--; bvec[k]->bs[i] &= ~(1LL << j); }
+ else if ( ((bvec[k]->bs[i] >> j) & 1LL) && !((bvec[k]->bs[i2] >> j2) & 1LL) ) { bvec[k]->bs[i2] |= (1LL << j2); bvec[k]->bs[i] &= ~(1LL << j); }
+ else if ( !((bvec[k]->bs[i] >> j) & 1LL) && ((bvec[k]->bs[i2] >> j2) & 1LL) ) { bvec[k]->bs[i2] &= ~(1LL << j2); bvec[k]->n_ones--; }
+ /* else do nothing (from zero to zero) */
+ }
+
+ else for (k = 0; k < n_b; k++) { // copy 0 or 1 values
+ if ( ((bvec[k]->bs[i] >> j) & 1LL) ) bipartition_set_lowlevel (bvec[k], i2, j2); // will check if n_ones change or not
+ else bipartition_unset_lowlevel (bvec[k], i2, j2);
+ }
+}
+
+void
+bipartition_resize_vector (bipartition *bvec, int n_b)
+{
+ int k, i = bvec[0]->n->ints - 1;
+ for (k = 0; k < n_b; k++) { bvec[k]->bs[i] &= bvec[0]->n->mask; bipartition_count_n_ones (bvec[k]); }
+}
+
+
diff --git a/tests/testthat/test_dist_tree.R b/tests/testthat/test_dist_tree.R
new file mode 100644
index 0000000..c6e3af0
--- /dev/null
+++ b/tests/testthat/test_dist_tree.R
@@ -0,0 +1,22 @@
+context("dist_tree")
+
+data(Laurasiatherian)
+dm <- dist.ml(Laurasiatherian)
+
+
+test_that("check nnls functions", {
+ tree_nj <- NJ(dm)
+ tree_unj <- UNJ(dm)
+ tree_nnls_unj <- nnls.phylo(tree_unj, dm)
+ tree_nnls_nj <- nnls.phylo(tree_nj, dm)
+
+ tree_upgma <- upgma(dm)
+ tree_wpgma <- wpgma(dm)
+ tree_nnls_upgma <- nnls.phylo(tree_upgma, dm, rooted=TRUE)
+ tree_nnls_wpgma <- nnls.phylo(tree_wpgma, dm, rooted=TRUE)
+
+ expect_equal(tree_upgma, tree_nnls_upgma)
+ expect_false(all.equal(tree_wpgma, tree_nnls_wpgma))
+ expect_equal(tree_unj, tree_nnls_unj)
+ expect_false(all.equal(tree_nj, tree_nnls_nj))
+})
diff --git a/tests/testthat/test_distances.R b/tests/testthat/test_distances.R
index f1ac117..63bb2a7 100644
--- a/tests/testthat/test_distances.R
+++ b/tests/testthat/test_distances.R
@@ -8,11 +8,27 @@ weights <- as.vector(1000*exp(fit$site))
attr(X, "weight") <- weights
dm <- cophenetic(tree)
+Y <- phyDat(matrix(c("A", "C", "G", "T", "A", "C", "G", "A"), 2, 4, dimnames=list(c("a", "b", NULL)), byrow=TRUE))
+
+fun <- function(s) - 3/4 * log(1 - 4/3 * s)
+
test_that("dist.ml works properly", {
- skip_on_cran()
+# skip_on_cran()
expect_that(dist.logDet(X), is_a("dist"))
expect_that(dist.hamming(X), is_a("dist"))
expect_that(dist.ml(X), is_a("dist"))
- all.equal(as.matrix(dist.ml(X, k=4, shape=.5)), dm)
-})
\ No newline at end of file
+ expect_equal(as.matrix(dist.ml(X, k=4, shape=.5)), dm)
+ expect_equal(as.matrix(dist.ml(Y)), as.matrix(fun(dist.hamming(Y))))
+})
+
+
+test_that("read/write of distances works", {
+ skip_on_cran()
+ writeDist(dm, "dm.txt")
+ expect_equal(as.dist(dm), readDist("dm.txt"))
+ unlink("dm.txt")
+})
+
+
+
diff --git a/tests/testthat/test_hadamard.R b/tests/testthat/test_hadamard.R
new file mode 100644
index 0000000..717acea
--- /dev/null
+++ b/tests/testthat/test_hadamard.R
@@ -0,0 +1,21 @@
+context("Hadamard conjugation")
+
+v <- 1:8
+data(yeast)
+dm <- dist.hamming(yeast)
+# RY-coding
+yeast_ry <- acgt2ry(yeast)
+# delete ambiguous states
+# dat4 <- phyDat(as.character(yeast), type="USER", levels=c("a","c", "g", "t"), ambiguity=NULL)
+# fit4 <- h4st(dat4)
+
+test_that("Hadamard conjugation works as expected", {
+ # fast fft like multiplication
+ expect_is(H <- hadamard(3), "matrix")
+ expect_equal(as.vector(H %*% v), fhm(v))
+ expect_is(spl_ry <- h2st(yeast_ry), "splits")
+ expect_is(spl_dm <- distanceHadamard(dm), "splits")
+})
+
+
+
diff --git a/tests/testthat/test_parsimony.R b/tests/testthat/test_parsimony.R
index ee2ac27..48542a9 100644
--- a/tests/testthat/test_parsimony.R
+++ b/tests/testthat/test_parsimony.R
@@ -1,5 +1,8 @@
context("parsimony")
+data(yeast)
+all_trees <- allTrees(8, tip.label = names(yeast))
+
tree1 = read.tree(text = "((t1,t2),t3,t4);")
tree2 = read.tree(text = "((t1,t3),t2,t4);")
trees = .compressTipLabel(c(tree1, tree2))
@@ -14,11 +17,40 @@ test_that("parsimony works properly", {
expect_that(fitch(trees, dat), equals(c(1,2)))
expect_that(sankoff(tree1, dat), equals(1))
expect_that(sankoff(tree2, dat), equals(2))
-
expect_that(parsimony(tree1, dat), equals(1))
})
+test_that("bab works properly", {
+ skip_on_cran()
+# all_trees <- allTrees(8, tip.label = names(yeast))
+ all_pars <- fitch(all_trees, yeast)
+ bab_tree <- bab(yeast)
+ expect_equal(min(all_pars), fitch(bab_tree, yeast))
+})
+
+test_that("rearrangements works properly", {
+ skip_on_cran()
+ tree <- all_trees[[1]]
+ start <- fitch(tree, yeast)
+ bab_tree <- bab(yeast)
+ best <- fitch(bab_tree, yeast)
+ best_fitch <- optim.parsimony(tree, yeast, rearrangements = "NNI")
+ best_sankoff <- optim.parsimony(tree, yeast, method="sankoff", rearrangements = "NNI")
+ expect_equal(attr(best_fitch, "pscore"), attr(best_sankoff, "pscore"))
+})
+
+test_that("tree length works properly", {
+ skip_on_cran()
+ tree <- nj(dist.hamming(yeast))
+ pscore <- fitch(tree, yeast)
+ tree1 <- acctran(tree, yeast)
+ expect_equal(sum(tree1$edge.length), pscore)
+ tree2 <- rtree(100)
+ dat <- simSeq(tree2)
+ tree2 <- acctran(tree2, dat)
+ expect_equal(sum(tree2$edge.length), fitch(tree2,dat))
+})
diff --git a/tests/testthat/test_phyDat.R b/tests/testthat/test_phyDat.R
index 820165a..b6b2fbb 100644
--- a/tests/testthat/test_phyDat.R
+++ b/tests/testthat/test_phyDat.R
@@ -1,19 +1,77 @@
context("conversion_and_subsetting")
-
+# to test: phyDat.codon,
data(Laurasiatherian)
+data(chloroplast)
+set.seed(42)
+tree <- rtree(10)
+codon_align <- simSeq(tree, l=100, type = "CODON")
+
+
phy_matrix <- as.character(Laurasiatherian)
phy_df <- as.data.frame(Laurasiatherian)
+phy_vec_dna <- phy_matrix[,1]
+phy_vec_user <- sample(c("0","1"), 26, replace=TRUE)
+names(phy_vec_user) <- letters
phy_dnabin <- as.DNAbin(Laurasiatherian)
phy_align <- phyDat2alignment(Laurasiatherian)
test_that("conversion work as expected", {
+## skip_on_cran()
+ expect_is(phy_matrix, "matrix")
+ expect_is(phy_df, "data.frame")
+ expect_is(phy_dnabin, "DNAbin")
+ expect_is(phy_align, "alignment")
+ expect_is(as.phyDat(phy_matrix), "phyDat")
+ expect_is(as.phyDat(phy_df), "phyDat")
+ expect_is(as.phyDat(phy_dnabin), "phyDat")
+ expect_is(phyDat(phy_vec_dna), "phyDat")
+ expect_is(phyDat(phy_vec_user, type="USER", levels = c("0","1")), "phyDat")
+ expect_is(as.phyDat(phy_dnabin), "phyDat")
+ expect_is(as.phyDat(phy_align), "phyDat")
+ expect_is(c2d <- codon2dna(codon_align), "phyDat")
+ expect_equal(dna2codon(c2d), codon_align)
+})
+
+
+#test_that("conversion with Biostrings work as expected", {
+# skip_on_cran()
+# library(Biostrings)
+# expect_is(MA_AA <- as.MultipleAlignment(chloroplast), "AAMultipleAlignment")
+# expect_equal(as.phyDat(MA_AA), chloroplast)
+# expect_is(MA_DNA <- as.MultipleAlignment(Laurasiatherian), "DNAMultipleAlignment")
+# expect_equal(as.phyDat(MA_DNA), Laurasiatherian)
+#})
+
+
+test_that("subsetting and combining work as expected", {
+## skip_on_cran()
+
+ expect_is(subset_1 <- subset(Laurasiatherian, select = 1:1000, site.pattern = FALSE), "phyDat")
+ expect_is(subset_2 <- subset(Laurasiatherian, select = 1001:3179, site.pattern = FALSE), "phyDat")
+ expect_is(lauraCbind1 <- cbind(subset_1, subset_2), "phyDat")
+ expect_equal(baseFreq(lauraCbind1), baseFreq(Laurasiatherian))
+
+ expect_is(subset_3 <- subset(Laurasiatherian, select = 1:100), "phyDat")
+ expect_is(subset_4 <- subset(Laurasiatherian, select = 101:1605), "phyDat")
+ expect_is(lauraCbind2 <- cbind(subset_3, subset_4), "phyDat")
+ expect_equal(baseFreq(lauraCbind2), baseFreq(Laurasiatherian))
+})
+
+
+test_that("read and write work as expected", {
skip_on_cran()
- expect_that(phy_matrix, is_a("matrix"))
- expect_that(phy_df, is_a("data.frame"))
- expect_that(phy_dnabin, is_a("DNAbin"))
- expect_that(phy_align, is_a("alignment"))
- expect_that(as.phyDat(phy_matrix), is_a("phyDat"))
- expect_that(as.phyDat(phy_df), is_a("phyDat"))
- expect_that(as.phyDat(phy_dnabin), is_a("phyDat"))
- expect_that(as.phyDat(phy_align), is_a("phyDat"))
-})
\ No newline at end of file
+
+ write.phyDat(Laurasiatherian, "tmp1.txt")
+ expect_is(laura <- read.phyDat("tmp1.txt"), "phyDat")
+ expect_equal(laura, Laurasiatherian)
+ unlink("tmp1.txt")
+ write.phyDat(chloroplast, "tmp2.txt")
+ expect_is(chloro <- read.phyDat("tmp2.txt", type="AA"), "phyDat")
+ expect_equal(chloro, chloroplast)
+ unlink("tmp2.txt")
+ write.phyDat(chloroplast, "tmp.fas", format="fasta")
+ expect_is(chloro_fas <- read.phyDat("tmp.fas", type="AA", format = "fasta"), "phyDat")
+ expect_equal(chloro_fas, chloroplast)
+ unlink("tmp.fas")
+})
+
diff --git a/tests/testthat/test_pml.R b/tests/testthat/test_pml.R
index 949c4c0..08cb387 100644
--- a/tests/testthat/test_pml.R
+++ b/tests/testthat/test_pml.R
@@ -16,7 +16,9 @@ dat = allSitePattern(5)
weights = as.vector(1000 * exp(pml(treeR1, dat)$siteLik))
attr(dat, "weight") = weights
-bf = c(.1,.2,.3,.4)
+dat_Mk <- subset(dat, select = -c(1,342,683, 1024))
+
+
Q = c(6:1)
@@ -33,17 +35,6 @@ pmlR3 = pml(treeR3, dat)
pmlR3.fitted = optim.pml(pmlR3, TRUE, optRooted = TRUE, control = pml.control(epsilon=1e-10, trace=0))
-
-# optim.pml:
-# bf F81
-# Q Sym
-# Gamma
-# Inv
-
-
-#tr_acctran = acctran(tree1, dat)
-#tr_ratchet = pratchet(dat, trace=0)
-#bab(dat)
test_that("edge length optimisation works properly", {
skip_on_cran()
expect_equal(logLik(pmlU2.fitted), logLik(pmlU1))
@@ -59,3 +50,81 @@ test_that("NNI optimisation works properly", {
expect_equal(pmlU3.fitted$tree, pmlU1$tree, tolerance=1e-6)
# expect_equal(pmlR3.fitted$tree, pmlR1$tree, tolerance=5e-6)
})
+
+
+
+test_that("bf optimisation works properly", {
+ skip_on_cran()
+ bf <- c(.1,.2,.3,.4)
+ fit_T <- pml(treeU1, dat, bf=bf)
+ weights <- as.vector(1000 * exp(fit_T$siteLik))
+ dat_tmp <- dat
+ attr(dat_tmp, "weight") <- weights
+ fit0 <- pml(treeU1, dat_tmp)
+ fit.bf <- optim.pml(fit0, optEdge=FALSE, optBf = TRUE, control = pml.control(epsilon=1e-10, trace=0))
+ expect_equal(logLik(fit.bf), logLik(pml(treeU1, dat_tmp, bf=bf)))
+ expect_equal(bf, fit.bf$bf, tolerance=5e-4)
+})
+
+
+test_that("Q optimisation works properly", {
+ skip_on_cran()
+ Q <- c(6:1)
+ fit_T <- pml(treeU1, dat, Q=Q)
+ weights <- as.vector(1000 * exp(fit_T$siteLik))
+ dat_tmp <- dat
+ attr(dat_tmp, "weight") <- weights
+ fit0 <- pml(treeU1, dat_tmp)
+ fit.Q <- optim.pml(fit0, optEdge=FALSE, optQ = TRUE, control = pml.control(epsilon=1e-10, trace=0))
+ expect_equal(logLik(fit.Q), logLik(pml(treeU1, dat_tmp, Q=Q)))
+ expect_equal(Q, fit.Q$Q, tolerance=5e-4)
+})
+
+
+test_that("Inv optimisation works properly", {
+ skip_on_cran()
+ inv <- 0.25
+ fit_T <- pml(treeU1, dat, inv=inv)
+ weights <- as.vector(1000 * exp(fit_T$siteLik))
+ dat_tmp <- dat
+ attr(dat_tmp, "weight") <- weights
+ fit0 <- pml(treeU1, dat_tmp)
+ fit.Inv <- optim.pml(fit0, optEdge=FALSE, optInv = TRUE, control = pml.control(epsilon=1e-10, trace=0))
+ expect_equal(logLik(fit.Inv), logLik(pml(treeU1, dat_tmp, inv=inv)))
+ expect_equal(inv, fit.Inv$inv, tolerance=5e-4)
+})
+
+
+test_that("Gamma optimisation works properly", {
+ skip_on_cran()
+ shape <- 2
+ fit_T <- pml(treeU1, dat, shape=shape, k=4)
+ weights <- as.vector(1000 * exp(fit_T$siteLik))
+ dat_tmp <- dat
+ attr(dat_tmp, "weight") <- weights
+ fit0 <- pml(treeU1, dat_tmp, k=4)
+ fit.Gamma <- optim.pml(fit0, optEdge=FALSE, optGamma = TRUE, control = pml.control(epsilon=1e-10, trace=0))
+ expect_equal(logLik(fit.Gamma), logLik(pml(treeU1, dat_tmp, shape = shape, k=4)))
+ expect_equal(shape, fit.Gamma$shape, tolerance=5e-4)
+})
+
+
+test_that("rate optimisation works properly", {
+ skip_on_cran()
+ rate <- 2
+ fit_T <- pml(treeU1, dat, rate=rate)
+ weights <- as.vector(1000 * exp(fit_T$siteLik))
+ dat_tmp <- dat
+ attr(dat_tmp, "weight") <- weights
+ fit0 <- pml(treeU1, dat_tmp)
+ fit.rate <- optim.pml(fit0, optEdge=FALSE, optRate = TRUE, control = pml.control(epsilon=1e-10, trace=0))
+ expect_equal(logLik(fit.rate), logLik(pml(treeU1, dat_tmp, rate=rate)))
+ expect_equal(rate, fit.rate$rate, tolerance=5e-4)
+})
+
+test_that("Mkv model works properly", {
+ skip_on_cran()
+ expect_equal(logLik(pmlU2.fitted), logLik(pmlU1))
+})
+
+## 31.91%
\ No newline at end of file
diff --git a/tests/testthat/test_splits.R b/tests/testthat/test_splits.R
index 1f2e20e..569e5ce 100644
--- a/tests/testthat/test_splits.R
+++ b/tests/testthat/test_splits.R
@@ -8,16 +8,55 @@ spl2tree <- as.phylo(tree2spl)
dm <- cophenetic(tree2spl)
mat <- as.matrix(tree2spl)
Mat <- as.Matrix(tree2spl)
+trees <- nni(tree)
test_that("splits ", {
## skip on CRAN
skip_on_cran()
## check classes
+ expect_is(as.splits(trees), "splits")
expect_is(tree2spl, "splits")
expect_is(spl2tree,"phylo")
expect_is(dm,"dist")
expect_is(mat, "matrix")
expect_is(Mat, "Matrix")
expect_equal(spl2tree , tree)
+
+ spl <- allCircularSplits(6)
+ spl <- phangorn:::oneWise(spl, 6)
+ write.nexus.splits(spl, "tmp.nex")
+ spl2 <- read.nexus.splits("tmp.nex")
+ attr(spl2, "splitlabels") <- NULL
+ attr(spl2, "weights") <- NULL
+ class(spl2) <- "splits"
+ expect_equal(spl2 , spl)
+ unlink("tmp.nex")
+})
+
+
+test_that("networx ", {
+ net1 <- neighborNet(dm)
+ write.nexus.networx(net1, "tmp.nex")
+ net2 <- read.nexus.networx("tmp.nex")
+ net3 <- as.networx(tree)
+ # delete some additional attributes
+ net2$.plot <- net2$translate <- NULL
+ attr(net1, "order") = NULL
+
+ expect_is(net1, "networx")
+ expect_is(net2, "networx")
+ expect_is(net3, "networx")
+# expect_equal(net1, net2, tolerance=1e-6)
+# expect_equal(net3, net2, tolerance=1e-6)
+ expect_equal(net1, net3)
+ unlink("tmp.nex")
+ cnet <- consensusNet(as.splits(trees))
+ expect_is(cnet, "networx")
+
+ net1$edge.length <- cnet$edge.length <- cnet$edge.labels <- NULL
+ attr(cnet, "order") = NULL
+ expect_equal(cnet, net1)
})
+
+
diff --git a/tests/testthat/test_superTree.R b/tests/testthat/test_superTree.R
new file mode 100644
index 0000000..40bd660
--- /dev/null
+++ b/tests/testthat/test_superTree.R
@@ -0,0 +1,26 @@
+context("superTree")
+
+tree <- rtree(50, rooted=FALSE)
+
+trees_simple <- nni(tree)
+
+trees <- rNNI(tree, sample(10, 100, replace = TRUE))
+trees <- .uncompressTipLabel(trees)
+labels <- paste0("t", 1:50)
+trees2 <- lapply(trees, function(x)drop.tip(x, sample(labels, 10)))
+class(trees2) <- "multiPhylo"
+
+test_that("superTree", {
+ ## check superTree
+ skip_on_cran()
+ simple_superTree <- superTree(trees_simple, rooted=FALSE)
+ difficult_superTree <- superTree(trees2, rooted=FALSE)
+ expect_equal(RF.dist(simple_superTree, tree) , 0)
+ expect_equal(RF.dist(difficult_superTree, tree) , 0)
+ rf_superTree <- superTree(trees, method="RF")
+ spr_superTree <- superTree(trees, method="SPR")
+ expect_gte(attr(rf_superTree, "score"), attr(spr_superTree, "score"))
+})
+
+
+
diff --git a/tests/testthat/test_treeManipulation.R b/tests/testthat/test_treeManipulation.R
index 006d785..f3aadff 100644
--- a/tests/testthat/test_treeManipulation.R
+++ b/tests/testthat/test_treeManipulation.R
@@ -1,7 +1,11 @@
context("treeManipulation")
set.seed(42)
-tree <- rtree(100)
+tree <- rtree(100, rooted=FALSE)
+
+
+
+trees <- lapply(sample(10:500,50), function(x)tree <- rtree(x, rooted=FALSE) )
#Ancestors(x, node, type=c("all","parent"))
#Children(x, node)
@@ -17,19 +21,40 @@ node_108 <- mrca.phylo(tree, node=desc_108)
test_that("ancestor, mrca, descendants", {
## check RF.dist and tree dist
- expect_equal(mrca.phylo(tree, node=desc_108), 108)
+ expect_equal(mrca.phylo(tree, node=desc_108), 108L)
+ kids_108 <- Descendants(tree, 108, "children")
+ expect_equal(Ancestors(tree, kids_108, "parent"), rep(108L, length(kids_108)))
+ expect_equal(Siblings(tree, kids_108[1], include.self=TRUE), kids_108)
+
})
-test_that("allTrees, nni", {
+test_that("allTrees", {
## allTrees
expect_is(allTrees(6), "multiPhylo")
expect_true(all(RF.dist(allTrees(6))>0))
+})
+
+
+test_that("nni", {
## nni
expect_is(nni(tree), "multiPhylo")
expect_true(all(RF.dist(nni(tree), tree)>0))
})
+# TODO: check why rooted trees give error in development version
+test_that("midpoint", {
+ # topology stays the same
+ expect_equal( max( sapply(trees, function(x)RF.dist(x,midpoint(x)))), 0)
+ # 2 * max(height) == max(cophenetic)
+ expect_equal( max( node.depth.edgelength(midpoint(tree)) *2) , max(cophenetic(tree)))
+})
+
+test_that("maxCladeCred", {
+ tree <- rcoal(100)
+ trees <- nni(tree)
+ expect_equal(maxCladeCred(c(tree, trees)), tree)
+})
diff --git a/tests/testthat/test_treedist.R b/tests/testthat/test_treedist.R
index 2f776ef..9f564d4 100644
--- a/tests/testthat/test_treedist.R
+++ b/tests/testthat/test_treedist.R
@@ -33,11 +33,14 @@ test_that("Kuhner-Felsenstein distance (branch score difference)", {
expect_that(treedist(tree1, tree1)[2], is_equivalent_to(0))
expect_that(treedist(tree1, tree2)[2], is_equivalent_to(sqrt(2)))
expect_that(treedist(tree1, tree3)[2], is_equivalent_to(1))
+
+ expect_equal(KF.dist(tree1, c(tree1, tree2, tree3)), c(0, sqrt(2), 1))
})
test_that("path distance", {
## check path dist
expect_that(path.dist(tree1, tree1), is_equivalent_to(0))
+ expect_equal(path.dist(trees)[1:13] , path.dist(trees[[1]], trees[2:14]))
expect_is(path.dist(trees),"dist")
})
@@ -73,5 +76,29 @@ test_that("When each tree has unit branch lengths, RF = wRF", {
max(abs(RF.dist(trees) - wRF.dist(trees))) # find maximum abs difference between RF and wRF distance (expect 0)
})), # find max of all these
0) # expect equal to 0
+ # now the same for comparison of one tree with many
+ expect_equal(
+ max(
+ sapply(sample(10:500,50), function(x) { # test some random numbers of tips between 10 and 500
+ trees <- rmtree(20, x, rooted=FALSE, br=1) # generate 20 unrooted trees with the given number of tips and unit branch lengths
+ max(abs(RF.dist(trees[[1]], trees) - wRF.dist(trees[[1]], trees))) # find maximum abs difference between RF and wRF distance (expect 0)
+ })), # find max of all these
+ 0) # expect equal to 0
})
+#############################
+# test sprdist from leomrtns
+#############################
+test_that("SPR distance", {
+ ## check spr dist
+ set.seed(123)
+ tree1 <- rtree(100, rooted = FALSE)
+ tree2 <- rSPR(tree1, 1)
+ trees <- rSPR(tree1, 1:5)
+ expect_equal(sprdist(tree1, tree2)[[1]], 1)
+ expect_equal(sprdist(tree1, tree2)[[3]], RF.dist(tree1, tree2))
+ expect_equal(SPR.dist(tree1, trees), 1:5)
+ expect_is(SPR.dist(trees), "dist")
+})
+
+
diff --git a/vignettes/Networx.Rmd b/vignettes/Networx.Rmd
index 129ba78..05c69ab 100644
--- a/vignettes/Networx.Rmd
+++ b/vignettes/Networx.Rmd
@@ -11,7 +11,7 @@ vignette: >
---
-This tutorial gives a basic introduction for constructing phylogenetic networks and adding parameters to trees or networx objects using [phangorn](http://cran.r-project.org/package=phangorn) [@Schliep2011] in R.
+This tutorial gives a basic introduction for constructing phylogenetic networks and adding parameters to trees or networx objects using [phangorn](https://cran.r-project.org/package=phangorn) [@Schliep2011] in R.
Splits graphs or phylogenetic networks are a useful way to display conflicting data or to summarize different trees. Here, we present two popular networks, consensus networks [@Holland2004]
and Neighbor-Net [@Bryant2004].
Trees or networks are often missing either edge weights or edge support values. We show how to improve a tree/networx object by adding support values or estimating edge weights using non-negative Least-Squares (nnls).
diff --git a/vignettes/phangorn.bib b/vignettes/phangorn.bib
index 2e6c227..5523152 100644
--- a/vignettes/phangorn.bib
+++ b/vignettes/phangorn.bib
@@ -12,6 +12,35 @@
}
+ at article{Leo2008,
+ author = {de Oliveira Martins, Leonardo AND Leal, Élcio AND Kishino, Hirohisa},
+ journal = {PLoS ONE},
+ publisher = {Public Library of Science},
+ title = {Phylogenetic Detection of Recombination with a Bayesian Prior on the Distance between Trees},
+ year = {2008},
+ month = {07},
+ volume = {3},
+ url = {http://dx.plos.org/10.1371%2Fjournal.pone.0002651},
+ pages = {1-13},
+ number = {7},
+ doi = {10.1371/journal.pone.0002651}
+}
+
+
+ at article{Leo2016,
+ author = {De Oliveira Martins, Leonardo and Mallo, Diego and Posada, David},
+ title = {A Bayesian Supertree Model for Genome-Wide Species Tree Reconstruction},
+ volume = {65},
+ number = {3},
+ pages = {397-416},
+ year = {2016},
+ oi = {10.1093/sysbio/syu082},
+ URL = {http://sysbio.oxfordjournals.org/content/65/3/397.abstract},
+ eprint = {http://sysbio.oxfordjournals.org/content/65/3/397.full.pdf+html},
+ journal = {Systematic Biology}
+}
+
+
@article{Schwarz1978,
author = "Schwarz, Gideon",
doi = "10.1214/aos/1176344136",
@@ -234,6 +263,32 @@ algorithms},
}
+ at article {Paradis2017,
+ author = {Paradis, Emmanuel and Gosselin, Thierry and Goudet, Jérôme and Jombart, Thibaut and Schliep, Klaus},
+ title = {Linking genomics and population genetics with R},
+ journal = {Molecular Ecology Resources},
+ issn = {1755-0998},
+ url = {http://dx.doi.org/10.1111/1755-0998.12577},
+ doi = {10.1111/1755-0998.12577},
+ pages = {n/a--n/a},
+ keywords = {multivariate analysis, NGS, R, SNP, VCF},
+ year = {2017},
+}
+
+
+ at article {Jombart2017,
+ author = {Jombart, Thibaut and Archer, Frederick and Schliep, Klaus and Kamvar, Zhian and Harris, Rebecca and Paradis, Emmanuel and Goudet, Jérome and Lapp, Hilmar},
+ title = {apex: phylogenetics with multiple genes},
+ journal = {Molecular Ecology Resources},
+ issn = {1755-0998},
+ url = {http://dx.doi.org/10.1111/1755-0998.12567},
+ doi = {10.1111/1755-0998.12567},
+ pages = {n/a--n/a},
+ keywords = {R, package, Software, Genetics, Phylogenies},
+ year = {2017},
+}
+
+
@article{Kuhner1994,
author = {Kuhner, M K and Felsenstein, J},
title = {A simulation comparison of phylogeny algorithms under equal and unequal evolutionary rates.},
@@ -515,6 +570,30 @@ journal = {Molecular Biology and Evolution}
Year = "2005"
}
+ at article{Pupko2000,
+author = {Pupko, Tal and Pe, Itsik and Shamir, Ron and Graur, Dan},
+title = {A Fast Algorithm for Joint Reconstruction of Ancestral Amino Acid Sequences},
+volume = {17},
+number = {6},
+pages = {890-896},
+year = {2000},
+eprint = {http://mbe.oxfordjournals.org/content/17/6/890.full.pdf+html},
+journal = {Molecular Biology and Evolution}
+}
+
+ at article{Pupko2002,
+author = {Pupko, Tal and Pe'er, Itsik and Hasegawa, Masami and Graur, Dan and Friedman, Nir},
+title = {A branch-and-bound algorithm for the inference of ancestral amino-acid sequences when the replacement rate varies among sites: Application to the evolution of five gene families},
+volume = {18},
+number = {8},
+pages = {1116-1123},
+year = {2002},
+doi = {10.1093/bioinformatics/18.8.1116},
+URL = {http://bioinformatics.oxfordjournals.org/content/18/8/1116.abstract},
+eprint = {http://bioinformatics.oxfordjournals.org/content/18/8/1116.full.pdf+html},
+journal = {Bioinformatics}
+}
+
@InCollection{ Swofford1996,
Author = "Swofford, D.L. and Olsen, G.J. and Waddell, P.J. and Hillis, D.M.",
Title = "Phylogenetic Inference",
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/r-cran-phangorn.git
More information about the debian-med-commit
mailing list