[med-svn] [r-cran-phangorn] 01/02: Imported Upstream version 1.99.14
Andreas Tille
tille at debian.org
Sat Aug 22 14:23:43 UTC 2015
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 27047e457b8859de5e21a38ffaf87d32616f11f4
Author: Andreas Tille <tille at debian.org>
Date: Sat Aug 22 16:20:10 2015 +0200
Imported Upstream version 1.99.14
---
DESCRIPTION | 16 +-
MD5 | 79 +--
NAMESPACE | 28 +-
NEWS | 46 ++
R/clanistic.R | 2 +-
R/distSeq.R | 19 +
R/distTree.R | 554 +++++++++++++++++++++
R/fitch.R | 2 +-
R/modelTest.R | 201 ++++++++
R/neighborNet.R | 9 +-
R/networx.R | 20 +-
R/parsimony.R | 10 +-
R/phyDat.R | 58 ++-
R/phylo.R | 1010 +++++++++++----------------------------
R/treeManipulation.R | 1 +
R/treedist.R | 39 +-
README.md | 6 +
TODO | 15 +
build/vignette.rds | Bin 335 -> 335 bytes
inst/doc/Networx.R | 14 +-
inst/doc/Networx.Rmd | 2 +-
inst/doc/Networx.html | 2 +-
inst/doc/Trees.R | 2 +-
inst/doc/Trees.Rnw | 2 +-
man/bootstrap.pml.Rd | 5 +-
man/densiTree.Rd | 3 +
man/dist.hamming.Rd | 6 +
man/modelTest.Rd | 7 +-
man/parsimony.Rd | 4 +-
man/phyDat.Rd | 17 +-
man/plot.networx.Rd | 2 +-
man/pml.fit.Rd | 3 +
man/splitsNetwork.Rd | 2 +-
man/treedist.Rd | 3 +-
src/fitch.c | 42 +-
src/ml.c | 495 +++++++++++--------
src/sankoff.c | 2 +-
tests/testthat.R | 4 +
tests/testthat/test_parsimony.R | 21 +
tests/testthat/test_phyDat.R | 19 +
tests/testthat/test_treedist.R | 14 +
vignettes/Networx.Rmd | 2 +-
vignettes/Trees.Rnw | 2 +-
vignettes/phangorn.bib | 203 +++++++-
44 files changed, 1882 insertions(+), 1111 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index d2cf18d..a10e037 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,19 +1,21 @@
Package: phangorn
Title: Phylogenetic Analysis in R
-Version: 1.99-13
-Date: 2015-04-07
+Version: 1.99.14
+Date: 2015-07-09
Author: Klaus Schliep, Emmanuel Paradis
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.2)
-Imports: quadprog, igraph (>= 0.6), Matrix, parallel, nnls
-Suggests: seqLogo, seqinr, xtable, flashClust, rgl, knitr
+Depends: R (>= 3.0.0), ape (>= 3.3)
+Imports: quadprog, igraph (>= 0.6), Matrix, parallel, nnls, methods,
+ utils, stats, graphics, grDevices
+Suggests: testthat, seqLogo, seqinr, xtable, flashClust, rgl, knitr
+ByteCompile: TRUE
License: GPL (>= 2)
VignetteBuilder: utils, knitr
URL: https://github.com/KlausVigo/phangorn
Repository: CRAN
-Packaged: 2015-04-07 21:11:16 UTC; klaus
NeedsCompilation: yes
-Date/Publication: 2015-04-08 00:44:28
+Packaged: 2015-07-09 14:57:03 UTC; klaus
+Date/Publication: 2015-07-09 18:11:29
diff --git a/MD5 b/MD5
index ab6d131..a2a4a4f 100644
--- a/MD5
+++ b/MD5
@@ -1,28 +1,31 @@
-ba46ad1e822741dcea4f2fd1269a5742 *DESCRIPTION
-c699d9d5eb46e6f75a81ce923cc3411b *NAMESPACE
-0c9b7550aefe626cca7d81bc5f24d59b *NEWS
+6935e0c843976a470143b72ada3f1785 *DESCRIPTION
+2d67e088ea29f95611ee42b755b256e0 *NAMESPACE
+d091257060fa258698275d7f93abfb41 *NEWS
115bc297835dcaab7d32527fcacb2061 *R/Coalescent.R
874ee7c405b1384b2f04828ac27574ef *R/Densi.R
f916be823b9a838d49145695b0a37aa1 *R/SOWH.R
188e8e2a81eb0c50cc6b7cda6d2c9b98 *R/cladePar.R
-db718843face57068ea5d842cccc22e2 *R/clanistic.R
+5a95cd8530b05882c63f5d48073d7d3f *R/clanistic.R
e4992df0ce2eaeb1380543eed22aa769 *R/dist.p.R
-eb92decd10b1af4636c9dfa3cc97b4b9 *R/distSeq.R
-4dc4a844ea48eb8a694cf8449b23b0af *R/fitch.R
+31db0490edd334101645270e922c3341 *R/distSeq.R
+0693ec904c923ea37447ec460d2a86f7 *R/distTree.R
+ecbb9305836be819d53745855d9ab3ff *R/fitch.R
86e1095878e3b46eca343d1804a6957e *R/hadamard.R
-0feeebf15c688ec1dd1a003e4c9e488c *R/neighborNet.R
-43346cdbdfd0a819f69a6fbcf8f10ff5 *R/networx.R
-0079fca0b3af71d1730ceab3aeae4f94 *R/parsimony.R
-1378c04926c78c226a3e34d5301f7fa3 *R/phyDat.R
-51b2f1bbf59dd590fba1bd33a5def284 *R/phylo.R
+75c515f0fcb9127a033df76e071eb2ca *R/modelTest.R
+c4c58c9e35b58e21176f482fd6150002 *R/neighborNet.R
+4fbe2e2c5b83e638962921d0b75d4203 *R/networx.R
+6c4a570cf453429cfe60b3f0d23ee69b *R/parsimony.R
+7d2d79ba2ac0f697caa372f487a0f06a *R/phyDat.R
+8579b368fcf48fa7fa0ffed2daf1576b *R/phylo.R
ef07cb79ed0ece0caa494ed7c566e17b *R/sankoff.R
bd5ce765336904ef9fd6ee286dc42220 *R/simSeq.R
2b4d6884b4ee2b92e617cea76f79d3da *R/sysdata.rda
-b10db1358eff8e0d1ed28aff3f318ce5 *R/treeManipulation.R
-d6c76453e0fba7e7691111839479c7f0 *R/treedist.R
+e59cce6bb3c64d7853ce58450f92f7fb *R/treeManipulation.R
+c7639eb4e9f2d2749b4e4aa54555bc75 *R/treedist.R
c925e4cd5ce95bf91968fc909d6eef76 *R/zzz.R
-fd3b3539a74c05bdcc9f368f3ee222ad *README.md
-aa1792f357c0e14cc50a2a2b4bb7facd *build/vignette.rds
+76c6c00aad6e70113fca211f6ee04c6d *README.md
+df6481d8fdef70e63bfc09054cbb7bab *TODO
+c0bb685510ddd82f07347d07b3468cfc *build/vignette.rds
4a92b07bcd170e85bd237e460acafa93 *data/Laurasiatherian.RData
4b269c2e640293341b9b1c04d4dd7f4e *data/chloroplast.RData
19f84425e8caaf9f605490cdff17dd81 *data/yeast.RData
@@ -31,11 +34,11 @@ b4d3fa3f40aae3a68846d1f28274e9a0 *inst/CITATION
cd23e5801a1f3632bf0cdd76e5242279 *inst/doc/Ancestral.R
0b18e27e3877621f05e6f78c89a0f5d2 *inst/doc/Ancestral.Rnw
bfd8e0dfe0c57e211e4600e8221b530d *inst/doc/Ancestral.pdf
-5ae7b5f00020d5544ffe34f27420749b *inst/doc/Networx.R
-efd64c8a6ee13e6443a20adf473d17e4 *inst/doc/Networx.Rmd
-555b48fe21566da8383c1ee92e7a9b81 *inst/doc/Networx.html
-e7c9012ad6d10f1c247564b2562b3eda *inst/doc/Trees.R
-a6ff44c68626895b856effc7e646446d *inst/doc/Trees.Rnw
+8d18b0e4a292379ba270a7a26f4d041e *inst/doc/Networx.R
+6cf134802b4b98bb286e71f135ead8bd *inst/doc/Networx.Rmd
+80f428b2d03d8874db4e9071225057b9 *inst/doc/Networx.html
+6e9e97e846668cbf1e077067764e6283 *inst/doc/Trees.R
+e3f16cf1e79426ba67016beed7ee6011 *inst/doc/Trees.Rnw
46d2167ad02c1ee4670cd07e6420876c *inst/doc/Trees.pdf
1d5112656a0fe77fd68081a8a81474b9 *inst/doc/phangorn-specials.R
6e60535c60981d4f5492b73ebf7344df *inst/doc/phangorn-specials.Rnw
@@ -66,52 +69,56 @@ b3418f878911164dd974046f0492c79c *man/SOWH.test.Rd
2c88d4ded0156a68c44cee9e8fc290b4 *man/ancestral.pml.Rd
52ca504c035177a8ac1dbe43fc8e90d5 *man/as.splits.Rd
3b53888e5f37fd54344286a9175d8dc2 *man/bab.Rd
-22ea01fa824c0992bd55a8a43de76b21 *man/bootstrap.pml.Rd
+1a6c4849b1cd2c672e93ad940d1e8652 *man/bootstrap.pml.Rd
22ad0ef60c2edd8b3d003f170f2fa15a *man/chloroplast.Rd
4361615c7d1dc7beb8a4e6118505423c *man/cladePar.Rd
8d3c6fd646d8a64e29dc79791692a388 *man/consensusNet.Rd
-e99c28d255d25b20b6c02b772e88bacc *man/densiTree.Rd
+5868068eb5f8be4302c0c3caadd25943 *man/densiTree.Rd
f44e43abbbf7707821627c8e88f6e9fa *man/designTree.Rd
924c4e906042d64ac4fc4907af8c6cf1 *man/dfactorial.Rd
-559e1e95a100d3393e508513ecce311d *man/dist.hamming.Rd
+da8a199fb9669c127cedfc3e58a650ad *man/dist.hamming.Rd
e7e01cc12be2fac182ff773826ad286f *man/dist.p.Rd
4f8fa46c2cf1f904d1b85f35fc4db884 *man/distanceHadamard.Rd
72abceeb247f22b3da6560df0c20468a *man/getClans.Rd
c1cacff95f20b574c03126d3f8e24106 *man/hadamard.Rd
9c306689d9d7cd539d8d1055b35e24a2 *man/lento.Rd
f24a38210a0e90d3241e27514a6f5b8a *man/midpoint.Rd
-970c755f7ac1a0d81ca7d38bf2119a9f *man/modelTest.Rd
+0d51482b4a6375d72d220b80a667f917 *man/modelTest.Rd
5a825bb21eeef78dccf760a5ee0aca26 *man/neighborNet.Rd
94cd3c21e4dd4a6919db7fde5146a711 *man/nni.Rd
-f5175ae5af6f9dea5bd0ee51e302e6a0 *man/parsimony.Rd
+f1ffe19a44e661d78c1bff049213d230 *man/parsimony.Rd
06d3770fff0d07f19393875735cfa004 *man/phangorn-package.Rd
-faa3f02d77384b128f4682819b55e04f *man/phyDat.Rd
-1b46f1588ec4e8706df8de61eaa63058 *man/plot.networx.Rd
+71f4c68026c05fc0467f8d2366e629de *man/phyDat.Rd
+d7eac237d2e8950f6aca364bdc349542 *man/plot.networx.Rd
cb484368c04bf20d286dedf974c83da1 *man/pml.Rd
-deadaacd3d541c707db8aab53ca1e35d *man/pml.fit.Rd
+f8182027284ceba292fc437e87af3fad *man/pml.fit.Rd
9506d15a81467d18f40bb4741b3d9d28 *man/pmlCluster.Rd
8b06bdee57c510a690ea04b40bac4844 *man/pmlMix.Rd
ff67d17fa29bf88cea645905334c5ecc *man/pmlPart.Rd
bfc83a71d7de4e1c0b994b8b9c853439 *man/read.aa.Rd
eaadcb69db59ca4c512eeb6258013714 *man/simSeq.Rd
-9bfb43eb7ed5a2a3a36496bfffa965dd *man/splitsNetwork.Rd
+2e7848b9bac018bde56302f70066a49a *man/splitsNetwork.Rd
2614e984b5c5aa40d9364d3542172654 *man/superTree.Rd
-9cac90a57bbdef5d40f016648d431df9 *man/treedist.Rd
+9f0dd1accd30d7b177038b787d3a3296 *man/treedist.Rd
64399c9fb7fb25d103c25aa925ab7a10 *man/upgma.Rd
b97649fe8b91a0632a1ded89a6f43125 *man/yeast.Rd
9a8672b7759d360a54c251d1866b3e05 *src/Makevars
56d140fd1f0fceb9031c4ad70dfb96f1 *src/dist.c
-a750f640db7ff1a3c56a98ccfe02041a *src/fitch.c
-70437536c7263fe33acf2beeaa96e461 *src/ml.c
+c658fdd18a82e58979b0105e0747e8b3 *src/fitch.c
+63c45a73f83920a12ab7f921555a5e11 *src/ml.c
a573d9543f40a6a32c5af7c14ad8993a *src/phangorn.c
d8e5c864aa0b122ce426a67479d3a610 *src/read_aa.c
-2256cf09434a7ff05e750f4a11124aba *src/sankoff.c
+fb760099aedbfcf9df0e638780df4cf3 *src/sankoff.c
+7d1eaef831fb2372b367f5be1b0a5790 *tests/testthat.R
+ee456463a48eac9fb661c9b0b00f8aba *tests/testthat/test_parsimony.R
+3b6b3635995a131b5beb515ae2897e28 *tests/testthat/test_phyDat.R
+a2f72d835a2eb1ad694dc19f1f48242f *tests/testthat/test_treedist.R
0b18e27e3877621f05e6f78c89a0f5d2 *vignettes/Ancestral.Rnw
-efd64c8a6ee13e6443a20adf473d17e4 *vignettes/Networx.Rmd
+6cf134802b4b98bb286e71f135ead8bd *vignettes/Networx.Rmd
9608eda76927e9438fa12a63ef8692cd *vignettes/Trees.RData
-a6ff44c68626895b856effc7e646446d *vignettes/Trees.Rnw
+e3f16cf1e79426ba67016beed7ee6011 *vignettes/Trees.Rnw
6af5f4d4c2e93469cc19d46a97ab5d0f *vignettes/exdna.txt
1ca8b1d97a011a9d0ecd9630d548dfb3 *vignettes/movie.gif
6e60535c60981d4f5492b73ebf7344df *vignettes/phangorn-specials.Rnw
-27f308df3a021c770a6b43e4c5ad9eff *vignettes/phangorn.bib
+baba59abae9999be96bcc36d17793385 *vignettes/phangorn.bib
d3069d1eff9e70bed655b8962bf4ee2b *vignettes/primates.dna
diff --git a/NAMESPACE b/NAMESPACE
index 8560741..9f83d9a 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,20 +1,29 @@
useDynLib(phangorn, .registration=TRUE)
-importFrom(ape, as.DNAbin, as.phylo, plot.phylo, as.prop.part, nj,old2new.phylo,
- new2old.phylo, is.rooted, unroot, root, is.binary.tree, as.alignment,
+importFrom(ape, as.DNAbin, as.phylo, plot.phylo, as.prop.part, nj,old2new.phylo, dist.dna,
+ new2old.phylo, cophenetic.phylo, is.rooted, unroot, root, is.binary.tree, as.alignment,
di2multi, multi2di, .uncompressTipLabel, .compressTipLabel, prop.part,
postprocess.prop.part, speciesTree, plotPhyloCoor, cladogram.plot, phylogram.plot,
node.depth.edgelength, drop.tip, stree, rtree, is.ultrametric, .PlotPhyloEnv,
nodelabels, edgelabels, BOTHlabels, read.nexus.data, write.nexus.data,
read.dna, write.dna, reorder.phylo, dist.nodes, collapse.singles, rotate)
-importFrom(Matrix, Matrix, sparseMatrix, spMatrix, crossprod, solve)
importFrom(igraph, graph, graph.adjacency, layout.kamada.kawai, plot.igraph,
get.shortest.paths, set.edge.attribute, vcount, graph.edgelist, topological.sort)
-importFrom(stats, AIC, logLik, reorder)
+importFrom(Matrix, Matrix, sparseMatrix, spMatrix, crossprod, solve) #, %*%
#importFrom(fastmatch, fmatch)
-#importFrom(rgl, open3d, segments3d, spheres3d, rgl.texts)
-importFrom(parallel, mclapply)
importFrom(nnls, nnls)
+importFrom(parallel, mclapply)
+importFrom(quadprog, solve.QP, solve.QP.compact)
+#importFrom(rgl, open3d, segments3d, spheres3d, rgl.texts)
+#importFrom(Rcpp, evalCpp)
+importFrom(stats, AIC, logLik, reorder, update, optim, optimize, constrOptim, na.omit,
+ cophenetic, hclust, as.dist, pchisq, reshape, qgamma, pgamma, model.matrix,
+ aggregate, lm.fit, xtabs, quantile)
+importFrom(graphics, plot.new, plot.window, text, par, plot, abline, strwidth, axis, title, segments,
+ points, image, matplot, legend, hist)
+importFrom(utils, read.table, download.file, citation, stack, installed.packages, write.table)
+importFrom(grDevices, rgb, rainbow, adjustcolor)
+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,
@@ -29,12 +38,14 @@ export(pml, optim.pml, pml.control, parsimony, optim.parsimony, pratchet, NJ, UN
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, countCycles,
- SOWH.test, plot.networx, presenceAbsence, as.networx.splits, addTrivialSplits)
+ SOWH.test, plot.networx, presenceAbsence, as.networx.splits, addTrivialSplits,
+ phyDat2alignment, AICc, readDist, writeDist, discrete.gamma)
S3method('[', splits)
S3method(addConfidences, phylo)
S3method(addConfidences, splits)
S3method(addConfidences, networx)
+S3method(AICc, pml)
S3method(as.character, phyDat)
S3method(as.data.frame, phyDat)
S3method(anova, pml)
@@ -57,10 +68,13 @@ S3method(as.splits, multiPhylo)
S3method(as.splits, prop.part)
S3method(as.splits, networx)
S3method(as.prop.part, splits)
+S3method(cophenetic, networx)
+S3method(cophenetic, splits)
S3method(logLik, pml)
S3method(logLik, pmlPart)
S3method(logLik, pmlMix)
S3method(plot, pml)
+S3method(plot, pmlPart)
S3method(plot, networx)
S3method(plot, pmlCluster)
S3method(print, phyDat)
diff --git a/NEWS b/NEWS
index bdae77e..a860529 100644
--- a/NEWS
+++ b/NEWS
@@ -1,3 +1,49 @@
+ CHANGES in PHANGORN VERSION 1.99-14
+
+NEW FEATURES
+
+ o phyDat2alignment to exports files to seqinr
+
+ o readDist, writeDist to import / export distance matrices
+
+ o cophenetic distance function for splits and networx
+
+ o added unit tests
+
+BUG FIXES
+
+ o as.splits.networx did not work properly for 4 taxa
+
+ (reported by Laurélène Faye)
+
+ o RF.dist returned sometimes wrong values
+
+ (reported by Andres Dajles)
+
+ o plotBS did sometimes not work if the tree had no edge lengths
+
+ o plotBS did not work propoerly if input trees were rooted
+
+ (reported by Quynh Quach)
+
+ o plot.networx ignored cex argument in "2D" plots
+
+ o Siblings ignored include.self argument if node is a vector
+
+OTHER CHANGES
+
+ o plotBS got an additional argument p to plot only support values
+
+ greater than p
+
+ o pml and optim.pml now uses more C-code (and is a bit faster)
+
+ o defaults in modelTest changed
+
+ o discrete.gamma is now exported in the NAMESPACE
+
+
+
CHANGES in PHANGORN VERSION 1.99-13
OTHER CHANGES
diff --git a/R/clanistic.R b/R/clanistic.R
index 0e28944..4d44847 100644
--- a/R/clanistic.R
+++ b/R/clanistic.R
@@ -341,7 +341,7 @@ compareSplits <- function(res, nam1, nam2){
splits = attr(res, "Perfect")
dat = attr(attr(res, "Perfect"), "data")
res = matrix(NA, length(ntrees), length(nam1)*length(nam2))
- for(m in 1:length(trees)){
+ for(m in 1:ntrees){
k=1
trnam=ntrees[m]
for(i in nam1){
diff --git a/R/distSeq.R b/R/distSeq.R
index 8ec34c2..ce8a825 100644
--- a/R/distSeq.R
+++ b/R/distSeq.R
@@ -191,3 +191,22 @@ dist.logDet = function (x)
return(d)
}
+
+readDist <- function(file){ #, format="phylip"
+ tmp <- read.table(file, skip=1, stringsAsFactors = FALSE)
+ labels = tmp[,1]
+ dm <- as.matrix(tmp[,-1])
+ dimnames(dm)=list(labels, labels)
+ as.dist(dm)
+}
+
+
+writeDist <- function(dm, file=""){ # , format="phylip"
+ dm <- as.matrix(dm)
+ cat(ncol(dm), "\n", file=file)
+ write.table(dm, file, append=TRUE, quote=FALSE, col.names=FALSE)
+}
+
+
+
+
diff --git a/R/distTree.R b/R/distTree.R
new file mode 100644
index 0000000..38d0142
--- /dev/null
+++ b/R/distTree.R
@@ -0,0 +1,554 @@
+#
+# UPGMA, NJ, UNJ, nnls
+#
+"upgma" <- function(D,method="average",...){
+ DD=as.dist(D)
+ hc = hclust(DD,method=method,...)
+ result = as.phylo(hc)
+ result = reorder(result, "postorder")
+ result
+}
+
+
+"wpgma" <- function(D,method="mcquitty",...){
+ DD=as.dist(D)
+ hc = hclust(DD,method=method,...)
+ result = as.phylo(hc)
+ result = reorder(result, "postorder")
+ result
+}
+
+
+NJ_old <- function(x)
+{
+ x = as.matrix(x)
+ labels <- attr(x, "Labels")[[1]]
+ edge.length = NULL
+ edge = NULL
+ d = as.matrix(x)
+ if (is.null(labels))
+ labels = colnames(d)
+ l = dim(d)[1]
+ m = l - 2
+ nam = 1L:l
+ k = 2L * l - 2L
+ while (l > 2) {
+ r = rowSums(d)/(l - 2)
+ i = 0
+ j = 0
+ tmp <- .C("out", as.double(d), as.double(r), as.integer(l),
+ as.integer(i), as.integer(j))
+ e2 = tmp[[5]]
+ e1 = tmp[[4]]
+ l1 = d[e1, e2]/2 + (r[e1] - r[e2])/(2)
+ l2 = d[e1, e2] - l1
+ edge.length = c(l1, l2, edge.length)
+ edge = rbind(c(k, nam[e2]), edge)
+ edge = rbind(c(k, nam[e1]), edge)
+ nam = c(nam[c(-e1, -e2)], k)
+ dnew = (d[e1, ] + d[e2, ] - d[e1, e2])/2
+ d = cbind(d, dnew)
+ d = rbind(d, c(dnew, 0))
+ d = d[-c(e1, e2), -c(e1, e2)]
+ k = k - 1L
+ l = l - 1L
+ }
+ edge.length = c(d[2, 1], edge.length)
+ attr(edge.length,"names") = NULL
+ result = list(edge = rbind(c(nam[2], nam[1]), edge), edge.length = edge.length,
+ tip.label = labels, Nnode = m)
+ class(result) <- "phylo"
+ reorder(result, "postorder")
+}
+
+
+NJ <- function(x) reorder(nj(x), "postorder")
+
+
+UNJ <- function(x)
+{
+ x = as.matrix(x)
+ labels <- attr(x, "Labels")[[1]]
+ edge.length = NULL
+ edge = NULL
+ d = as.matrix(x)
+ if (is.null(labels))
+ labels = colnames(d)
+ l = dim(d)[1]
+ n = l
+ nam = as.character(1:l)
+ m=l-2
+ nam = 1:l
+ k = 2*l-2
+ w = rep(1,l)
+ while (l > 2) {
+ r = rowSums(d)/(l - 2)
+ i = 0
+ j = 0
+ tmp <- .C("out", as.double(d), as.double(r), as.integer(l), as.integer(i), as.integer(j))
+ e2 = tmp[[5]]
+ e1 = tmp[[4]]
+ l1 = d[e1, e2]/2 + sum((d[e1,-c(e1,e2)] - d[e2,-c(e1,e2)])*w[-c(e1,e2)])/(2*(n-w[e1]-w[e2]))
+ l2 = d[e1, e2]/2 + sum((d[e2,-c(e1,e2)] - d[e1,-c(e1,e2)])*w[-c(e1,e2)])/(2*(n-w[e1]-w[e2]))
+ edge.length = c(l1, l2, edge.length)
+ edge = rbind(c(k, nam[e2]), edge)
+ edge = rbind(c(k, nam[e1]), edge)
+ nam = c(nam[c(-e1, -e2)], k)
+
+ dnew = (w[e1]*d[e1, ] + w[e2]*d[e2, ] - w[e1]*l1 - w[e2]*l2)/(w[e1] + w[e2])
+ d = cbind(d, dnew)
+ d = rbind(d, c(dnew, 0))
+ d = d[-c(e1, e2), -c(e1, e2)]
+ w = c(w, w[e1] + w[e2])
+ w = w[-c(e1, e2)]
+ k = k - 1
+ l = l - 1
+ }
+ edge.length=c(d[2,1],edge.length)
+ result = list(edge = rbind(c(nam[2], nam[1]), edge),
+ edge.length=edge.length, tip.label = labels, Nnode=m)
+ class(result) <- "phylo"
+ reorder(result)
+}
+
+
+PNJ <- function (data)
+{
+ q <- l <- r <- length(data)
+ weight <- attr(data,"weight")
+
+ height = NULL
+ parentNodes <- NULL
+ childNodes <- NULL
+ nam <- names(data)
+ tip.label <- nam
+ edge = 1:q
+
+ z = 0
+ D = matrix(0, q, q)
+
+ for (i in 1:(l - 1)) {
+ for (j in (i + 1):l) {
+ w = (data[[i]] * data[[j]]) %*% c(1, 1, 1, 1)
+ D[i, j] = sum(weight[w==0])
+ }
+ }
+
+ while (l > 1) {
+ l = l - 1
+ z = z + 1
+ d = D + t(D)
+ if(l>1) r = rowSums(d)/(l-1)
+ if(l==1) r = rowSums(d)
+ M = d - outer(r,r,"+")
+ diag(M) = Inf
+
+ e=which.min(M)
+ e0=e%%length(r)
+ e1 = ifelse(e0==0, length(r), e0)
+ e2= ifelse(e0==0, e%/%length(r), e%/%length(r) + 1)
+
+ ind = c(e1,e2)
+ len = d[e]/2
+ nam = c(nam[-ind], as.character(-l))
+
+ parentNodes = c(parentNodes,-l,-l)
+ childNodes = c(childNodes,edge[e1],edge[e2])
+
+ height = c(height, len, len)
+ edge = c(edge[-ind], -l)
+ w = (data[[e1]] * data[[e2]]) %*% c(1, 1, 1, 1)
+ w = which(w == 0)
+ newDat = data[[e1]] * data[[e2]]
+ newDat[w, ] = data[[e1]][w, ] + data[[e2]][w, ]
+ data = data[-c(e1,e2)]
+ data[[l]] = newDat
+ if (l > 1) {
+ D = as.matrix(D[, -ind])
+ D = D[-ind, ]
+ dv = numeric(l - 1)
+ for (i in 1:(l - 1)) {
+ w = (data[[i]] * data[[l]]) %*% c(1, 1, 1, 1)
+ dv[i] = sum(weight[w==0])
+ }
+ D = cbind(D, dv)
+ D = rbind(D, 0)
+ }
+ }
+ tree <- list(edge = cbind(as.character(parentNodes),as.character(childNodes)),tip.label=tip.label)
+ class(tree) <- "phylo"
+ tree <- old2new.phylo(tree)
+ reorder(tree)
+}
+
+
+#
+# Distance Matrix methods
+#
+
+# as.Matrix, sparse = TRUE,
+designTree <- function(tree, method="unrooted", sparse=FALSE, ...){
+ if (!is.na(pmatch(method, "all")))
+ method <- "unrooted"
+ METHOD <- c("unrooted", "rooted")
+ method <- pmatch(method, METHOD)
+ if (is.na(method)) stop("invalid method")
+ if (method == -1) stop("ambiguous method")
+ if(!is.rooted(tree) & method==2) stop("tree has to be rooted")
+ if(method==1){ X <- designUnrooted(tree,...)
+ if(sparse) X = Matrix(X)
+ }
+ if(method==2) X <- designUltra(tree, sparse=sparse,...)
+ X
+}
+
+
+# splits now work
+designUnrooted = function (tree, order = NULL)
+{
+ if(inherits(tree, "phylo")){
+ if (is.rooted(tree))
+ tree = unroot(tree)
+ p = bipartition(tree)
+ }
+ if(inherits(tree, "splits")) p <- as.matrix(tree)
+ if (!is.null(order))
+ p = p[, order]
+
+ m = dim(p)[2]
+ ind = rowSums(p)
+ p=p[ind!=m,]
+ n = dim(p)[1]
+ res = matrix(0, (m - 1) * m/2, n)
+ k = 1
+ for (i in 1:(m - 1)) {
+ for (j in (i + 1):m) {
+ res[k, ] = p[, i] != p[, j]
+ k = k + 1
+ }
+ }
+ if(inherits(tree, "phylo"))colnames(res) = paste(tree$edge[, 1], tree$edge[, 2], sep = "<->")
+ res
+}
+
+
+designUltra <- function (tree, sparse=TRUE)
+{
+ if (is.null(attr(tree, "order")) || attr(tree, "order") == "cladewise")
+ tree = reorder(tree, "postorder")
+ leri = allChildren(tree)
+ bp = bip(tree)
+ n = length(tree$tip)
+ l = tree$Nnode
+ nodes = integer(l)
+ k = 1L
+ u=numeric( n * (n - 1)/2)
+ v=numeric( n * (n - 1)/2)
+ m = 1L
+ for (i in 1:length(leri)) {
+ if (!is.null(leri[[i]])) {
+ if(length(leri[[i]])==2)ind = getIndex(bp[[leri[[i]][1] ]], bp[[leri[[i]][2] ]], n)
+ else {
+ ind=NULL
+ le=leri[[i]]
+ nl = length(le)
+ for(j in 1:(nl-1)) ind =c(ind, getIndex(bp[[le[j] ]], unlist(bp[ le[(j+1):nl] ]), n))
+ }
+ li = length(ind)
+ v[m: (m+li-1)]=k
+ u[m: (m+li-1)]=ind
+ nodes[k]=i
+ m = m+li
+ k = k + 1L
+ }
+ }
+ if(sparse) X = sparseMatrix(i=u,j=v, x=2L)
+ else{
+ X = matrix(0L, n * (n - 1)/2, l)
+ X[cbind(u,v)]=2L
+ }
+ colnames(X) = nodes
+ attr(X, "nodes") = nodes
+ X
+}
+
+
+designUnrooted2 <- function (tree, sparse=TRUE)
+{
+ if (is.null(attr(tree, "order")) || attr(tree, "order") == "cladewise")
+ tree = reorder(tree, "postorder")
+ leri = allChildren(tree)
+ bp = bip(tree)
+ n = length(tree$tip)
+ l = tree$Nnode
+ nodes = integer(l)
+ nTips = as.integer(length(tree$tip))
+ k = nTips
+ u=numeric( n * (n - 1)/2)
+ v=numeric( n * (n - 1)/2)
+ z=numeric( n * (n - 1)/2)
+ y=numeric( n * (n - 1)/2)
+ p=1L
+ m = 1L
+ for (i in 1:length(leri)) {
+ if (!is.null(leri[[i]])) {
+
+ 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] ]]))
+ ztmp = rep(bp[[leri[[i]][2] ]], length(bp[[leri[[i]][1] ]]))
+ }
+ else {
+ # browser()
+ ind=NULL
+ le=leri[[i]]
+ nl = length(le)
+ ytmp=NULL
+ ztmp=NULL
+ for(j in 1:(nl-1)){
+ bp1 = bp[[le[j] ]]
+ bp2 = unlist(bp[le[(j+1):nl] ])
+ ind =c(ind, getIndex(bp1, unlist(bp2), n))
+ ytmp = c(ytmp, rep(bp1, each = length(bp2)))
+ ztmp = c(ztmp, rep(bp2, length(bp1)))
+ }
+ }
+
+ li = length(ind)
+ v[m: (m+li-1)]=k
+ u[m: (m+li-1)]=ind
+ y[m: (m+li-1)]=ytmp
+ z[m: (m+li-1)]=ztmp
+
+ nodes[p]=i
+ m = m+li
+ k = k + 1L
+ p = p + 1L
+ }
+ }
+ jj = c(y,z) #[ind],v)
+ ii = c(u,u) #[ind],u)
+ ind = (jj < nTips)
+ jj = c(jj[ind], v)
+ ii = c(ii[ind], u)
+ l1 = length(u)
+ l2 = sum(ind)
+ x= rep(c(-1L,2L), c(l2, l1))
+
+ X = sparseMatrix(i=ii,j=jj, x=x)
+ if(!sparse){
+ X = as.matrix(X)
+ }
+ nodes = c(1:(nTips-1L), nodes)
+ colnames(X) = nodes
+ attr(X, "nodes") = nodes
+ X
+}
+
+
+nnls.tree <- function(dm, tree, rooted=FALSE, trace=1){
+ if(is.rooted(tree) & rooted==FALSE){
+ tree = unroot(tree)
+ warning("tree was rooted, I unrooted the tree!")
+ }
+ tree = reorder(tree, "postorder")
+ dm = as.matrix(dm)
+ k = dim(dm)[1]
+ labels = tree$tip
+ dm = dm[labels,labels]
+ y = dm[lower.tri(dm)]
+ #computing the design matrix from the tree
+ if(rooted) X = designUltra(tree)
+ else X = designUnrooted2(tree)
+
+ lab = attr(X, "nodes")
+
+ # na.action
+ if(any(is.na(y))){
+ ind = which(is.na(y))
+ X = X[-ind,,drop=FALSE]
+ y= y[-ind]
+ }
+ # LS solution
+ Dmat <- crossprod(X) # cross-product computations
+ dvec <- crossprod(X, y)
+ betahat <- as.vector(solve(Dmat, dvec))
+
+ bhat = numeric(max(tree$edge))
+ bhat[as.integer(lab)] = betahat
+ betahat = bhat[tree$edge[,1]] - bhat[tree$edge[,2]]
+
+ if(!any(betahat<0)){
+ if(!rooted){
+ RSS = sum((y-(X%*%betahat))^2)
+ if(trace)print(paste("RSS:", RSS))
+ attr(tree, "RSS") = RSS
+ }
+ tree$edge.length = betahat
+ return(tree)
+ }
+
+ # 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
+
+ # quadratic programing solving
+ 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]]
+ tree$edge.length = betahat
+ tree
+}
+
+
+nnls.phylo <- function(x, dm, rooted=FALSE, trace=0){
+ nnls.tree(dm, x, rooted, trace=trace)
+}
+
+
+nnls.splits <- 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)
+ 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)
+ 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]
+ }
+
+ Dmat <- crossprod(X) # cross-product computations
+ dvec <- crossprod(X, y)
+ betahat <- as.vector(solve(Dmat, dvec))
+
+ if(!any(betahat<0)){
+ RSS = sum((y-(X%*%betahat))^2)
+ if(trace)print(paste("RSS:", RSS))
+ attr(x, "RSS") = RSS
+ attr(x, "weights") = betahat
+ return(x)
+ }
+ n = dim(X)[2]
+
+ 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]
+ attr(x, "weights") = betahat[ind]
+ if(trace)print(paste("RSS:", RSS))
+ attr(x, "RSS") = RSS
+ x
+}
+
+nnls.networx <- function(x, dm){
+ spl <- attr(x, "splits")
+ spl2 <- nnls.splits(spl, dm)
+ weight <- attr(spl, "weight")
+ weight[] <- 0
+ weight[match(spl2, spl)] = attr(spl2, "weight")
+ attr(attr(x, "splits"), "weight") <- weight
+ x$edge.length = weight[x$splitIndex]
+ x
+}
+
+
+designSplits <- function (x, splits = "all", ...)
+{
+ if (!is.na(pmatch(splits, "all")))
+ splits <- "all"
+ if(inherits(x, "splits")) return(designUnrooted(x))
+ SPLITS <- c("all", "star") #,"caterpillar")
+ splits <- pmatch(splits, SPLITS)
+ 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)
+ return(X)
+}
+
+# add return splits=FALSE
+designAll <- function(n, add.split=FALSE){
+ Y = matrix(0L, n*(n-1)/2, n)
+ k = 1
+ for(i in 1:(n-1)){
+ for(j in (i+1):n){
+ Y[k,c(i,j)]=1L
+ k=k+1L
+ }
+ }
+ m <- n-1L
+ X <- matrix(0L, m+1, 2^m)
+ for(i in 1:m)
+ X[i, ] <- rep(rep(c(0L,1L), each=2^(i-1)),2^(m-i))
+ X <- X[,-1]
+ if(!add.split) return((Y%*%X)%%2)
+ list(X=(Y%*%X)%%2,Splits=t(X))
+}
+
+
+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)))
+ res
+}
+
+
diff --git a/R/fitch.R b/R/fitch.R
index de52340..16bae9e 100644
--- a/R/fitch.R
+++ b/R/fitch.R
@@ -172,7 +172,7 @@ fitch.spr <- function(tree, data){
tree
}
-# raus
+# raus oder richten
fitch.spr2 <- function(tree, data){
nTips = as.integer(length(tree$tip))
nr = attr(data, "nr")
diff --git a/R/modelTest.R b/R/modelTest.R
new file mode 100644
index 0000000..2fae133
--- /dev/null
+++ b/R/modelTest.R
@@ -0,0 +1,201 @@
+modelTest2 <- function (object, tree = NULL, model = c("JC", "F81", "K80",
+ "HKY", "SYM", "GTR"), G = TRUE, I = TRUE, k = 4, freq=FALSE,
+ control = pml.control(epsilon = 1e-08, maxit = 10, trace = 1),
+ multicore = FALSE, mc.cores = getOption("mc.cores", 1L))
+{
+ if (class(object) == "phyDat")
+ data = object
+ if (class(object) == "pml") {
+ data = object$data
+ if (is.null(tree))
+ tree = object$tree
+ }
+
+ if(attr(data, "type")=="DNA") type = 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(attr(data, "type")=="AA") type = .aamodels
+ model = match.arg(model, type, TRUE)
+
+ env = new.env()
+ assign("data", data, envir=env)
+
+ if (is.null(tree))
+ tree = NJ(dist.hamming(data))
+ else{
+ tree <- nnls.phylo(tree, dist.ml(data))
+ # may need something faster for trees > 500 taxa
+ }
+ trace <- control$trace
+ control$trace = trace - 1
+ fit = pml(tree, data)
+ fit = optim.pml(fit, control = control)
+ l = length(model)
+ if(attr(fit$data, "type")=="DNA")freq=FALSE
+ n = 1L + sum(I + G + (G & I) + freq + (freq & I) + (freq & G) + (freq & G & I))
+ nseq = sum(attr(data, "weight"))
+
+
+
+ fitPar = function(model, fit, G, I, k, freq) {
+ m = 1
+ res = matrix(NA, n, 6)
+ res = as.data.frame(res)
+ colnames(res) = c("Model", "df", "logLik", "AIC", "AICc", "BIC")
+ data.frame(c("Model", "df", "logLik", "AIC", "AICc", "BIC"))
+ calls = vector("list", n)
+ trees = vector("list", n)
+ fittmp = optim.pml(fit, model = model, control = control)
+ res[m, 1] = model
+ res[m, 2] = fittmp$df
+ res[m, 3] = fittmp$logLik
+ res[m, 4] = AIC(fittmp)
+ res[m, 5] = AICc(fittmp)
+ res[m, 6] = AIC(fittmp, k = log(nseq))
+ calls[[m]] = fittmp$call
+
+ trees[[m]] = fittmp$tree
+ m = m + 1
+ if (I) {
+ if(trace>0)print(paste(model, "+I", sep = ""))
+ fitI = optim.pml(fittmp, model = model, optInv = TRUE,
+ control = control)
+ res[m, 1] = paste(model, "+I", sep = "")
+ res[m, 2] = fitI$df
+ res[m, 3] = fitI$logLik
+ res[m, 4] = AIC(fitI)
+ res[m, 5] = AICc(fitI)
+ res[m, 6] = AIC(fitI, k = log(nseq))
+ calls[[m]] = fitI$call
+ trees[[m]] = fitI$tree
+ m = m + 1
+ }
+ if (G) {
+ if(trace>0)print(paste(model, "+G", sep = ""))
+ fitG = update(fittmp, k = k)
+ fitG = optim.pml(fitG, model = model, optGamma = TRUE,
+ control = control)
+ res[m, 1] = paste(model, "+G", sep = "")
+ res[m, 2] = fitG$df
+ res[m, 3] = fitG$logLik
+ res[m, 4] = AIC(fitG)
+ res[m, 5] = AICc(fitG)
+ res[m, 6] = AIC(fitG, k = log(nseq))
+ calls[[m]] = fitG$call
+ trees[[m]] = fitG$tree
+ m = m + 1
+ }
+ if (G & I) {
+ if(trace>0)print(paste(model, "+G+I", sep = ""))
+ fitGI = update(fitI, k = k)
+ fitGI = optim.pml(fitGI, model = model, optGamma = TRUE,
+ optInv = TRUE, control = control)
+ res[m, 1] = paste(model, "+G+I", sep = "")
+ res[m, 2] = fitGI$df
+ res[m, 3] = fitGI$logLik
+ res[m, 4] = AIC(fitGI)
+ res[m, 5] = AICc(fitGI)
+ res[m, 6] = AIC(fitGI, k = log(nseq))
+ calls[[m]] = fitGI$call
+ trees[[m]] = fitGI$tree
+ m = m + 1
+ }
+ if (freq) {
+ if(trace>0)print(paste(model, "+F", sep = ""))
+ fitF = optim.pml(fittmp, model = model, optBf = TRUE,
+ control = control)
+ res[m, 1] = paste(model, "+F", sep = "")
+ res[m, 2] = fitF$df
+ res[m, 3] = fitF$logLik
+ res[m, 4] = AIC(fitF)
+ res[m, 5] = AICc(fitF)
+ res[m, 6] = AIC(fitF, k = log(nseq))
+ calls[[m]] = fitF$call
+ trees[[m]] = fitF$tree
+ m = m + 1
+ }
+ if (freq & I) {
+ if(trace>0)print(paste(model, "+I+F", sep = ""))
+ fitIF <- update(fitF, inv = fitI$inv)
+ fitIF = optim.pml(fitIF, model = model, optBf = TRUE, optInv = TRUE,
+ control = control)
+ res[m, 1] = paste(model, "+I+F", sep = "")
+ res[m, 2] = fitIF$df
+ res[m, 3] = fitIF$logLik
+ res[m, 4] = AIC(fitIF)
+ res[m, 5] = AICc(fitIF)
+ res[m, 6] = AIC(fitIF, k = log(nseq))
+ calls[[m]] = fitIF$call
+ trees[[m]] = fitIF$tree
+ m = m + 1
+ }
+ if (freq & G) {
+ if(trace>0)print(paste(model, "+G+F", sep = ""))
+ fitGF <- update(fitF, k=k, shape=fitG$shape)
+ fitGF = optim.pml(fitGF, model = model, optBf = TRUE, optGamma = TRUE,
+ control = control)
+ res[m, 1] = paste(model, "+G+F", sep = "")
+ res[m, 2] = fitGF$df
+ res[m, 3] = fitGF$logLik
+ res[m, 4] = AIC(fitGF)
+ res[m, 5] = AICc(fitGF)
+ res[m, 6] = AIC(fitGF, k = log(nseq))
+ calls[[m]] = fitGF$call
+ trees[[m]] = fitGF$tree
+ m = m + 1
+ }
+ if (freq & G & I) {
+ if(trace>0)print(paste(model, "+G+I+F", sep = ""))
+ fitGIF <- update(fitIF, k=k)
+ fitGIF = optim.pml(fitGIF, model = model, optBf = TRUE, optInv = TRUE, , optGamma = TRUE,
+ control = control)
+ res[m, 1] = paste(model, "+G+I+F", sep = "")
+ res[m, 2] = fitGIF$df
+ res[m, 3] = fitGIF$logLik
+ res[m, 4] = AIC(fitGIF)
+ res[m, 5] = AICc(fitGIF)
+ res[m, 6] = AIC(fitGIF, k = log(nseq))
+ calls[[m]] = fitGIF$call
+ trees[[m]] = fitGIF$tree
+ m = m + 1
+ }
+
+
+ list(res, trees, calls)
+ }
+ eval.success <- FALSE
+ if (!eval.success & multicore) {
+ # !require(parallel) ||
+ # if (.Platform$GUI != "X11") {
+ # warning("package 'parallel' not found or GUI is used, \n analysis is performed in serial")
+ # }
+ # else {
+ RES <- mclapply(model, fitPar, fit, G, I, k, freq, mc.cores=mc.cores)
+ eval.success <- TRUE
+ # }
+ }
+ if (!eval.success)
+ RES <- lapply(model, fitPar, fit, G, I, k, freq)
+ # res <- RES <- lapply(model, fitPar, fit, G, I, k, freq)
+ RESULT = matrix(NA, n * l, 6)
+ RESULT = as.data.frame(RESULT)
+ colnames(RESULT) = c("Model", "df", "logLik", "AIC", "AICc", "BIC")
+ for (i in 1:l) RESULT[((i - 1) * n + 1):(n * i), ] = RES[[i]][[1]]
+ for(i in 1:l){
+ for(j in 1:n){
+ mo = RES[[i]][[1]][j,1]
+ tname = paste("tree_", mo, sep = "")
+ tmpmod = RES[[i]][[3]][[j]]
+ tmpmod["tree"] = call(tname)
+ if(!is.null(tmpmod[["k"]]))tmpmod["k"] = k
+ if(attr(data, "type")=="AA") tmpmod["model"] = RES[[i]][[1]][1,1]
+ assign(tname, RES[[i]][[2]][[j]], envir=env)
+ assign(mo, tmpmod, envir=env)
+ }
+ }
+ attr(RESULT, "env") = env
+ RESULT
+}
+
+
diff --git a/R/neighborNet.R b/R/neighborNet.R
index a2d3fe5..023d295 100644
--- a/R/neighborNet.R
+++ b/R/neighborNet.R
@@ -10,9 +10,11 @@ cyclicSplits <- 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>4){
+ 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
@@ -21,6 +23,7 @@ cyclicSplits <- function(k, labels=NULL){
}
if(is.null(labels)) labels=(as.character(1:k))
attr(res, 'labels') =labels
+ attr(res, "cycle") = 1:k
class(res)="splits"
res
}
diff --git a/R/networx.R b/R/networx.R
index 6b5eb6f..896bbb8 100644
--- a/R/networx.R
+++ b/R/networx.R
@@ -337,7 +337,6 @@ compatible3 <- function(x, y=NULL)
ylabels = attr(y, "labels")
if(identical(xlabels, ylabels)) labels = xlabels
else labels = intersect(xlabels, ylabels)
-# if(length(labels) maybe warning
nx = length(x)
ny = length(y)
bp1 = as.matrix(x)[,labels, drop=FALSE]
@@ -993,7 +992,10 @@ plot.networx = function(x, type="3D", use.edge.length = TRUE, show.tip.label=TRU
if(show.tip.label)node.label[1:nTips] = ""
chk <- FALSE
- if(type=="3D") chk <- .check.pkg("rgl")
+
+
+
+ if(type=="3D") chk <- requireNamespace("rgl", quietly = TRUE) #.check.pkg("rgl")
if(!chk && type=="3D"){
warning("type=\"3D\" requires the package \"rgl\"\n, plotting =\"2D\" instead!\n")
type="2D"
@@ -1020,9 +1022,14 @@ plotRGL <- function(coords, net, show.tip.label=TRUE,
show.nodes=FALSE, tip.color = "blue", edge.color="grey",
edge.width = 3, font = 3, cex = par("cex"), ...){
- chk <- .check.pkg("rgl")
- if(!chk) open3d <- segments3d <- spheres3d <- rgl.texts <- function(...)NULL
-
+# chk <- .check.pkg("rgl")
+# if(!chk) open3d <- segments3d <- spheres3d <- rgl.texts <- function(...)NULL
+
+ open3d <- rgl::open3d
+ segments3d <- rgl::segments3d
+ spheres3d <- rgl::spheres3d
+ rgl.texts <- rgl::rgl.texts
+
edge = net$edge
x = coords[,1]
@@ -1063,7 +1070,7 @@ plot2D <- function(coords, net, show.tip.label=TRUE,
yy = coords[,2]
nTips = length(label)
- cex=1
+# cex=1
xlim <- range(xx)
ylim <- range(yy)
@@ -1301,7 +1308,6 @@ read.nexus.splits <- function(file)
tmp = X[i]
tmp = sub("\\s+", "", tmp)
tmp = strsplit(tmp, "\t")[[1]]
- if(length(tmp)!=(mformat+1)) warning("blub")
if(flab)labels[j] = as.numeric(tmp[ind[1]])
if(fwei)weights[j] = as.numeric(tmp[ind[2]])
if(fcon)confidences[j] = as.numeric(tmp[ind[3]])
diff --git a/R/parsimony.R b/R/parsimony.R
index c16adb1..ebab78c 100644
--- a/R/parsimony.R
+++ b/R/parsimony.R
@@ -261,8 +261,6 @@ lowerBound <- function(x, cost=NULL){
if(is.null(cost))cost <- 1 - diag(nc)
-
-
if(any(ust>1)){
ust = ust[ust>1]
m <- max(ust)
@@ -287,7 +285,7 @@ upperBound <- function(x, cost=NULL){
CI <- function (tree, data, cost=NULL){
- pscore = sankoff(tree, data, cost=cost)
+ pscore = sankoffNew(tree, data, cost=cost)
weight = attr(data, "weight")
data = subset(data, tree$tip.label)
m = lowerBound(data, cost=cost)
@@ -429,7 +427,7 @@ old2new.phyDat <- function(obj){
X = matrix(rep(rowSums(contrast), each=nr),nrow=nr)
res <- vector("list", l)
for(i in 1:l){
- browser()
+# browser()
tmp = X - tcrossprod(obj[[i]], contrast)
res[[i]] = unlist(apply(tmp, 1, function(x)which(x<1e-6)[1]))
}
@@ -643,7 +641,7 @@ optim.parsimony <- function(tree,data, method='fitch', cost=NULL, trace=1, rearr
}
-pratchet <- function (data, start=NULL, method="fitch", maxit=100, k=5, trace=1, all=FALSE, rearrangements="SPR", ...)
+pratchet <- function (data, start=NULL, method="fitch", maxit=1000, k=10, trace=1, all=FALSE, rearrangements="SPR", ...)
{
eps = 1e-08
# if(method=="fitch" && (is.null(attr(data, "compressed")) || attr(data, "compressed") == FALSE))
@@ -762,7 +760,7 @@ 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")
+ tree <- reorder(tree, "pruningwise")
# if (!is.binary.tree(tree))
# stop("Tree must be binary!")
tmp = fitch(tree, data, site = "data")
diff --git a/R/phyDat.R b/R/phyDat.R
index 40141d2..bb93cf5 100644
--- a/R/phyDat.R
+++ b/R/phyDat.R
@@ -26,11 +26,14 @@ phyDat.default <- function (data, levels = NULL, return.index = TRUE, contrast =
if (is.matrix(data))
nam = row.names(data)
else nam = names(data)
+ if(is.null(nam))stop("data object must contain taxa names")
if (class(data) == "DNAbin")
data = as.character(data)
if (is.matrix(data))
data = as.data.frame(t(data), stringsAsFactors = FALSE)
+ if (is.vector(data))data = as.data.frame(t(data), stringsAsFactors = FALSE)
else data = as.data.frame(data, stringsAsFactors = FALSE)
+ if(length(data[[1]])==1) compress=FALSE
if(compress){
ddd = fast.table(data)
data = ddd$data
@@ -63,22 +66,18 @@ phyDat.default <- function (data, levels = NULL, return.index = TRUE, contrast =
contrast = rbind(contrast, matrix(1, k, l))
}
}
- d = dim(data)
att = attributes(data)
- data = match(unlist(data), all.levels)
- attr(data, "dim") = d
- data = as.data.frame(data, stringsAsFactors=FALSE)
+ data = lapply(data, match, all.levels) # avoid unlist
attributes(data) = att
-
+
row.names(data) = as.character(1:p)
data = na.omit(data)
-
+
aaa = match(index, attr(data, "na.action"))
index = index[is.na(aaa)]
index = match(index, unique(index))
rn = as.numeric(rownames(data))
- attr(data, "na.action") = NULL
-
+ attr(data, "na.action") = NULL
weight = weight[rn]
p = dim(data)[1]
names(data) = nam
@@ -328,7 +327,20 @@ as.phyDat.alignment <- function (x, type="DNA",...)
}
-as.alignment.phyDat <- function(x, ...) as.alignment(as.character(x))
+#as.alignment.phyDat <- function(x, ...) as.alignment(as.character(x))
+
+phyDat2alignment <- function(x){
+ z = as.character(x)
+ nam = rownames(z)
+ type = attr(x, "type")
+ seq <- switch(type,
+ DNA = tolower(apply(z, 1, paste, collapse="")),
+ AA = toupper(apply(z, 1, paste, collapse="")))
+ names(seq) <- NULL
+ res <- list(nb=length(seq), nam=nam, seq=seq, com=NA)
+ class(res) = "alignment"
+ res
+}
as.phyDat.matrix <- function (x, ...) phyDat(data=x, ...)
@@ -766,13 +778,23 @@ unique.phyDat <- function(x, incomparables=FALSE, ...) getCols(x, !duplicated(x)
allSitePattern <- function(n,levels=c("a","c","g","t"), names=NULL){
l=length(levels)
- X=matrix(0, l^n,n)
+ X=vector("list", n)
+ if(is.null(names))names(X) = paste("t",1:n, sep="")
+ else names(X)=names
for(i in 1:n)
- X[, i] = rep(rep(c(1:l), each=l^(i-1)),l^(n-i))
- for(i in 1:l)X[X==i] = levels[i]
- if(is.null(names))colnames(X) = paste("t",1:n, sep="")
- else colnames(X)=names
- phyDat.default(t(X), levels)
+ X[[i]] = rep(rep(levels, each=l^(i-1)),l^(n-i))
+ X = as.data.frame(X)
+ phyDat.default(X, levels, compress=FALSE)
+}
+
+
+constSitePattern <- function(n,levels=c("a","c","g","t"), names=NULL){
+ l=length(levels)
+ X=matrix(0, l,n)
+ X = matrix(rep(levels, each=n), n, l)
+ if(is.null(names))rownames(X) = paste("t",1:n, sep="")
+ else rownames(X)=names
+ phyDat.default(X, levels)
}
@@ -935,3 +957,9 @@ read.aa <- function (file, format = "interleaved", skip = 0, nlines = 0,
obj
}
+
+genlight2phyDat <- function(x, ambiguity=NA){
+ tmp <- as.matrix(x)
+ lev <- na.omit(unique(as.vector(tmp)))
+ phyDat(tmp, "USER", levels=lev, ambiguity=ambiguity)
+}
diff --git a/R/phylo.R b/R/phylo.R
index 2f47dc5..6eb75da 100644
--- a/R/phylo.R
+++ b/R/phylo.R
@@ -1,189 +1,4 @@
#
-# UPGMA, NJ and UNJ
-#
-"upgma" <- function(D,method="average",...){
- DD=as.dist(D)
- hc = hclust(DD,method=method,...)
- result = as.phylo(hc)
- result = reorder(result, "postorder")
- result
-}
-
-
-"wpgma" <- function(D,method="mcquitty",...){
- DD=as.dist(D)
- hc = hclust(DD,method=method,...)
- result = as.phylo(hc)
- result = reorder(result, "postorder")
- result
-}
-
-
-NJ_old <- function(x)
-{
- x = as.matrix(x)
- labels <- attr(x, "Labels")[[1]]
- edge.length = NULL
- edge = NULL
- d = as.matrix(x)
- if (is.null(labels))
- labels = colnames(d)
- l = dim(d)[1]
- m = l - 2
- nam = 1L:l
- k = 2L * l - 2L
- while (l > 2) {
- r = rowSums(d)/(l - 2)
- i = 0
- j = 0
- tmp <- .C("out", as.double(d), as.double(r), as.integer(l),
- as.integer(i), as.integer(j))
- e2 = tmp[[5]]
- e1 = tmp[[4]]
- l1 = d[e1, e2]/2 + (r[e1] - r[e2])/(2)
- l2 = d[e1, e2] - l1
- edge.length = c(l1, l2, edge.length)
- edge = rbind(c(k, nam[e2]), edge)
- edge = rbind(c(k, nam[e1]), edge)
- nam = c(nam[c(-e1, -e2)], k)
- dnew = (d[e1, ] + d[e2, ] - d[e1, e2])/2
- d = cbind(d, dnew)
- d = rbind(d, c(dnew, 0))
- d = d[-c(e1, e2), -c(e1, e2)]
- k = k - 1L
- l = l - 1L
- }
- edge.length = c(d[2, 1], edge.length)
- attr(edge.length,"names") = NULL
- result = list(edge = rbind(c(nam[2], nam[1]), edge), edge.length = edge.length,
- tip.label = labels, Nnode = m)
- class(result) <- "phylo"
- reorder(result, "postorder")
-}
-
-
-NJ <- function(x) reorder(nj(x), "postorder")
-
-
-UNJ <- function(x)
-{
- x = as.matrix(x)
- labels <- attr(x, "Labels")[[1]]
- edge.length = NULL
- edge = NULL
- d = as.matrix(x)
- if (is.null(labels))
- labels = colnames(d)
- l = dim(d)[1]
- n = l
- nam = as.character(1:l)
- m=l-2
- nam = 1:l
- k = 2*l-2
- w = rep(1,l)
- while (l > 2) {
- r = rowSums(d)/(l - 2)
- i = 0
- j = 0
- tmp <- .C("out", as.double(d), as.double(r), as.integer(l), as.integer(i), as.integer(j))
- e2 = tmp[[5]]
- e1 = tmp[[4]]
- l1 = d[e1, e2]/2 + sum((d[e1,-c(e1,e2)] - d[e2,-c(e1,e2)])*w[-c(e1,e2)])/(2*(n-w[e1]-w[e2]))
- l2 = d[e1, e2]/2 + sum((d[e2,-c(e1,e2)] - d[e1,-c(e1,e2)])*w[-c(e1,e2)])/(2*(n-w[e1]-w[e2]))
- edge.length = c(l1, l2, edge.length)
- edge = rbind(c(k, nam[e2]), edge)
- edge = rbind(c(k, nam[e1]), edge)
- nam = c(nam[c(-e1, -e2)], k)
-
- dnew = (w[e1]*d[e1, ] + w[e2]*d[e2, ] - w[e1]*l1 - w[e2]*l2)/(w[e1] + w[e2])
- d = cbind(d, dnew)
- d = rbind(d, c(dnew, 0))
- d = d[-c(e1, e2), -c(e1, e2)]
- w = c(w, w[e1] + w[e2])
- w = w[-c(e1, e2)]
- k = k - 1
- l = l - 1
- }
- edge.length=c(d[2,1],edge.length)
- result = list(edge = rbind(c(nam[2], nam[1]), edge),
- edge.length=edge.length, tip.label = labels, Nnode=m)
- class(result) <- "phylo"
- reorder(result)
-}
-
-
-PNJ <- function (data)
-{
- q <- l <- r <- length(data)
- weight <- attr(data,"weight")
-
- height = NULL
- parentNodes <- NULL
- childNodes <- NULL
- nam <- names(data)
- tip.label <- nam
- edge = 1:q
-
- z = 0
- D = matrix(0, q, q)
-
- for (i in 1:(l - 1)) {
- for (j in (i + 1):l) {
- w = (data[[i]] * data[[j]]) %*% c(1, 1, 1, 1)
- D[i, j] = sum(weight[w==0])
- }
- }
-
- while (l > 1) {
- l = l - 1
- z = z + 1
- d = D + t(D)
- if(l>1) r = rowSums(d)/(l-1)
- if(l==1) r = rowSums(d)
- M = d - outer(r,r,"+")
- diag(M) = Inf
-
- e=which.min(M)
- e0=e%%length(r)
- e1 = ifelse(e0==0, length(r), e0)
- e2= ifelse(e0==0, e%/%length(r), e%/%length(r) + 1)
-
- ind = c(e1,e2)
- len = d[e]/2
- nam = c(nam[-ind], as.character(-l))
-
- parentNodes = c(parentNodes,-l,-l)
- childNodes = c(childNodes,edge[e1],edge[e2])
-
- height = c(height, len, len)
- edge = c(edge[-ind], -l)
- w = (data[[e1]] * data[[e2]]) %*% c(1, 1, 1, 1)
- w = which(w == 0)
- newDat = data[[e1]] * data[[e2]]
- newDat[w, ] = data[[e1]][w, ] + data[[e2]][w, ]
- data = data[-c(e1,e2)]
- data[[l]] = newDat
- if (l > 1) {
- D = as.matrix(D[, -ind])
- D = D[-ind, ]
- dv = numeric(l - 1)
- for (i in 1:(l - 1)) {
- w = (data[[i]] * data[[l]]) %*% c(1, 1, 1, 1)
- dv[i] = sum(weight[w==0])
- }
- D = cbind(D, dv)
- D = rbind(D, 0)
- }
- }
- tree <- list(edge = cbind(as.character(parentNodes),as.character(childNodes)),tip.label=tip.label)
- class(tree) <- "phylo"
- tree <- old2new.phylo(tree)
- reorder(tree)
-}
-
-
-
-#
# Maximum likelihood estimation
#
discrete.gamma <- function (alpha, k)
@@ -282,8 +97,8 @@ subsChoice <- function(type=c("JC", "F81", "K80", "HKY", "TrNe", "TrN", "TPM1",
modelTest <- function (object, tree = NULL, model = c("JC", "F81", "K80",
"HKY", "SYM", "GTR"), G = TRUE, I = TRUE, k = 4, control = pml.control(epsilon = 1e-08,
- maxit = 3, trace = 1), multicore = FALSE)
-{
+ maxit = 10, trace = 1), multicore = FALSE)
+{
if (class(object) == "phyDat")
data = object
if (class(object) == "pml") {
@@ -291,20 +106,23 @@ modelTest <- function (object, tree = NULL, model = c("JC", "F81", "K80",
if (is.null(tree))
tree = object$tree
}
-
+
if(attr(data, "type")=="DNA") type = 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(attr(data, "type")=="AA") type = .aamodels
-# type = c("WAG", "JTT", "Dayhoff", "LG", "cpREV", "mtmam", "mtArt", "MtZoa", "mtREV24")
+ "K81", "TPM1u", "TPM2", "TPM2u", "TPM3", "TPM3u", "TIM1e",
+ "TIM1", "TIM2e", "TIM2", "TIM3e", "TIM3", "TVMe", "TVM",
+ "SYM", "GTR")
+ if(attr(data, "type")=="AA") type = .aamodels
model = match.arg(model, type, TRUE)
-
+
env = new.env()
assign("data", data, envir=env)
-
+
if (is.null(tree))
tree = NJ(dist.hamming(data))
+ else{
+ tree <- nnls.phylo(tree, dist.ml(data))
+ # may need something faster for trees > 500 taxa
+ }
trace <- control$trace
control$trace = trace - 1
fit = pml(tree, data)
@@ -314,10 +132,10 @@ modelTest <- function (object, tree = NULL, model = c("JC", "F81", "K80",
nseq = sum(attr(data, "weight"))
fitPar = function(model, fit, G, I, k) {
m = 1
- res = matrix(NA, n, 5)
+ res = matrix(NA, n, 6)
res = as.data.frame(res)
- colnames(res) = c("Model", "df", "logLik", "AIC", "BIC")
- data.frame(c("Model", "df", "logLik", "AIC", "BIC"))
+ colnames(res) = c("Model", "df", "logLik", "AIC", "AICc", "BIC")
+ data.frame(c("Model", "df", "logLik", "AIC", "AICc", "BIC"))
calls = vector("list", n)
trees = vector("list", n)
fittmp = optim.pml(fit, model = model, control = control)
@@ -325,20 +143,22 @@ modelTest <- function (object, tree = NULL, model = c("JC", "F81", "K80",
res[m, 2] = fittmp$df
res[m, 3] = fittmp$logLik
res[m, 4] = AIC(fittmp)
- res[m, 5] = AIC(fittmp, k = log(nseq))
+ res[m, 5] = AICc(fittmp)
+ res[m, 6] = AIC(fittmp, k = log(nseq))
calls[[m]] = fittmp$call
-
+
trees[[m]] = fittmp$tree
m = m + 1
if (I) {
if(trace>0)print(paste(model, "+I", sep = ""))
fitI = optim.pml(fittmp, model = model, optInv = TRUE,
- control = control)
+ control = control)
res[m, 1] = paste(model, "+I", sep = "")
res[m, 2] = fitI$df
res[m, 3] = fitI$logLik
res[m, 4] = AIC(fitI)
- res[m, 5] = AIC(fitI, k = log(nseq))
+ res[m, 5] = AICc(fitI)
+ res[m, 6] = AIC(fitI, k = log(nseq))
calls[[m]] = fitI$call
trees[[m]] = fitI$tree
m = m + 1
@@ -347,12 +167,13 @@ modelTest <- function (object, tree = NULL, model = c("JC", "F81", "K80",
if(trace>0)print(paste(model, "+G", sep = ""))
fitG = update(fittmp, k = k)
fitG = optim.pml(fitG, model = model, optGamma = TRUE,
- control = control)
+ control = control)
res[m, 1] = paste(model, "+G", sep = "")
res[m, 2] = fitG$df
res[m, 3] = fitG$logLik
res[m, 4] = AIC(fitG)
- res[m, 5] = AIC(fitG, k = log(nseq))
+ res[m, 5] = AICc(fitG)
+ res[m, 6] = AIC(fitG, k = log(nseq))
calls[[m]] = fitG$call
trees[[m]] = fitG$tree
m = m + 1
@@ -360,12 +181,13 @@ modelTest <- function (object, tree = NULL, model = c("JC", "F81", "K80",
if (G & I) {
if(trace>0)print(paste(model, "+G+I", sep = ""))
fitGI = optim.pml(fitG, model = model, optGamma = TRUE,
- optInv = TRUE, control = control)
+ optInv = TRUE, control = control)
res[m, 1] = paste(model, "+G+I", sep = "")
res[m, 2] = fitGI$df
res[m, 3] = fitGI$logLik
res[m, 4] = AIC(fitGI)
- res[m, 5] = AIC(fitGI, k = log(nseq))
+ res[m, 5] = AICc(fitGI)
+ res[m, 6] = AIC(fitGI, k = log(nseq))
calls[[m]] = fitGI$call
trees[[m]] = fitGI$tree
m = m + 1
@@ -374,7 +196,7 @@ modelTest <- function (object, tree = NULL, model = c("JC", "F81", "K80",
}
eval.success <- FALSE
if (!eval.success & multicore) {
-# !require(parallel) ||
+ # !require(parallel) ||
if (.Platform$GUI != "X11") {
warning("package 'parallel' not found or GUI is used, \n analysis is performed in serial")
}
@@ -385,10 +207,10 @@ modelTest <- function (object, tree = NULL, model = c("JC", "F81", "K80",
}
if (!eval.success)
res <- RES <- lapply(model, fitPar, fit, G, I, k)
-
- RESULT = matrix(NA, n * l, 5)
+
+ RESULT = matrix(NA, n * l, 6)
RESULT = as.data.frame(RESULT)
- colnames(RESULT) = c("Model", "df", "logLik", "AIC", "BIC")
+ colnames(RESULT) = c("Model", "df", "logLik", "AIC", "AICc", "BIC")
for (i in 1:l) RESULT[((i - 1) * n + 1):(n * i), ] = RES[[i]][[1]]
for(i in 1:l){
for(j in 1:n){
@@ -398,7 +220,7 @@ modelTest <- function (object, tree = NULL, model = c("JC", "F81", "K80",
tmpmod["tree"] = call(tname)
if(!is.null(tmpmod[["k"]]))tmpmod["k"] = k
if(attr(data, "type")=="AA") tmpmod["model"] = RES[[i]][[1]][1,1]
- assign(tname, RES[[i]][[2]][[j]], envir=env)
+ assign(tname, RES[[i]][[2]][[j]], envir=env)
assign(mo, tmpmod, envir=env)
}
}
@@ -487,6 +309,20 @@ logLik.pml <- function(object,...){
res
}
+
+AICc <- function (object, ...)
+ UseMethod("AICc")
+
+
+AICc.pml <- function(object, ...){
+ n = sum(object$weight)
+ k = object$df
+# if(k>=(n-1))return(NULL)
+ res = AIC(object)
+ res + (2*k*(k+1))/(n-k-1)
+}
+
+
anova.pml <- function (object, ...)
{
X <- c(list(object), list(...))
@@ -568,13 +404,14 @@ getP <- function(el, eig=edQt(), g=1.0){
}
-lli = function (data, tree, ...)
+lli <- function (data, tree=NULL, ...)
{
- contrast = attr(data, "contrast")
- nr = attr(data, "nr")
- nc = attr(data, "nc")
- nco = as.integer(dim(contrast)[1])
- .Call("invSites", data[tree$tip.label], as.integer(nr), as.integer(nc), contrast, as.integer(nco))
+ contrast <- attr(data, "contrast")
+ nr <- attr(data, "nr")
+ nc <- attr(data, "nc")
+ nco <- as.integer(dim(contrast)[1])
+ if(!is.null(tree)) data <- subset(data, tree$tip.label)
+ .Call("invSites", data, as.integer(nr), as.integer(nc), contrast, as.integer(nco))
}
@@ -620,18 +457,26 @@ edQ2 <- function(Q){
}
-
-pml.free <- function(){.C("ll_free")}
+pml.free <- function(){
+ .C("ll_free")
+# rm(.INV, .iind, envir = parent.frame())
+}
pml.init <- function(data, k=1L){
- nTips <- length(data)
- nr <- attr(data, "nr")
- nc <- attr(data, "nc")
- .C("ll_init", as.integer(nr), as.integer(nTips), as.integer(nc), as.integer(k))
+ nTips <- length(data)
+ nr <- attr(data, "nr")
+ nc <- attr(data, "nc")
+ .C("ll_init", as.integer(nr), as.integer(nTips), as.integer(nc), as.integer(k))
+ INV <- lli(data) #, tree
+# .iind <<- which((INV %*% rep(1, nc)) > 0)
+# .INV <<- Matrix(INV, sparse=TRUE)
+ assign(".iind", which((INV %*% rep(1, nc)) > 0), envir=parent.frame())
+ assign(".INV", Matrix(INV, sparse=TRUE), envir=parent.frame())
}
+
pml.free2 <- function(){.C("ll_free2")}
pml.init2 <- function(data, k=1L){
@@ -643,36 +488,6 @@ pml.init2 <- function(data, k=1L){
}
-ll <- function(dat1, tree, bf = c(0.25, 0.25, 0.25, 0.25), g = 1,
- eig = NULL, assign.dat = FALSE, ...)
-{
- q = length(tree$tip.label)
- node <- tree$edge[, 1]
- edge <- tree$edge[, 2]
- m = length(edge) + 1
-# if (is.null(eig)) eig = edQt(bf = bf, Q = Q) # raus
- el <- tree$edge.length
- P <- getP(el, eig, g)
- nr <- as.integer(attr(dat1,"nr"))
- nc <- as.integer(attr(dat1,"nc"))
- node = as.integer(node-min(node))
- edge = as.integer(edge-1L)
- nTips = as.integer(length(tree$tip))
- mNodes = as.integer(max(node) + 1L)
- contrast = attr(dat1, "contrast")
- nco = as.integer(dim(contrast)[1])
- res <- .Call("LogLik2", dat1[tree$tip.label], P, nr, nc, node, edge, nTips, mNodes, contrast, nco)
- result = res[[1]] %*% bf # root statt 1
- if (assign.dat){
- dat = vector(mode = "list", length = m)
- dat[(q+1):m] <- res
- attr(dat, "names") = c(tree$tip.label, as.character((q + 1):m))
- assign("asdf", dat, envir = parent.frame(n = 1))
- }
- result
-}
-
-
fn.quartet <- function(old.el, eig, bf, dat, g=1, w=1, weight, ll.0) {
l= length(dat[,1])
ll = ll.0
@@ -774,7 +589,6 @@ pml.nni <- function (tree, data, w, g, eig, bf, ll.0, ll, ...)
INDEX <- indexNNI(tree)
rootEdges <- attr(INDEX,"root")
.dat <- NULL
-
data = getCols(data, tree$tip)
parent = tree$edge[,1]
@@ -955,40 +769,6 @@ optim.quartet <- function (old.el, eig, bf, dat, g = 1, w = 1, weight, ll.0 = we
}
-optim.quartet2 <- function (old.el, eig, bf, dat1, dat2, dat3, dat4, g = 1, w = 1, weight, ll.0 = weight *
- 0, control = list(eps = 1e-08, maxit = 5, trace = 0), llcomp=-Inf, evi, contrast, contrast2, ext=c(FALSE, FALSE, FALSE, FALSE))
-{
- eps = 1
- iter = 0
- while (eps > control$eps && iter < control$maxit) {
- tmp <- fn.quartet2(old.el = old.el, eig = eig, bf = bf, dat1 = dat1, dat2 = dat2, dat3 = dat3, dat4 = dat4,
- g = g, w = w, weight = weight, ll.0 = ll.0, contrast=contrast, ext=ext)
- old.ll = tmp$ll
- el1 <- fs3(old.el[1], eig, tmp$res[, 1], dat1, weight,
- g = g, w = w, bf = bf, ll.0 = ll.0, contrast=contrast, contrast2=contrast2, evi=evi, ext = ext[1], getA=TRUE, getB=FALSE)
- el2 <- fs3(old.el[2], eig, el1[[2]], dat2, weight,
- g = g, w = w, bf = bf, ll.0 = ll.0, contrast=contrast, contrast2=contrast2, evi=evi, ext = ext[2], getA=TRUE, getB=FALSE)
- el5 <- fs3(old.el[5], eig, el2[[2]], tmp$res[, 2], weight,
- g = g, w = w, bf = bf, ll.0 = ll.0, contrast=contrast, contrast2=contrast2, evi=evi, ext = 0L, getA=FALSE, getB=TRUE)
- el3 <- fs3(old.el[3], eig, el5[[3]], dat3, weight,
- g = g, w = w, bf = bf, ll.0 = ll.0, contrast=contrast, contrast2=contrast2, evi=evi, ext = ext[3], getA=TRUE, getB=FALSE)
- el4 <- fs3(old.el[4], eig, el3[[2]], dat4, weight,
- g = g, w = w, bf = bf, ll.0 = ll.0, contrast=contrast, contrast2=contrast2, evi=evi, ext = ext[4], getA=FALSE, getB=FALSE)
- old.el[1] = el1[[1]]
- old.el[2] = el2[[1]]
- old.el[3] = el3[[1]]
- old.el[4] = el4[[1]]
- old.el[5] = el5[[1]]
- iter = iter + 1
- ll = el4[[4]]
- eps = (old.ll - ll) / ll
- if(ll<llcomp)return(list(old.el, ll))
- old.ll = ll
- }
- list(old.el, ll)
-}
-
-
plot.pml<-function(x,...)plot.phylo(x$tree,...)
@@ -1110,7 +890,7 @@ optim.pml <- function (object, optNni = FALSE, optBf = FALSE, optQ = FALSE,
k = object$k
if(k==1 & optGamma){
optGamma = FALSE
- warning('only one rate class, ignored optGamma')
+ message('only one rate class, ignored optGamma')
}
shape = object$shape
w = object$w
@@ -1137,9 +917,15 @@ optim.pml <- function (object, optNni = FALSE, optBf = FALSE, optQ = FALSE,
nc <- as.integer(attr(data, "nc"))
nTips <- as.integer(length(tree$tip.label))
- on.exit(.C("ll_free"))
- .C("ll_init", nr, nTips, nc, as.integer(k))
-
+# on.exit(.C("ll_free"))
+# .C("ll_init", nr, nTips, nc, as.integer(k))
+ .INV <- .iind <- NULL
+ on.exit({
+ pml.free()
+# rm(.INV, .iind)
+ })
+ pml.init(data, k)
+
if (optEdge) {
res <- optimEdge(tree, data, eig=eig, w=w, g=g, bf=bf, rate=rate, ll.0=ll.0, INV=INV,
control = pml.control(epsilon = 1e-07, maxit = 5, trace=trace - 1))
@@ -1295,6 +1081,7 @@ optim.pml <- function (object, optNni = FALSE, optBf = FALSE, optQ = FALSE,
else{
tmp <- rooted.nni(tree, data, eig=eig, w=w, g=g, bf=bf, rate=rate, ll.0=ll.0, INV=INV, ...)
swap = swap + tmp$swap
+ res <- optimRooted(tmp$tree, data, eig=eig, w=w, g=g, bf=bf, rate=rate, ll.0=ll.0, INV=INV, control = pml.control(epsilon = 1e-07, maxit = 5, trace = trace-1))
tree <- tmp$tree
ll2 = tmp$logLik
}
@@ -1417,7 +1204,7 @@ optimEdge <- function (tree, data, eig=eig, w=w, g=g, bf=bf, rate=rate, ll.0=ll.
tree <- reorder(tree, "postorder")
nTips <- length(tree$tip)
el <- tree$edge.length
- tree$edge.length[el < 0] <- 1e-08
+ tree$edge.length[el < 1e-08] <- 1e-08
oldtree = tree
k = length(w)
data = subset(data, tree$tip)
@@ -1432,21 +1219,17 @@ optimEdge <- function (tree, data, eig=eig, w=w, g=g, bf=bf, rate=rate, ll.0=ll.
treeP = tree
tree = reorder(tree)
- el <- tree$edge.length
child = tree$edge[, 2]
parent = tree$edge[, 1]
- loli = parent[1]
-
- pvec <- integer(max(tree$edge))
+ m <- max(tree$edge)
+ pvec <- integer(m)
pvec[child] <- parent
- EL = numeric(max(child))
+ EL = numeric(m)
EL[child] = tree$edge.length
- nTips = min(parent) - 1
n = length(tree$edge.length)
- lt = length(tree$tip)
nr = as.integer(length(weight))
nc = as.integer(length(bf))
@@ -1455,63 +1238,37 @@ optimEdge <- function (tree, data, eig=eig, w=w, g=g, bf=bf, rate=rate, ll.0=ll.
lg = k
rootNode = getRoot(tree)
ScaleEPS = 1.0/4294967296.0
- anc = Ancestors(tree, 1:max(tree$edge), "parent")
+ anc = Ancestors(tree, 1:m, "parent")
+ anc0 = as.integer(c(0L, anc))
while (eps > control$eps && iter < control$maxit) {
blub3 <- .Call("extractScale", as.integer(rootNode), w, g, as.integer(nr), as.integer(nc), as.integer(nTips))
rowM = apply(blub3, 1, min)
blub3 = (blub3-rowM)
blub3 = ScaleEPS ^ (blub3)
-
- for(j in 1:n) {
- ch = child[j]
- pa = parent[j]
-
- while(loli != pa){
- blub <- .Call("moveloli", as.integer(loli), as.integer(anc[loli]), eig, EL[loli], w, g, as.integer(nr), as.integer(nc), as.integer(nTips))
- loli=anc[loli]
- }
-
- old.el = tree$edge.length[j]
- if (old.el < 1e-8) old.el <- 1e-8
-
- X <- .Call("moveDad", data, as.integer(pa), as.integer(ch), eig, evi, old.el, w, g, as.integer(nr), as.integer(nc), as.integer(nTips), as.double(contrast), as.double(contrast2), nco)
- # in moveDad, das sollte teuer sein in Kopien
- for(i in 1:k)X[[i]] <- X[[i]] * blub3[,i]
-
- newEL <- .Call("FS5", eig, nc, as.double(old.el), as.double(w), as.double(g), X, as.integer(length(w)), as.integer(length(weight)), as.double(bf), as.double(weight), as.double(ll.0))
-
- el[j] = newEL[[1]] # nur EL?
- EL[ch] = newEL[[1]]
- if (child[j] > nTips) loli = child[j]
- else loli = parent[j]
-
- blub <- .Call("updateLL", data, as.integer(pa), as.integer(ch), eig, newEL[[1]], w, g,
- as.integer(nr), as.integer(nc), as.integer(nTips), as.double(contrast), nco)
-
- tree$edge.length = el
- }
- tree$edge.length = el
+ EL <- .Call("optE", as.integer(parent), as.integer(child),
+ as.integer(anc0), eig, evi, EL, w, g, as.integer(nr), as.integer(nc),
+ as.integer(nTips), as.double(contrast),
+ as.double(contrast2), nco, blub3, 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.fit2(treeP, data, bf=bf, g=g, w=w, eig=eig, ll.0=ll.0, k=k)
eps = ( old.ll - newll ) / newll
if( eps <0 ) return(list(oldtree, old.ll))
- oldtree = tree
+ oldtree = treeP
if(control$trace>1) cat(old.ll, " -> ", newll, "\n")
old.ll = newll
- loli = parent[1]
-
+ # loli = parent[1]
}
if(control$trace>0) cat(start.ll, " -> ", newll, "\n")
list(tree=treeP, logLik=newll, c(eps, iter))
}
-
# bf raus C naeher
+# data=data, k=k, g=g, w=w, eig=eig, bf=bf, ll.0=ll.0, INV=INV)
pml.move <- function(EDGE, el, data, g, w, eig, k, nTips, bf){
node <- EDGE[, 1]
edge <- EDGE[, 2]
@@ -1610,7 +1367,7 @@ optimPartInv <- function (object, inv = 0.01, ...)
}
res = optimize(f = fn, interval = c(0, 1), lower = 0, upper = 1,
maximum = TRUE, tol = 1e-04, object, ...)
- print(res[[2]])
+# print(res[[2]])
res[[1]]
}
@@ -1791,6 +1548,11 @@ pmlPart2multiPhylo <- function(x){
}
+plot.pmlPart<- function(x, ...){
+ plot(pmlPart2multiPhylo(x), ...)
+}
+
+
pmlPart <- function (formula, object, control=pml.control(epsilon=1e-8, maxit=10, trace=1), model=NULL, rooted=FALSE, ...)
{
call <- match.call()
@@ -1912,9 +1674,6 @@ pmlPart <- function (formula, object, control=pml.control(epsilon=1e-8, maxit=10
}
-#
-# Distance Matrix methods
-#
bip <- function (obj)
{
@@ -1928,286 +1687,6 @@ bip <- function (obj)
}
-# as.Matrix, sparse = TRUE,
-designTree <- function(tree, method="unrooted", sparse=FALSE, ...){
- if (!is.na(pmatch(method, "all")))
- method <- "unrooted"
- METHOD <- c("unrooted", "rooted")
- method <- pmatch(method, METHOD)
- if (is.na(method)) stop("invalid method")
- if (method == -1) stop("ambiguous method")
- if(!is.rooted(tree) & method==2) stop("tree has to be rooted")
- if(method==1){ X <- designUnrooted(tree,...)
- if(sparse) X = Matrix(X)
- }
- if(method==2) X <- designUltra(tree, sparse=sparse,...)
- X
-}
-
-
-# splits now work
-designUnrooted = function (tree, order = NULL)
-{
- if(inherits(tree, "phylo")){
- if (is.rooted(tree))
- tree = unroot(tree)
- p = bipartition(tree)
- }
- if(inherits(tree, "splits")) p <- as.matrix(tree)
- if (!is.null(order))
- p = p[, order]
-
- m = dim(p)[2]
- ind = rowSums(p)
- p=p[ind!=m,]
- n = dim(p)[1]
- res = matrix(0, (m - 1) * m/2, n)
- k = 1
- for (i in 1:(m - 1)) {
- for (j in (i + 1):m) {
- res[k, ] = p[, i] != p[, j]
- k = k + 1
- }
- }
- if(inherits(tree, "phylo"))colnames(res) = paste(tree$edge[, 1], tree$edge[, 2], sep = "<->")
- res
-}
-
-
-designUltra <- function (tree, sparse=FALSE)
-{
- if (is.null(attr(tree, "order")) || attr(tree, "order") == "cladewise")
- tree = reorder(tree, "postorder")
- leri = allChildren(tree)
- bp = bip(tree)
- n = length(tree$tip)
- l = tree$Nnode
- nodes = integer(l)
- k = 1L
- u=numeric( n * (n - 1)/2)
- v=numeric( n * (n - 1)/2)
- m = 1L
- for (i in 1:length(leri)) {
- if (!is.null(leri[[i]])) {
- ind = getIndex(bp[[leri[[i]][1] ]], bp[[leri[[i]][2] ]], n)
- li = length(ind)
- v[m: (m+li-1)]=k
- u[m: (m+li-1)]=ind
- nodes[k]=i
- m = m+li
- k = k + 1L
- }
- }
- if(sparse) X = sparseMatrix(i=u,j=v, x=2L)
- else{
- X = matrix(0L, n * (n - 1)/2, l)
- X[cbind(u,v)]=2L
- }
- colnames(X) = nodes
- attr(X, "nodes") = nodes
- X
-}
-
-
-nnls.tree <- function(dm, tree, rooted=FALSE, trace=1){
- if(is.rooted(tree) & rooted==FALSE){
- tree = unroot(tree)
- warning("tree was rooted, I unrooted the tree!")
- }
- tree = reorder(tree, "postorder")
- dm = as.matrix(dm)
- k = dim(dm)[1]
- labels = tree$tip
- dm = dm[labels,labels]
- y = dm[lower.tri(dm)]
-#computing the design matrix from the tree
- if(rooted) X = designUltra(tree)
- else X = designUnrooted(tree)
-# na.action
- if(any(is.na(y))){
- ind = which(is.na(y))
- X = X[-ind,,drop=FALSE]
- y= y[-ind]
- }
-# LS solution
- fit = lm.fit(X,y)
- betahat = as.vector(fit$coefficients) # added as.vector
- if(rooted){
- bhat = numeric(max(tree$edge))
- bhat[as.integer(colnames(X))] = betahat
- betahat = bhat[tree$edge[,1]] - bhat[tree$edge[,2]]
- }
- if(!any(betahat<0)){
- RSS = sum(fit$residuals^2)
- if(trace)print(paste("RSS:", RSS))
- attr(tree, "RSS") = RSS
- tree$edge.length = betahat # deleted []
- return(tree)
- }
-# non-negative LS
- n = dim(X)[2]
- Dmat <- crossprod(X) # cross-product computations
- dvec <- crossprod(X, y)
- if(rooted){
- l = nrow(tree$edge)
- Amat = matrix(0, n, l)
- ind = match(tree$edge[,1], colnames(X))
- Amat[cbind(ind, 1:l)] = 1
- ind = match(tree$edge[,2], colnames(X))
- Amat[cbind(ind, 1:l)] = -1
- }
- else Amat <- diag(n)
- betahat <- quadprog::solve.QP(Dmat,dvec,Amat)$sol # quadratic programing solving
- RSS = sum((y-(X%*%betahat))^2)
-
- if(rooted){
- bhat = numeric(max(tree$edge))
- bhat[as.integer(colnames(X))] = betahat
- betahat = bhat[tree$edge[,1]] - bhat[tree$edge[,2]]
- }
- tree$edge.length = betahat # deleted []
- if(trace)print(paste("RSS:", RSS))
- attr(tree, "RSS") = RSS
- tree
-}
-
-
-nnls.phylo <- function(x, dm, rooted=FALSE, trace=0){
- nnls.tree(dm, x, rooted, trace=trace)
-}
-
-
-nnls.splits <- 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)
- 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)
- 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]
- }
-
- Dmat <- crossprod(X) # cross-product computations
- dvec <- crossprod(X, y)
- betahat <- as.vector(solve(Dmat, dvec))
-
- if(!any(betahat<0)){
- RSS = sum((y-(X%*%betahat))^2)
- if(trace)print(paste("RSS:", RSS))
- attr(x, "RSS") = RSS
- attr(x, "weights") = betahat
- return(x)
- }
- n = dim(X)[2]
-
- 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]
- attr(x, "weights") = betahat[ind]
- if(trace)print(paste("RSS:", RSS))
- attr(x, "RSS") = RSS
- x
-}
-
-nnls.networx <- function(x, dm){
- spl <- attr(x, "splits")
- spl2 <- nnls.splits(spl, dm)
- weight <- attr(spl, "weight")
- weight[] <- 0
- weight[match(spl2, spl)] = attr(spl2, "weight")
- attr(attr(x, "splits"), "weight") <- weight
- x$edge.length = weight[x$splitIndex]
- x
-}
-
-
-designSplits <- function (x, splits = "all", ...)
-{
- if (!is.na(pmatch(splits, "all")))
- splits <- "all"
- if(inherits(x, "splits")) return(designUnrooted(x))
- SPLITS <- c("all", "star") #,"caterpillar")
- splits <- pmatch(splits, SPLITS)
- 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)
- return(X)
-}
-
-# add return splits=FALSE
-designAll <- function(n, add.split=FALSE){
- Y = matrix(0L, n*(n-1)/2, n)
- k = 1
- for(i in 1:(n-1)){
- for(j in (i+1):n){
- Y[k,c(i,j)]=1L
- k=k+1L
- }
- }
- m <- n-1L
- X <- matrix(0L, m+1, 2^m)
- for(i in 1:m)
- X[i, ] <- rep(rep(c(0L,1L), each=2^(i-1)),2^(m-i))
- X <- X[,-1]
- if(!add.split) return((Y%*%X)%%2)
- list(X=(Y%*%X)%%2,Splits=t(X))
-}
-
-
-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)))
- res
-}
-
-
bipart <- function(obj){
if (is.null(attr(obj, "order")) || attr(obj, "order") == "cladewise")
obj <- reorder(obj, "postorder")
@@ -2788,7 +2267,7 @@ optimMixBf <- function(object, bf=c(.25,.25,.25,.25), omega,...){
}
res = optim(par=lbf, fn=fn, gr=NULL, method="Nelder-Mead",
control=list(fnscale=-1, maxit=500), object, omega=omega,...)
- print(res[[2]])
+# print(res[[2]])
bf = exp(c(res[[1]],0))
bf = bf/sum(bf)
}
@@ -2806,7 +2285,7 @@ optimMixInv <- function(object, inv=0.01, omega,...){
}
res = optimize(f=fn, interval = c(0,1), lower = 0, upper = 1, maximum = TRUE,
tol = .0001, object, omega=omega,...)
- print(res[[2]])
+# print(res[[2]])
res[[1]]
}
@@ -3626,7 +3105,7 @@ checkLabels <- function(tree, tip){
plotBS <- function (tree, BStrees, type = "unrooted", bs.col = "black",
- bs.adj = NULL, ...)
+ bs.adj = NULL, p=80, ...)
{
# prop.clades raus??
prop.clades <- function(phy, ..., part = NULL, rooted = FALSE) {
@@ -3658,7 +3137,7 @@ plotBS <- function (tree, BStrees, type = "unrooted", bs.col = "black",
type <- match.arg(type, c("phylogram", "cladogram", "fan",
"unrooted", "radial"))
if (type == "phylogram" | type == "cladogram") {
- if (!is.rooted(tree))
+ if (!is.rooted(tree) & !is.null(tree$edge.length))
tree2 = midpoint(tree)
else tree2=tree
plot(tree2, type = type, ...)
@@ -3668,7 +3147,7 @@ plotBS <- function (tree, BStrees, type = "unrooted", bs.col = "black",
x = prop.clades(tree, BStrees)
x = round((x/length(BStrees)) * 100)
tree$node.label = x
- label = c(rep("", length(tree$tip)), x)
+ label = c(rep(0, length(tree$tip)), x)
ind <- get("last_plot.phylo", envir = .PlotPhyloEnv)$edge[,
2]
if (type == "phylogram" | type == "cladogram") {
@@ -3677,16 +3156,18 @@ plotBS <- function (tree, BStrees, type = "unrooted", bs.col = "black",
label[root] = 0
ind2 = matchEdges(tree2, tree)
label = label[ind2]
- ind = which(label > 0)
+ ind = which(label > p)
+# browser()
if (is.null(bs.adj))
bs.adj = c(1, 1)
- nodelabels(text = label[ind], node = ind, frame = "none",
+ if(length(ind)>0)nodelabels(text = label[ind], node = ind, frame = "none",
col = bs.col, adj = bs.adj, ...)
}
else {
if (is.null(bs.adj))
bs.adj = c(0.5, 0.5)
- edgelabels(label[ind], frame = "none", col = bs.col,
+ ind2 = which(label[ind]>p)
+ if(length(ind2>0))edgelabels(label[ind][ind2],ind2, frame = "none", col = bs.col,
adj = bs.adj, ...)
}
invisible(tree)
@@ -3773,6 +3254,80 @@ pml.fit2 <- function (tree, data, bf = rep(1/length(levels), length(levels)),
}
+pml.fit4 <- function (tree, data, bf = rep(1/length(levels), length(levels)),
+ shape = 1, k = 1, Q = rep(1, length(levels) * (length(levels) - 1)/2),
+ levels = attr(data, "levels"), inv = 0, rate = 1, g = NULL, w = NULL,
+ eig = NULL, INV = NULL, ll.0 = NULL, llMix = NULL, wMix = 0, ..., site=FALSE)
+{
+ weight <- as.double(attr(data, "weight"))
+ nr <- as.integer(attr(data, "nr"))
+ nc <- as.integer(attr(data, "nc"))
+ nTips <- as.integer(length(tree$tip.label))
+ k <- as.integer(k)
+ m = 1
+ if (is.null(eig))
+ eig = edQt(bf = bf, Q = Q)
+ if (is.null(w)) {
+ w = rep(1/k, k)
+ if (inv > 0)
+ w <- (1 - inv) * w
+ if (wMix > 0)
+ w <- (1 - wMix) * w
+ }
+ if (is.null(g)) {
+ g = discrete.gamma(shape, k)
+ if (inv > 0)
+ g <- g/(1 - inv)
+ g <- g * rate
+ }
+ # inv0 <- inv
+ if(any(g<.gEps)){
+ for(i in 1:length(g)){
+ if(g[i]<.gEps){
+ inv <- inv+w[i]
+ }
+ }
+ w <- w[g>.gEps]
+ g <- g[g>.gEps]
+ # kmax <- k
+ k <- length(w)
+ }
+# .iind <- get(".iind", parent.frame())
+# .INV <- get(".INV", parent.frame())
+# if(is.null(ll.0))
+
+ if (is.null(ll.0)){
+ ll.0 <- numeric(attr(data,"nr"))
+ }
+ if(inv>0)
+ ll.0 <- as.matrix(INV %*% (bf * inv))
+# if(inv>0)
+# ll.0 <- as.matrix(.INV %*% (bf * inv))
+
+ node <- tree$edge[, 1]
+ edge <- tree$edge[, 2]
+ root <- as.integer(node[length(node)])
+ el <- as.double(tree$edge.length)
+ node = as.integer(node - nTips - 1L) # min(node))
+ edge = as.integer(edge - 1L)
+
+ contrast = attr(data, "contrast")
+ nco = as.integer(dim(contrast)[1])
+
+ siteLik <- .Call("PML4", dlist=data, el, as.double(w), as.double(g), nr, nc, k, eig, as.double(bf),
+ node, edge, nTips, root, nco, contrast, N=as.integer(length(edge)))
+# if(inv>0) siteLik[.iind] = log(exp(siteLik[.iind]) + ll.0[.iind])
+ ind = which(ll.0>0)
+# if(!is.null(ll.0)) siteLik[.iind] = log(exp(siteLik[.iind]) + ll.0[.iind])
+ if(!is.null(ll.0)) siteLik[ind] = log(exp(siteLik[ind]) + ll.0[ind])
+ if(wMix >0) siteLik <- log(exp(siteLik) * (1-wMix) + llMix )
+ loglik <- sum(weight * siteLik)
+ if(!site) return(loglik)
+ return(list(loglik=loglik, siteLik=siteLik)) #, resll=resll
+}
+
+
+
pml.fit <- function (tree, data, bf = rep(1/length(levels), length(levels)),
shape = 1, k = 1, Q = rep(1, length(levels) * (length(levels) - 1)/2),
levels = attr(data, "levels"), inv = 0, rate = 1, g = NULL, w = NULL,
@@ -3869,7 +3424,7 @@ pml <- function (tree, data, bf = NULL, Q = NULL, inv = 0, k = 1, shape = 1,
tree <- reorder(tree, "postorder")
if(any(tree$edge.length < 0)) {
tree$edge.length[tree$edge.length < 0] <- 1e-08
- warning("negative edges length changed to 0!")
+ message("negative edges length changed to 0!")
}
if (class(data)[1] != "phyDat") stop("data must be of class phyDat")
if(is.null(tree$edge.length)) stop("tree must have edge weights")
@@ -3947,47 +3502,43 @@ pml <- function (tree, data, bf = NULL, Q = NULL, inv = 0, k = 1, shape = 1,
optimRooted <- function(tree, data, eig=eig, w=w, g=g, bf=bf, rate=rate, ll.0=ll.0, INV=INV,
control = pml.control(epsilon = 1e-08, maxit = 25, trace=0), ...){
- ind0 = which(ll.0>0)
- contrast = attr(data, "contrast")
- tree$edge.length[tree$edge.length < 1e-08] <- 1e-08
+
+ tree$edge.length[tree$edge.length < 1e-08] <- 1e-08 # nicht richtig
nTips = as.integer(length(tree$tip.label))
k = length(w)
- weight = attr(data , "weight")
- nr = as.integer(attr(data , "nr"))
- nc = as.integer(attr(data , "nc"))
- # optimising rooted triplets
- optRoot0 <- function(t, tree, data, g, w, eig, bf, ll.0, k, INV){
+ # optimising rooted triplets
+ optRoot0 <- function(t, tree, data, g, w, eig, bf, ll.0, k){
l = length(tree$edge.length)
tree$edge.length[1:(l-1)] = tree$edge.length[1:(l-1)] + t
tree$edge.length[l] = tree$edge.length[l] - t
- loglik = pml.fit(tree, data, bf=bf, g=g, w=w, eig=eig, ll.0=ll.0, k=k, INV=INV)
+ loglik = pml.fit4(tree, data, bf=bf, g=g, w=w, eig=eig, INV=INV, ll.0=ll.0, k=k) #
loglik
}
# optim edges leading to the root
- optRoot2 <- function(t, tree, data, g, w, eig, bf, ll.0, k, INV){
+ optRoot2 <- function(t, tree, data, g, w, eig, bf, ll.0, k){
tree$edge.length = tree$edge.length + t #c(el1+t, el2-t)
- loglik = pml.fit(tree, data, bf=bf, g=g, w=w, eig=eig, ll.0=ll.0, k=k, INV=INV)
+ loglik = pml.fit4(tree, data, bf=bf, g=g, w=w, eig=eig, INV=INV, ll.0=ll.0, k=k) #, INV=INV
loglik
}
-
+ # scale whole tree
scaleEdges = function(t=1, trace=0, tree, data, ...){
fn = function(t, tree, data,...){
tree$edge.length = tree$edge.length*t
- pml.fit(tree, data, ...)
+ pml.fit4(tree, data, ...)
}
optimize(f=fn, interval=c(0.25,4), tree=tree, data=data, ..., maximum = TRUE,
tol = .00001)
}
parent = tree$edge[, 1]
child = tree$edge[, 2]
-
+
anc <- Ancestors(tree, 1:max(tree$edge), "parent")
sibs <- Siblings(tree, 1:max(tree$edge))
allKids <- cvector <- allChildren(tree)
rootNode = getRoot(tree)
- child2 = orderNNI(cvector, rootNode, nTips, TRUE)
+ child2 = orderNNI(tree, nTips) #(cvector, rootNode, nTips, TRUE)
lengthList <- edgeList <- vector("list", max(tree$edge))
for(i in tree$edge[,2]){
@@ -4003,28 +3554,35 @@ optimRooted <- function(tree, data, eig=eig, w=w, g=g, bf=bf, rate=rate, ll.0=ll
}
}
- ll <- pml.fit(tree, data, bf=bf, k=k, eig=eig, INV=INV, ll.0=ll.0, w=w, g=g)
- if(control$trace>2)cat("ll", ll, "\n")
+ treeList <- vector("list", max(child2))
+ for(i in child2){
+ pa <- anc[i]
+ kids <- allKids[[i]]
+ treeList[[i]] <- cbind(i, c(kids, pa))
+ }
+
+ ll <- pml.fit4(tree, data, bf=bf, k=k, eig=eig, ll.0=ll.0, INV=INV, w=w, g=g) #, INV=INV
+# if(control$trace>2)cat("ll", ll, "\n")
eps=10
iter = 1
-
+
EL = numeric(max(tree$edge))
EL[tree$edge[, 2]] = tree$edge.length
-
+
eps0 =1e-8
- tmp <- scaleEdges(t, trace=0, tree, data, bf = bf, k=k, ll.0=ll.0, INV=INV, eig = eig, w=w, g=g)
- if(control$trace>2)cat("scale", tmp[[2]], "\n")
+ tmp <- scaleEdges(t, trace=0, tree, data, bf = bf, k=k, ll.0=ll.0, eig = eig, w=w, g=g)
+# if(control$trace>2)cat("scale", tmp[[2]], "\n")
t = tmp[[1]]
tree$edge.length = tree$edge.length*t
el = tree$edge.length
EL[tree$edge[, 2]] = tree$edge.length
- ll2 <- pml.fit(tree, data, bf=bf, k=k, eig=eig, ll.0=ll.0, INV=INV, w=w, g=g)
+ ll2 <- pml.fit4(tree, data, bf=bf, k=k, eig=eig, INV=INV, ll.0=ll.0, w=w, g=g)
tmptree = tree
while(eps>control$eps && iter < control$maxit){
- ll2 <- pml.fit(tree, data, bf=bf, k=k, eig=eig, ll.0=ll.0, INV=INV, w=w, g=g)
+ ll2 <- pml.fit4(tree, data, bf=bf, k=k, eig=eig, INV=INV, ll.0=ll.0, w=w, g=g)
loli <- rootNode
children <- allKids[[rootNode]]
@@ -4032,64 +3590,59 @@ optimRooted <- function(tree, data, eig=eig, w=w, g=g, bf=bf, rate=rate, ll.0=ll
minEl = min(kidsEl)
kidsEl = kidsEl - minEl
-# EDGE = cbind(rootNode, children)
tmptree$edge = cbind(rootNode, children)
tmptree$edge.length = kidsEl
- t <- optimize(f=optRoot2,interval=c(1e-8,3), tmptree, data=data, k=k, g=g, w=w, eig=eig, bf=bf, ll.0=ll.0, INV=INV, maximum=TRUE)
- optRoot2(t[[1]], tmptree, data=data, k=k, g=g, w=w, eig=eig, bf=bf, ll.0=ll.0, INV=INV)
-
- if(control$trace>2)cat("optRoot", t[[2]], "\n")
+ t <- optimize(f=optRoot2,interval=c(1e-8,3), tmptree, data=data, k=k, g=g, w=w, eig=eig, bf=bf, ll.0=ll.0, maximum=TRUE)
+ optRoot2(t[[1]], tmptree, data=data, k=k, g=g, w=w, eig=eig, bf=bf, ll.0=ll.0)
+# if(control$trace>2)cat("optRoot", t[[2]], "\n")
ll3 = t[[2]]
EL[children] = kidsEl + t[[1]]
-
+
tree$edge.length = EL[tree$edge[, 2]]
- ll2 <- pml.fit(tree, data, bf=bf, k=k, eig=eig, ll.0=ll.0, INV=INV, w=w, g=g)
-
- for(i in 1:length(child2)){ # length(child2)
+ ll2 <- pml.fit4(tree, data, bf=bf, k=k, eig=eig, INV=INV, ll.0=ll.0, w=w, g=g)
+ for(i in 1:length(child2)){
dad = child2[i]
- if(dad>nTips ){ # kann raus
- pa = anc[dad]
- while(loli != pa){
- tmpKids= cvector[[loli]]
- tmpEdge = cbind(loli, tmpKids)
- pml.move(tmpEdge, EL[tmpKids], data, g, w, eig, k, nTips, bf)
- loli=anc[loli]
- }
- pml.move(edgeList[[dad]], EL[lengthList[[dad]]], data, g, w, eig, k, nTips, bf)
- children <- allKids[[dad]]
- kidsEl <- EL[children]
- minEl = min(kidsEl)
- kidsEl = kidsEl - minEl
- maxEl = minEl + EL[dad]
- kids=c(children, pa)
- EDGE = cbind(dad, kids)
- tmptree$edge = EDGE
- tmptree$edge.length = c(kidsEl, maxEl)
-
- tmptree2 = tmptree
- tmptree2$edge.length = c(EL[children], EL[dad])
-
- t0 = optRoot0(0, tmptree2, data, g, w, eig, bf, ll.0, k, INV)
-
- t = optimize(f=optRoot0, interval=c(0+eps0,maxEl-eps0), tmptree, data=data, g=g, w=w, eig=eig, bf=bf, ll.0=ll.0, k=k, INV=INV, maximum=TRUE)
-
- if(control$trace>2) cat("edge", t[[2]], "\n")
- if(!is.nan(t[[2]]) & t[[2]] > ll3){
- optRoot0(t[[1]], tmptree, data=data, g=g, w=w, eig=eig, bf=bf, ll.0=ll.0, k=k, INV=INV)
- EL[children] = kidsEl+t[[1]]
- EL[dad] = maxEl-t[[1]]
- ll3 = t[[2]]
- }
- else optRoot0(0, tmptree2, data, g, w, eig, bf, ll.0, k, INV)
- loli = dad
- tree$edge.length = EL[tree$edge[, 2]]
+
+ # if(dad>nTips ){ # kann raus
+ pa = anc[dad]
+ while(loli != pa){
+ tmpKids= cvector[[loli]]
+ tmpEdge = cbind(loli, tmpKids)
+ pml.move(tmpEdge, EL[tmpKids], data, g, w, eig, k, nTips, bf)
+ loli=anc[loli]
+ }
+ pml.move(edgeList[[dad]], EL[lengthList[[dad]]], data, g, w, eig, k, nTips, bf)
+
+ children <- allKids[[dad]]
+ kidsEl <- EL[children]
+ minEl = min(kidsEl)
+ maxEl = EL[dad]
+ EDGE = treeList[[dad]]
+
+ tmptree$edge = EDGE
+ tmptree$edge.length = c(kidsEl, maxEl)
+
+ t0 = optRoot0(0, tmptree, data, g, w, eig, bf, ll.0, k)
+
+ t = optimize(f=optRoot0, interval=c(-minEl+eps0,maxEl-eps0), tmptree, data=data, g=g, w=w, eig=eig, bf=bf, ll.0=ll.0, k=k, maximum=TRUE)
+
+# if(control$trace>2) cat("edge", t[[2]], "\n")
+ if(!is.nan(t[[2]]) & t[[2]] > ll3){
+ optRoot0(t[[1]], tmptree, data=data, g=g, w=w, eig=eig, bf=bf, ll.0=ll.0, k=k)
+ EL[children] = kidsEl+t[[1]]
+ EL[dad] = maxEl-t[[1]]
+ ll3 = t[[2]]
}
- }
+ else optRoot0(0, tmptree, data, g, w, eig, bf, ll.0, k)
+ loli = dad
+ # }
+ }
+ tree$edge.length = EL[tree$edge[, 2]]
ll2 <- pml.fit(tree, data, bf=bf, k=k, eig=eig, ll.0=ll.0, INV=INV, w=w, g=g)
eps = (ll - ll2) / ll2
- if(control$trace>1) cat(ll, " -> ", ll2, "\n")
+ if(control$trace>1) cat("optimRooted: ", ll, " -> ", ll2, "\n")
ll=ll2
iter = iter+1
}
@@ -4097,6 +3650,7 @@ optimRooted <- function(tree, data, eig=eig, w=w, g=g, bf=bf, rate=rate, ll.0=ll
}
+
# copy node likelihoods from C to R
getNodeLogLik = function(data, i, j=1L){
nr = attr(data, "nr")
@@ -4131,39 +3685,13 @@ index.nni <- function (ch, cvector, pvector, root)
}
-# nicer traversal
-orderNNI <- function (cvector, root, nTips, nni=TRUE)
-{
- #if(nni) l = sum(sapply(cvector, function(x, nTips) sum(x > nTips), nTips))
- # eleganter
- if(nni) l = sum(unlist(cvector) > nTips)
- else l=length(unlist(cvector))
- i = 1L
- j = 1L
- res = integer(l + 1L)
- tmp = integer(l)
- res[1] = root
- while (i < (l + 1)) {
- x = cvector[[res[i]]]
- if(nni)x = x[x > nTips]
- if (length(x) == 0) {
- res[i + 1] = tmp[j - 1]
- j = j - 1L
- }
- else {
- res[i + 1] = x[1]
- m = length(x)
- if (m > 1) {
- tmp[j:(j + m - 2L)] = x[-1]
- j = j + m - 1L
- }
- }
- i = i + 1
- }
- res[-1]
+orderNNI <- function (tree, nTips){
+ res = reorder(tree)$edge[,2]
+ res = res[res>nTips]
+ res
}
-
+
rooted.nni <- function(tree, data, eig, w, g, bf, rate, ll.0, INV,
control = pml.control(epsilon = 1e-08, maxit = 25, trace=0), ...){
ind0 = which(ll.0>0)
@@ -4199,7 +3727,8 @@ rooted.nni <- function(tree, data, eig, w, g, bf, rate, ll.0, INV,
optRootU <- function(t, tree, data, bf, g, w, eig, ll.0, k, INV, nh){
tree$edge.length = getEL1(t, nh)
- pml.fit(tree, data, bf=bf, g=g, w=w, eig=eig, ll.0=ll.0, k=k, INV=INV)
+ pml.fit4(tree, data, bf=bf, k=k, eig=eig, ll.0=ll.0, INV=INV, w=w, g=g)
+# pml.fit(tree, data, bf=bf, g=g, w=w, eig=eig, ll.0=ll.0, k=k, INV=INV)
}
@@ -4223,13 +3752,15 @@ rooted.nni <- function(tree, data, eig, w, g, bf, rate, ll.0, INV,
optEdgeU <- function(t, tree, data, bf, g, w, eig, ll.0, k, INV, nh){
tree$edge.length = getEL2(t, nh)
- pml.fit(tree, data, bf=bf, g=g, w=w, eig=eig, ll.0=ll.0, k=k, INV=INV)
+# pml.fit(tree, data, bf=bf, g=g, w=w, eig=eig, ll.0=ll.0, k=k, INV=INV)
+ pml.fit4(tree, data, bf=bf, k=k, eig=eig, ll.0=ll.0, INV=INV, w=w, g=g)
}
child = tree$edge[, 2]
parent = tree$edge[, 1]
- ll <- pml.fit(tree, data, bf=bf, k=k, eig=eig, ll.0=ll.0, INV=INV, w=w, g=g)
+# ll <- pml.fit(tree, data, bf=bf, k=k, eig=eig, ll.0=ll.0, INV=INV, w=w, g=g)
+ ll <- pml.fit4(tree, data, bf=bf, k=k, eig=eig, ll.0=ll.0, INV=INV, w=w, g=g)
llstart <- ll
eps=.00001
iter = 1
@@ -4241,11 +3772,10 @@ rooted.nni <- function(tree, data, eig, w, g, bf, rate, ll.0, INV,
cvector = allChildren(tree)
sibs <- Siblings(tree, 1:max(tree$edge))
- child2 = orderNNI(cvector, rootNode, nTips)
+ child2 = orderNNI(tree, nTips)
while(iter < 2){
ll2 <- pml.fit(tree, data, bf=bf, k=k, eig=eig, ll.0=ll.0, INV=INV, w=w, g=g)
-# browser()
nh=nodeHeight(tree)
@@ -4317,7 +3847,8 @@ rooted.nni <- function(tree, data, eig, w, g, bf, rate, ll.0, INV,
if(loli!=rootNode){
tree1$edge <- X1
tree1$edge.length = abs(nh[X1[,1]] - nh[X1[,2]])
- ll0 = pml.fit(tree1, data, bf=bf, g=g, w=w, eig=eig, ll.0=ll.0, k=k, INV=INV)
+# ll0 = pml.fit(tree1, data, bf=bf, g=g, w=w, eig=eig, ll.0=ll.0, k=k, INV=INV)
+ ll0 <- pml.fit4(tree1, data, bf=bf, k=k, eig=eig, ll.0=ll.0, INV=INV, w=w, g=g)
# cat("quartet", ll0, ch, dad, "\n")
}
@@ -4390,7 +3921,8 @@ rooted.nni <- function(tree, data, eig, w, g, bf, rate, ll.0, INV,
if( (ll3 - 1e-5*ll3) < ll2){
loli = rootNode
- ll2 <- pml.fit(tree, data, bf=bf, eig=eig, ll.0=ll.0, w=w, g=g)
+# ll2 <- pml.fit(tree, data, bf=bf, eig=eig, ll.0=ll.0, w=w, g=g)
+ ll2 <- pml.fit4(tree, data, bf=bf, k=k, eig=eig, ll.0=ll.0, INV=INV, w=w, g=g)
nh=nodeHeight(tree)
EL[tree$edge[,2]] = tree$edge.length
ind=0
@@ -4439,7 +3971,9 @@ rooted.nni <- function(tree, data, eig, w, g, bf, rate, ll.0, INV,
}
}
- ll2 <- pml.fit(tree, data, bf=bf, g=g, w=w, eig=eig, ll.0=ll.0, k=k, INV=INV)
+# ll2 <- pml.fit(tree, data, bf=bf, g=g, w=w, eig=eig, ll.0=ll.0, k=k, INV=INV)
+ ll2 <- pml.fit4(tree, data, bf=bf, k=k, eig=eig, ll.0=ll.0, INV=INV, w=w, g=g)
+
eps = (ll - ll2) / ll2
if(control$trace>1) cat(ll, " -> ", ll2, "\n")
if(control$trace>1) cat("swap:", nchanges)
diff --git a/R/treeManipulation.R b/R/treeManipulation.R
index 050bb63..d9d21ef 100644
--- a/R/treeManipulation.R
+++ b/R/treeManipulation.R
@@ -1040,6 +1040,7 @@ Siblings = function (x, node, include.self = FALSE)
root <- as.integer(parents[!match(parents, child, 0)][1])
res = vector("list", l)
ch = allChildren(x)
+ if(include.self) return(ch[ pvector[node] ])
k = 1
for(i in node){
if(i != root){
diff --git a/R/treedist.R b/R/treedist.R
index 1d68789..d4f2040 100644
--- a/R/treedist.R
+++ b/R/treedist.R
@@ -32,6 +32,25 @@ coph <- function(x){
}
+cophenetic.splits <- function(x){
+ labels <- attr(x, "labels")
+ X <- splits2design(x)
+ dm <- as.vector(X%*%attr(x, "weight"))
+ attr(dm, "Size") <- length(labels)
+ attr(dm, "Labels") <- labels
+ attr(dm, "Diag") <- FALSE
+ attr(dm, "Upper") <- FALSE
+ class(dm) <- "dist"
+ dm
+}
+
+
+cophenetic.networx <- function(x){
+ spl <- attr(x, "splits")
+ cophenetic.splits(spl)
+}
+
+
SHORTwise <- function (x, nTips, delete=FALSE)
{
v <- 1:nTips
@@ -232,7 +251,7 @@ mRF<-function(trees){
}
-RF.dist <- function (tree1, tree2=NULL, check.labels = TRUE)
+RF.dist <- function (tree1, tree2=NULL, check.labels = TRUE, rooted=FALSE)
{
if(class(tree1)=="multiPhylo" && is.null(tree2))return(mRF(tree1))
if(class(tree1)=="phylo" && class(tree2)=="multiPhylo")return(mRF2(tree1, tree2, check.labels))
@@ -241,7 +260,12 @@ RF.dist <- function (tree1, tree2=NULL, check.labels = TRUE)
r2 = is.rooted(tree2)
if(r1 != r2){
warning("one tree is unrooted, unrooted both")
+ }
+ if(!rooted){
+ if(r1) tree1<-unroot(tree1)
+ if(r2) tree2<-unroot(tree2)
}
+
if (check.labels) {
ind <- match(tree1$tip.label, tree2$tip.label)
if (any(is.na(ind)) | length(tree1$tip.label) !=
@@ -256,18 +280,15 @@ RF.dist <- function (tree1, tree2=NULL, check.labels = TRUE)
if(!r1 | !r2){
if(r1) tree1 = unroot(tree1)
if(r2) tree2 = unroot(tree2)
-# ref1 <- Ancestors(tree1, 1, "parent")
-# tree1 <- reroot(tree1, ref1)
-# ref2 <- Ancestors(tree2, 1, "parent")
-# tree2 <- reroot(tree2, ref2)
}
if(!is.binary.tree(tree1) | !is.binary.tree(tree2))warning("Trees are not binary!")
bp1 = bipart(tree1)
bp2 = bipart(tree2)
-
- bp1 <- SHORTwise(bp1, length(tree1$tip))
- bp2 <- SHORTwise(bp2, length(tree2$tip))
-
+ if(!rooted){
+ bp1 <- SHORTwise(bp1, length(tree1$tip))
+ bp2 <- SHORTwise(bp2, length(tree2$tip))
+ }
RF = sum(match(bp1, bp2, nomatch=0L)==0L) + sum(match(bp2, bp1, nomatch=0L)==0L)
RF
}
+
diff --git a/README.md b/README.md
index ef6aaef..e427f93 100644
--- a/README.md
+++ b/README.md
@@ -1,3 +1,9 @@
+[](https://travis-ci.org/KlausVigo/phangorn)
+[](http://cran.r-project.org/package=phangorn)
+[](http://www.gnu.org/licenses/gpl-2.0.html)
+[](http://cran.r-project.org/package=phangorn)
+
+
phangorn
========================================================
diff --git a/TODO b/TODO
new file mode 100644
index 0000000..190248f
--- /dev/null
+++ b/TODO
@@ -0,0 +1,15 @@
+1.99-4
+
+modelTest * check optimisation
+ * AIC, AICc, BIC
+
+optim.pml * start tree optimisation
+ * unique tree (multifurcations)
+
+2.0.0
+Rcpp, RcppArmadillo
+imprevements pmlPart, pmlCluster
+
+
+
+
diff --git a/build/vignette.rds b/build/vignette.rds
index 8dfe031..b5a2678 100644
Binary files a/build/vignette.rds and b/build/vignette.rds differ
diff --git a/inst/doc/Networx.R b/inst/doc/Networx.R
index c8bce36..dcd4933 100644
--- a/inst/doc/Networx.R
+++ b/inst/doc/Networx.R
@@ -1,9 +1,9 @@
-## ----, eval=TRUE---------------------------------------------------------
+## ---- eval=TRUE----------------------------------------------------------
library(phangorn)
data(Laurasiatherian)
data(yeast)
-## ----, eval=TRUE---------------------------------------------------------
+## ---- eval=TRUE----------------------------------------------------------
set.seed(1)
bs <- bootstrap.phyDat(yeast, FUN = function(x)nj(dist.hamming(x)),
bs=100)
@@ -13,31 +13,31 @@ tree <- plotBS(tree, bs, "phylogram")
cnet <- consensusNet(bs, .3)
plot(cnet, "2D", show.edge.label=TRUE)
-## ----, eval=FALSE--------------------------------------------------------
+## ---- eval=FALSE---------------------------------------------------------
# plot(cnet)
# # rotate 3d plot
# play3d(spin3d(axis=c(0,1,0), rpm=6), duration=10)
# # create animated gif file
# movie3d(spin3d(axis=c(0,1,0), rpm=6), duration=10)
-## ----, eval=TRUE---------------------------------------------------------
+## ---- eval=TRUE----------------------------------------------------------
dm <- dist.hamming(yeast)
nnet <- neighborNet(dm)
par("mar" = rep(2, 4))
plot(nnet, "2D")
-## ----, eval=TRUE---------------------------------------------------------
+## ---- eval=TRUE----------------------------------------------------------
nnet <- addConfidences(nnet, tree)
par("mar" = rep(2, 4))
plot(nnet, "2D", show.edge.label=TRUE)
-## ----, eval=TRUE---------------------------------------------------------
+## ---- eval=TRUE----------------------------------------------------------
tree2 <- rNNI(tree, 2)
tree2 <- addConfidences(tree2, tree)
# several support values are missing
plot(tree2, show.node.label=TRUE)
-## ----, eval=TRUE---------------------------------------------------------
+## ---- eval=TRUE----------------------------------------------------------
cnet <- nnls.networx(cnet, dm)
par("mar" = rep(2, 4))
plot(cnet, "2D", show.edge.label=TRUE)
diff --git a/inst/doc/Networx.Rmd b/inst/doc/Networx.Rmd
index 8acfaf9..4086c5a 100644
--- a/inst/doc/Networx.Rmd
+++ b/inst/doc/Networx.Rmd
@@ -11,7 +11,7 @@ vignette: >
---
-This tutorial gives a basic introduction on constructing phylogenetic networks and to add parameter to trees or networx using [phangorn](http://cran.r-project.org/web/packages/phangorn/) [@Schliep2011] in R.
+This tutorial gives a basic introduction on constructing phylogenetic networks and to add parameter to trees or networx using [phangorn](http://cran.r-project.org/package=phangorn) [@Schliep2011] in R.
Splits graph or phylogenetic networks are a nice way to display conflict data or summarize different trees. Here we present to popular networks, consensus networks [@Holland2004]
and neighborNet [@Bryant2004].
Often trees or networks are missing either edge weights or support values about the edges. We show how to improve a tree/networx 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 bc15b66..a589175 100644
--- a/inst/doc/Networx.html
+++ b/inst/doc/Networx.html
@@ -58,7 +58,7 @@ code > span.er { color: #ff0000; font-weight: bold; }
</div>
-<p>This tutorial gives a basic introduction on constructing phylogenetic networks and to add parameter to trees or networx using <a href="http://cran.r-project.org/web/packages/phangorn/">phangorn</a> <span class="citation">(Schliep 2011)</span> in R. Splits graph or phylogenetic networks are a nice way to display conflict data or summarize different trees. Here we present to popular networks, consensus networks <span class="citation">(Holland et al. 2004)</span> and neighborNet <span cl [...]
+<p>This tutorial gives a basic introduction on constructing phylogenetic networks and to add parameter to trees or networx using <a href="http://cran.r-project.org/package=phangorn">phangorn</a> <span class="citation">(Schliep 2011)</span> in R. Splits graph or phylogenetic networks are a nice way to display conflict data or summarize different trees. Here we present to popular networks, consensus networks <span class="citation">(Holland et al. 2004)</span> and neighborNet <span class="c [...]
<p>We first load the phangorn package and a few data sets we use in this vignette.</p>
<pre class="sourceCode r"><code class="sourceCode r"><span class="kw">library</span>(phangorn)</code></pre>
<pre><code>## Loading required package: ape</code></pre>
diff --git a/inst/doc/Trees.R b/inst/doc/Trees.R
index c307472..7bb7665 100644
--- a/inst/doc/Trees.R
+++ b/inst/doc/Trees.R
@@ -17,7 +17,7 @@ primates = read.phyDat("primates.dna", format="phylip", type="DNA")
###################################################
### code chunk number 3: Trees.Rnw:75-78
###################################################
-dm = dist.dna(as.DNAbin(primates))
+dm = dist.ml(primates)
treeUPGMA = upgma(dm)
treeNJ = NJ(dm)
diff --git a/inst/doc/Trees.Rnw b/inst/doc/Trees.Rnw
index 6af3b5a..8fba8fb 100644
--- a/inst/doc/Trees.Rnw
+++ b/inst/doc/Trees.Rnw
@@ -73,7 +73,7 @@ After reading in the alignment we can build a first tree with distance based met
After constructing a distance matrix we reconstruct a rooted tree with UPGMA and alternatively an unrooted tree using Neighbor Joining \cite{Saitou1987,Studier1988}.
<<echo=TRUE>>=
-dm = dist.dna(as.DNAbin(primates))
+dm = dist.ml(primates)
treeUPGMA = upgma(dm)
treeNJ = NJ(dm)
@
diff --git a/man/bootstrap.pml.Rd b/man/bootstrap.pml.Rd
index 123ab15..89b74c9 100644
--- a/man/bootstrap.pml.Rd
+++ b/man/bootstrap.pml.Rd
@@ -11,7 +11,7 @@ Bootstrap }
\usage{
bootstrap.pml(x, bs = 100, trees = TRUE, multicore=FALSE, ...)
bootstrap.phyDat(x, FUN, bs = 100, mc.cores = 1L, ...)
-plotBS(tree, BStrees, type="unrooted", bs.col="black", bs.adj=NULL, ...)
+plotBS(tree, BStrees, type="unrooted", bs.col="black", bs.adj=NULL, p=80, ...)
}
\arguments{
\item{x}{
@@ -50,6 +50,9 @@ color of bootstrap support labels.
\item{bs.adj}{
one or two numeric values specifying the horizontal and vertical justification of the bootstrap labels.
}
+ \item{p}{
+only plot support values higher than this percentage number (default is 80).
+}
}
\details{
It is possible that the bootstrap is performed in parallel, with help of the multicore package.
diff --git a/man/densiTree.Rd b/man/densiTree.Rd
index ac4faa9..c094ed5 100644
--- a/man/densiTree.Rd
+++ b/man/densiTree.Rd
@@ -50,6 +50,9 @@ Trees should be rooted, other wise the output may not make sense.
}
\references{
densiTree is inspired from the great \href{www.cs.auckland.ac.nz/~remco/DensiTree}{DensiTree} program of Remco Bouckaert.
+
+Remco R. Bouckaert (2010) DensiTree: making sense of sets of phylogenetic trees
+\emph{Bioinformatics}, \bold{26 (10)}, 1372-1373.
}
\author{
Klaus Schliep \email{klaus.schliep at gmail.com}
diff --git a/man/dist.hamming.Rd b/man/dist.hamming.Rd
index c275e48..df3dcce 100644
--- a/man/dist.hamming.Rd
+++ b/man/dist.hamming.Rd
@@ -2,6 +2,8 @@
\alias{dist.hamming}
\alias{dist.logDet}
\alias{dist.ml}
+\alias{readDist}
+\alias{writeDist}
%- Also NEED an '\alias' for EACH other topic documented here.
\title{Pairwise Distances from Sequences}
\description{
@@ -13,6 +15,8 @@ for nucleotide and amino acid models.
dist.hamming(x, ratio = TRUE, exclude="none")
dist.logDet(x)
dist.ml(x, model="JC69", exclude="none", bf=NULL, Q=NULL, ...)
+readDist(file)
+writeDist(dm, file="")
}
\arguments{
\item{x}{An object of class \code{phyDat}}
@@ -22,6 +26,8 @@ dist.ml(x, model="JC69", exclude="none", bf=NULL, Q=NULL, ...)
\item{bf}{A vector of base frequencies.}
\item{Q}{A vector containing the lower triangular part of the rate matrix.}
\item{\dots}{Further arguments passed to or from other methods.}
+ \item{file}{A file name.}
+ \item{dm}{A \code{dist} object.}
}
\value{
an object of class \code{dist}
diff --git a/man/modelTest.Rd b/man/modelTest.Rd
index ac9a8c0..4ee89ac 100644
--- a/man/modelTest.Rd
+++ b/man/modelTest.Rd
@@ -1,5 +1,6 @@
\name{modelTest}
\alias{modelTest}
+\alias{AICc}
\title{
ModelTest
}
@@ -8,7 +9,7 @@ Comparison of different substition models
}
\usage{
modelTest(object, tree=NULL, model = c("JC", "F81", "K80", "HKY", "SYM", "GTR"),
- G = TRUE, I = TRUE, k = 4, control = pml.control(epsilon = 1e-08, maxit = 3,
+ G = TRUE, I = TRUE, k = 4, control = pml.control(epsilon = 1e-08, maxit = 10,
trace = 1), multicore = FALSE)
}
%- maybe also 'usage' for other objects documented here.
@@ -29,10 +30,12 @@ This is only possible without GUI interface and under linux.
Only nucleotide models are tested so far.
}
\value{
-A data.frame containing the log-likelihood, AIC and BIC all tested models.
+A data.frame containing the log-likelihood, AIC, AICc and BIC all tested models.
The data.frame has an attributes "env" which is an environment which contains all the trees, the data and the calls to allow get the estimated models, e.g. as a starting point for further analysis (see example).
}
\references{
+Burnham, K. P. and Anderson, D. R (2002) \emph{Model selection and multimodel inference: a practical information-theoretic approach}. 2nd ed. Springer, New York
+
Posada, D. and Crandall, K.A. (1998) MODELTEST: testing the model of DNA substitution. \emph{Bioinformatics} \bold{14(9)}: 817-818
Posada, D. (2008) jModelTest: Phylogenetic Model Averaging. \emph{Molecular Biology and Evolution} \bold{25}: 1253-1256
diff --git a/man/parsimony.Rd b/man/parsimony.Rd
index de790ef..fc5180d 100644
--- a/man/parsimony.Rd
+++ b/man/parsimony.Rd
@@ -21,7 +21,7 @@ rearrangements or sub tree pruning and regrafting (SPR). \code{pratchet} impleme
\usage{
parsimony(tree, data, method="fitch", ...)
optim.parsimony(tree, data, method="fitch", cost=NULL, trace=1, rearrangements="SPR", ...)
-pratchet(data, start=NULL, method="fitch", maxit=100, k=5, trace=1, all=FALSE,
+pratchet(data, start=NULL, method="fitch", maxit=1000, k=10, trace=1, all=FALSE,
rearrangements="SPR", ...)
fitch(tree, data, site = "pscore")
sankoff(tree, data, cost = NULL, site = "pscore")
@@ -69,7 +69,7 @@ tree = NJ(dm)
parsimony(tree, Laurasiatherian)
treeRA <- random.addition(Laurasiatherian)
treeNNI <- optim.parsimony(tree, Laurasiatherian)
-treeRatchet <- pratchet(Laurasiatherian, start=tree)
+treeRatchet <- pratchet(Laurasiatherian, start=tree, maxit=100, k=5)
# assign edge length
treeRatchet <- acctran(treeRatchet, Laurasiatherian)
diff --git a/man/phyDat.Rd b/man/phyDat.Rd
index ab2e200..dacd2d5 100644
--- a/man/phyDat.Rd
+++ b/man/phyDat.Rd
@@ -1,6 +1,7 @@
\name{phyDat}
\alias{phyDat}
\alias{as.phyDat.DNAbin}
+\alias{as.phyDat.alignment}
\alias{as.data.frame.phyDat}
\alias{as.character.phyDat}
\alias{as.DNAbin.phyDat}
@@ -13,6 +14,7 @@
\alias{baseFreq}
\alias{cbind.phyDat}
\alias{c.phyDat}
+\alias{phyDat2alignment}
\title{Conversion among Sequence Formats}
\description{
These functions transform several DNA formats into the \code{phyDat} format.
@@ -23,10 +25,12 @@ phyDat(data, type = "DNA", levels = NULL, return.index=TRUE, ...)
read.phyDat(file, format="phylip", type="DNA", ...)
write.phyDat(x, file, format="phylip",...)
\method{as.phyDat}{DNAbin}(x, ...)
+\method{as.phyDat}{alignment}(x, type="DNA", ...)
\method{as.character}{phyDat}(x, allLevels = TRUE, ...)
\method{as.data.frame}{phyDat}(x, ...)
\method{as.DNAbin}{phyDat}(x, ...)
\method{subset}{phyDat}(x, subset, select, site.pattern = TRUE, ...)
+phyDat2alignment(x)
allSitePattern(n, levels=c("a","c","g","t"), names=NULL)
acgt2ry(obj)
baseFreq(obj, freq=FALSE, drop.unused.levels=FALSE)
@@ -55,12 +59,15 @@ baseFreq(obj, freq=FALSE, drop.unused.levels=FALSE)
If \code{type} "USER" a vector has to be give to \code{levels}.
For example c("a", "c", "g", "t", "-") would create a data object that
can be used in phylogenetic analysis with gaps as fifth state.
+There is a more detailed example for specifying "USER" defined data
+formats in the vignette "phangorn-specials".
+
\code{allSitePattern} returns all possible site patterns and can
be useful in simulation studies. For further details see the vignette
phangorn-specials.
\code{write.phyDat} calls the function write.dna or write.nexus.data and
-\code{read.phyDat} calls the function read.dna, read.aa or read.nexus.data
+\code{read.phyDat} calls the function \code{read.dna}, \code{read.aa} or \code{read.nexus.data}
see for more details over there.
You may import data directly with \code{\link[ape]{read.dna}} or \code{\link[ape]{read.nexus.data}}
@@ -70,8 +77,6 @@ The generic function \code{c} can be used to to combine sequences and \code{uniq
all unique sequences or unique haplotypes.
\code{acgt2ry} converts a \code{phyDat} object of nucleotides into an binary ry-coded dataset.
-
-There is a more detailed example for specifying USER defined data formats in the vignette advanced features.
}
\value{
The functions return an object of class \code{phyDat}.
@@ -79,8 +84,10 @@ The functions return an object of class \code{phyDat}.
\author{Klaus Schliep \email{klaus.schliep at gmail.com}}
-\seealso{ \code{\link{DNAbin}}, \code{\link{as.DNAbin}}, \code{\link{read.dna}}, \code{\link{read.aa}} and \code{\link{read.nexus.data}} and
-the example of \code{\link{pmlMix}} for the use of \code{allSitePattern}}
+\seealso{ \code{\link{DNAbin}}, \code{\link{as.DNAbin}}, \code{\link{read.dna}}, \code{\link{read.aa}}, \code{\link{read.nexus.data}}
+and the chapter 1 in the \code{vignette("phangorn-specials", package="phangorn")}
+and the example of \code{\link{pmlMix}} for the use of \code{allSitePattern}
+}
\examples{
data(Laurasiatherian)
class(Laurasiatherian)
diff --git a/man/plot.networx.Rd b/man/plot.networx.Rd
index ca2bc2a..faf69f7 100644
--- a/man/plot.networx.Rd
+++ b/man/plot.networx.Rd
@@ -105,7 +105,7 @@ set.seed(1)
tree1 = rtree(20, rooted=FALSE)
sp = as.splits(rNNI(tree1, n=10))
net = as.networx(sp)
-plot(net)
+plot(net, "2D")
\dontrun{
# also see example in consensusNet
example(consensusNet)
diff --git a/man/pml.fit.Rd b/man/pml.fit.Rd
index 8db390e..4c88e66 100644
--- a/man/pml.fit.Rd
+++ b/man/pml.fit.Rd
@@ -3,6 +3,7 @@
\alias{edQt}
\alias{pml.init}
\alias{pml.free}
+\alias{discrete.gamma}
\alias{lli}
%- Also NEED an '\alias' for EACH other topic documented here.
\title{
@@ -20,6 +21,7 @@ pml.init(data, k)
pml.free()
edQt(Q = c(1, 1, 1, 1, 1, 1), bf = c(0.25, 0.25, 0.25, 0.25))
lli(data, tree, ...)
+discrete.gamma(alpha, k)
}
%- maybe also 'usage' for other objects documented here.
\arguments{
@@ -27,6 +29,7 @@ lli(data, tree, ...)
\item{data}{An alignment, object of class \code{phyDat}.}
\item{bf}{Base frequencies.}
\item{shape}{Shape parameter of the gamma distribution.}
+ \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}{
diff --git a/man/splitsNetwork.Rd b/man/splitsNetwork.Rd
index c719184..979eca8 100644
--- a/man/splitsNetwork.Rd
+++ b/man/splitsNetwork.Rd
@@ -48,7 +48,7 @@ data(yeast)
dm = dist.ml(yeast)
fit = splitsNetwork(dm)
net = as.networx(fit)
-plot(net)
+plot(net, "2D")
write.nexus.splits(fit)
}
\keyword{ cluster }% at least one, from doc/KEYWORDS
diff --git a/man/treedist.Rd b/man/treedist.Rd
index 9b40ebf..3fab22d 100644
--- a/man/treedist.Rd
+++ b/man/treedist.Rd
@@ -8,13 +8,14 @@
}
\usage{
treedist(tree1, tree2, check.labels = TRUE)
-RF.dist(tree1, tree2=NULL, check.labels=TRUE)
+RF.dist(tree1, tree2=NULL, check.labels=TRUE, rooted=FALSE)
}
\arguments{
\item{tree1}{ A phylogenetic tree (class \code{phylo})
or vector of trees (an object of class \code{multiPhylo}). See details }
\item{tree2}{ A phylogenetic tree. }
\item{check.labels}{compares labels of the trees.}
+ \item{rooted}{take bipartitions for rooted trees into account, default is unrooting the trees.}
}
\value{
\code{treedist} returns a vector containing the following tree distance methods
diff --git a/src/fitch.c b/src/fitch.c
index 1280729..e34c8cb 100644
--- a/src/fitch.c
+++ b/src/fitch.c
@@ -31,31 +31,6 @@ void fitch_init(int *data, int *m, int *n, double *weights, int *nr)
}
-SEXP getData(SEXP n, SEXP k){
- int i, m=INTEGER(n)[0], l=INTEGER(k)[0];
- SEXP DAT, DAT2, RESULT;
- PROTECT(RESULT = allocVector(VECSXP, 2L));
- PROTECT(DAT = allocMatrix(INTSXP, m, l));
- PROTECT(DAT2 = allocMatrix(INTSXP, m, l));
- for(i=0; i< m*l; i++) INTEGER(DAT)[i] = data1[i];
- for(i=0; i< m*l; i++) INTEGER(DAT2)[i] = data2[i];
- SET_VECTOR_ELT(RESULT, 0, DAT);
- SET_VECTOR_ELT(RESULT, 1, DAT2);
- UNPROTECT(3);
- return(RESULT);
-}
-
-
-SEXP getWeight(SEXP n){
- int i, m=INTEGER(n)[0];
- SEXP RESULT;
- PROTECT(RESULT = allocVector(REALSXP, m));
- for(i=0; i<m; i++) REAL(RESULT)[i] = weight[i];
- UNPROTECT(1);
- return(RESULT);
-}
-
-
int bitcount(int x){
int count;
for (count=0; x != 0; x>>=1)
@@ -389,13 +364,6 @@ void fitchN(int *dat1, int *dat2, int *nr){
}
}
-// raus
-void fitchN2(int *res, int *dat, int *node, int *edge, int *nr, int *nl) {
- int i;
- for (i=0; i< *nl; i++) {
- fitchN(&res[(node[i]-1L) * (*nr)], &dat[(edge[i]-1L) * (*nr)], nr);
- }
-}
// MPR reconstruction nicht immer gleiches ergebnis
void fitchTriplet(int *res, int *dat1, int *dat2, int *dat3, int *nr)
@@ -673,13 +641,11 @@ SEXP FNALL5(SEXP nrx, SEXP node, SEXP edge, SEXP l, SEXP mx, SEXP my, SEXP root)
edge2 = (int *) R_alloc(2L * *n, sizeof(int));
node2 = (int *) R_alloc(2L * *n, sizeof(int));
pc = (int *) R_alloc(2L * *n, sizeof(int));
-
-// pars = (int *) R_alloc(*nr, sizeof(int)); // raus
+
pvtmp2 = (double *) R_alloc(m, sizeof(double));
-
- PROTECT(pvec = allocVector(REALSXP, m));
-
+ PROTECT(pvec = allocVector(REALSXP, m));
pvtmp = REAL(pvec);
+
for(i=0; i<m; i++){
pvtmp[i] = 0.0;
pvtmp2[i] = 0.0;
@@ -717,7 +683,7 @@ void fitchquartet(int *dat1, int *dat2, int *dat3, int *dat4, int *nr, double *w
}
}
-// weight raus double *weight,
+// weight raus
void fitchQuartet(int *index, int *n, int *nr, double *psc1, double *psc2, double *weight, double *res){
int i, e1, e2, e3, e4;
for(i=0; i<*n; i++){
diff --git a/src/ml.c b/src/ml.c
index 5e094b8..3c20579 100644
--- a/src/ml.c
+++ b/src/ml.c
@@ -28,18 +28,16 @@ double one = 1.0, zero = 0.0;
int ONE = 1L;
const double ScaleEPS = 1.0/4294967296.0;
const double ScaleMAX = 4294967296.0;
+const double LOG_SCALE_EPS = -22.18070977791824915926;
// 2^64 = 18446744073709551616
-static double *LL, *ROOT;
+static double *LL;
static int *SCM, *XXX;
-
void ll_free(){
free(LL);
free(SCM);
- free(ROOT);
-// free(XX);
}
@@ -52,7 +50,6 @@ void ll_init(int *nr, int *nTips, int *nc, int *k)
{
int i;
LL = (double *) calloc(*nr * *nc * *k * *nTips, sizeof(double));
- ROOT = (double *) calloc(*nr * *nc * *k, sizeof(double));
SCM = (int *) calloc(*nr * *k * *nTips, sizeof(int)); // * 2L
for(i =0; i < (*nr * *k * *nTips); i++) SCM[i] = 0L;
}
@@ -62,7 +59,6 @@ void ll_init(int *nr, int *nTips, int *nc, int *k)
void ll_free2(){
free(LL);
free(SCM);
- free(ROOT);
free(XXX);
}
@@ -71,7 +67,6 @@ 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));
- ROOT = (double *) calloc(*nr * *nc * *k, 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;
@@ -179,12 +174,10 @@ SEXP getPM(SEXP eig, SEXP nc, SEXP el, SEXP w){
void lll(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 ni, ei, j, i, rc; // R_len_t i, n = length(node);
double *rtmp, *P;
-
ni = -1;
rc = *nr * *nc;
rtmp = (double *) R_alloc(*nr * *nc, sizeof(double));
P = (double *) R_alloc(*nc * *nc, sizeof(double));
-
for(j=0; j < *nr; j++) scaleTmp[j] = 0L;
for(i = 0; i < n; i++) {
getP(eva, eve, evei, *nc, el[i], g, P);
@@ -215,12 +208,10 @@ void lll(SEXP dlist, double *eva, double *eve, double *evei, double *el, double
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);
double *rtmp, *P;
-
ni = -1;
rc = *nr * *nc;
rtmp = (double *) R_alloc(*nr * *nc, sizeof(double));
P = (double *) R_alloc(*nc * *nc, sizeof(double));
-
for(j=0; j < *nr; j++) scaleTmp[j] = 0L;
for(i = 0; i < n; i++) {
getP(eva, eve, evei, *nc, el[i], g, P);
@@ -245,8 +236,6 @@ void lll0(int *X, double *eva, double *eve, double *evei, double *el, double g,
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){
@@ -256,9 +245,7 @@ void lll3(SEXP dlist, double *eva, double *eve, double *evei, double *el, double
rc = *nr * *nc;
rtmp = (double *) R_alloc(*nr * *nc, sizeof(double));
P = (double *) R_alloc(*nc * *nc, sizeof(double));
-
for(j=0; j < *nr; j++) scaleTmp[j] = 0L;
-
for(i = 0; i < n; i++) {
getP(eva, eve, evei, *nc, el[i], g, P);
ei = edge[i];
@@ -328,15 +315,11 @@ SEXP PML_NEW(SEXP EL, SEXP W, SEXP G, SEXP NR, SEXP NC, SEXP K, SEXP eig, SEXP b
double *el=REAL(EL), *bfreq=REAL(bf), *contr=REAL(contrast);
double *g=REAL(G), *tmp, logScaleEPS;
SEXP TMP;
-
double *eva, *eve, *evei;
-
eva = REAL(VECTOR_ELT(eig, 0));
eve = REAL(VECTOR_ELT(eig, 1));
- evei = REAL(VECTOR_ELT(eig, 2));
-
+ evei = REAL(VECTOR_ELT(eig, 2));
SC = (int *) R_alloc(nr * k, sizeof(int));
-
PROTECT(TMP = allocMatrix(REALSXP, nr, k)); // changed
tmp=REAL(TMP);
for(i=0; i<(k*nr); i++)tmp[i]=0.0;
@@ -355,8 +338,7 @@ if(*nthreads <= 1){ *nthreads=1; }else{ *nthreads=(*nthreads < omp_get_max_threa
}
logScaleEPS = log(ScaleEPS);
- for(i=0; i<(k*nr); i++) tmp[i] = logScaleEPS * SC[i] + log(tmp[i]); //log
-
+ for(i=0; i<(k*nr); i++) tmp[i] = logScaleEPS * SC[i] + log(tmp[i]); //log
UNPROTECT(1);
return TMP;
}
@@ -368,13 +350,10 @@ SEXP PML3(SEXP dlist, SEXP EL, SEXP W, SEXP G, SEXP NR, SEXP NC, SEXP K, SEXP ei
double *g=REAL(G), *tmp, logScaleEPS;
SEXP TMP;
double *eva, *eve, *evei;
-
eva = REAL(VECTOR_ELT(eig, 0));
eve = REAL(VECTOR_ELT(eig, 1));
evei = REAL(VECTOR_ELT(eig, 2));
-
SC = (int *) R_alloc(nr * k, sizeof(int));
-
PROTECT(TMP = allocMatrix(REALSXP, nr, k)); // changed
tmp=REAL(TMP);
for(i=0; i<(k*nr); i++)tmp[i]=0.0;
@@ -394,12 +373,10 @@ SEXP PML0(SEXP dlist, SEXP EL, SEXP W, SEXP G, SEXP NR, SEXP NC, SEXP K, SEXP ei
int nTips = INTEGER(NTips)[0], *SC;
double *g=REAL(G), *tmp, logScaleEPS;
SEXP TMP;
- double *eva, *eve, *evei;
-
+ double *eva, *eve, *evei;
eva = REAL(VECTOR_ELT(eig, 0));
eve = REAL(VECTOR_ELT(eig, 1));
- evei = REAL(VECTOR_ELT(eig, 2));
-
+ evei = REAL(VECTOR_ELT(eig, 2));
SC = (int *) R_alloc(nr * k, sizeof(int));
PROTECT(TMP = allocMatrix(REALSXP, nr, k)); // changed
@@ -408,60 +385,13 @@ SEXP PML0(SEXP dlist, SEXP EL, SEXP W, SEXP G, SEXP NR, SEXP NC, SEXP K, SEXP ei
indLL = nr * nc * nTips;
for(i=0; i<k; i++){
lll(dlist, eva, eve, evei, REAL(EL), g[i], &nr, &nc, INTEGER(node), INTEGER(edge), nTips, REAL(contrast), INTEGER(nco)[0], INTEGER(N)[0], &SC[nr * i], REAL(bf), &tmp[i*nr], &LL[indLL *i]);
- }
-
+ }
logScaleEPS = log(ScaleEPS);
for(i=0; i<(k*nr); i++) tmp[i] = logScaleEPS * SC[i] + log(tmp[i]);
-
- UNPROTECT(1);
- return TMP;
-}
-
-
-// replace child with LL
-void moveLLNew(double *LL, double *child, double *P, int *nr, int *nc, double *tmp, int *LLSC, int *CSC){
- int j, a;
- F77_CALL(dgemm)(transa, transb, nr, nc, nc, &one, child, nr, P, nc, &zero, tmp, nr);
- for(j=0; j<(*nc * *nr); j++) LL[j]/=tmp[j]; // new child
- F77_CALL(dgemm)(transa, transb, nr, nc, nc, &one, LL, nr, P, nc, &zero, tmp, nr);
- for(j=0; j<(*nc * *nr); j++) child[j] *= tmp[j];
- for(j=0; j<*nr; j++){
- a = LLSC[j];
- LLSC[j] -= CSC[j];
- CSC[j] = a;
- }
-}
-
-
-void moveLL0(double *LL, double *child, double *P, int *nr, int *nc, double *tmp){
- int j;
- F77_CALL(dgemm)(transa, transb, nr, nc, nc, &one, child, nr, P, nc, &zero, tmp, nr);
- for(j=0; j<(*nc * *nr); j++) LL[j]/=tmp[j]; // new child
- F77_CALL(dgemm)(transa, transb, nr, nc, nc, &one, LL, nr, P, nc, &zero, tmp, nr);
- for(j=0; j<(*nc * *nr); j++) child[j] *= tmp[j];
-}
-
-
-void moveLL2(int *loli, int *nloli, double *eva, double *eve, double *evi, double *el, double *g, int *nr, int *nc, int *k, int *ntips){
- double *tmp, *P;
- int j;
- tmp = (double *) R_alloc(*nr * *nc, sizeof(double));
- P = (double *) R_alloc(*nc * *nc, sizeof(double));
- for(j = 0; j < *k; j++){
- getP(eva, eve, evi, *nc, el[0], g[j], P);
- moveLLNew(&LL[LINDEX2(*loli, j)], &LL[LINDEX2(*nloli, j)], P, nr, nc, tmp,
- &SCM[LINDEX3(*loli, j)], &SCM[LINDEX3(*nloli, j)]);
- }
+ UNPROTECT(1);
+ return TMP;
}
-void moveLLtmp(double *LL, double *child, double *P, int *nr, int *nc, double *tmp){
- int j;
- F77_CALL(dgemm)(transa, transb, nr, nc, nc, &one, child, nr, P, nc, &zero, tmp, nr);
- for(j=0; j<(*nc * *nr); j++) LL[j]/=tmp[j];
- F77_CALL(dgemm)(transa, transb, nr, nc, nc, &one, LL, nr, P, nc, &zero, tmp, nr);
- for(j=0; j<(*nc * *nr); j++) LL[j]=tmp[j] * child[j];
-}
-
void moveLL5(double *LL, double *child, double *P, int *nr, int *nc, double *tmp){
int j;
@@ -482,25 +412,15 @@ SEXP moveloli(SEXP CH, SEXP PA, SEXP eig, SEXP EL, SEXP W, SEXP G,
double *eva, *eve, *evei, *tmp, *P;
tmp = (double *) R_alloc(nr * nc, sizeof(double));
P = (double *) R_alloc(nc * nc, sizeof(double));
-
-// SEXP X, RESULT;
-// PROTECT(RESULT = allocVector(VECSXP, k));
-
+
eva = REAL(VECTOR_ELT(eig, 0));
eve = REAL(VECTOR_ELT(eig, 1));
evei = REAL(VECTOR_ELT(eig, 2));
for(i = 0; i < k; i++){
-// PROTECT(X = allocMatrix(REALSXP, nr, nc));
getP(eva, eve, evei, nc, el, g[i], P);
moveLL5(&LL[LINDEX(ch, i)], &LL[LINDEX(pa, i)], P, &nr, &nc, tmp);
-// blub = LINDEX(ch, i);
-// for(j=0; j< (nr*nc); j++) REAL(X)[j] = LL[blub+j];
-// SET_VECTOR_ELT(RESULT, i, X);
-// UNPROTECT(1);
}
-// UNPROTECT(1); //RESULT
-// return(RESULT);
return ScalarReal(1L);
}
@@ -603,10 +523,8 @@ SEXP moveDad(SEXP dlist, SEXP PA, SEXP CH, SEXP eig, SEXP EVI, SEXP EL, SEXP W,
else{
for(i = 0; i < k; i++){
PROTECT(X = allocMatrix(REALSXP, nr, nc));
- getP(eva, eve, evei, nc, el, g[i], P);
-// helpDAD2(double *dad, int *child, double *contrast, double *P, int nr, int nc, int nco, double *res)
+ getP(eva, eve, evei, nc, el, g[i], P);
helpDAD5(&LL[LINDEX(pa, i)], INTEGER(VECTOR_ELT(dlist, ch-1L)), contrast, P, nr, nc, nco, tmp);
-// void helpPrep2(double *dad, int *child, double *contrast, double *evi, int nr, int nc, int nrs, double *res)
helpPrep2(&LL[LINDEX(pa, i)], INTEGER(VECTOR_ELT(dlist, ch-1L)), contrast2, evi, nr, nc, nco, REAL(X)); //;
SET_VECTOR_ELT(RESULT, i, X);
UNPROTECT(1);
@@ -629,10 +547,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];
}
-
+// in optimEdgeOld
SEXP updateLL(SEXP dlist, SEXP PA, SEXP CH, SEXP eig, SEXP EL, SEXP W, SEXP G, SEXP NR,
- SEXP NC, SEXP NTIPS, SEXP CONTRAST, SEXP NCO){
-// SEXP RESULTS, X;
+ SEXP NC, SEXP NTIPS, SEXP CONTRAST, SEXP NCO){
int i, k=length(W);
int nc=INTEGER(NC)[0], nr=INTEGER(NR)[0], ntips=INTEGER(NTIPS)[0]; //, j, blub
int pa=INTEGER(PA)[0], ch=INTEGER(CH)[0], nco =INTEGER(NCO)[0];
@@ -642,8 +559,6 @@ SEXP updateLL(SEXP dlist, SEXP PA, SEXP CH, SEXP eig, SEXP EL, SEXP W, SEXP G, S
tmp = (double *) R_alloc(nr * nc, sizeof(double));
P = (double *) R_alloc(nc * nc, sizeof(double));
-// PROTECT(RESULT = allocVector(VECSXP, k))
-
eva = REAL(VECTOR_ELT(eig, 0));
eve = REAL(VECTOR_ELT(eig, 1));
evei = REAL(VECTOR_ELT(eig, 2));
@@ -660,7 +575,27 @@ SEXP updateLL(SEXP dlist, SEXP PA, SEXP CH, SEXP eig, SEXP EL, SEXP W, SEXP G, S
}
}
return ScalarReal(1L);
-// return(NULL);
+}
+
+
+void updateLL2(SEXP dlist, int pa, int ch, double *eva, double *eve, double*evei,
+ double el, double *w, double *g, int nr,
+ int nc, int ntips, double *contrast, int nco, int k,
+ double *tmp, double *P){
+ int i;
+
+ if(ch>ntips){
+ for(i = 0; i < k; i++){
+ getP(eva, eve, evei, nc, el, g[i], P);
+ goDown(&LL[LINDEX(pa, i)], &LL[LINDEX(ch, i)], P, nr, nc, tmp);
+ }
+ }
+ else{
+ for(i = 0; i < k; i++){
+ getP(eva, eve, evei, nc, el, g[i], P);
+ goUp(&LL[LINDEX(pa, i)], INTEGER(VECTOR_ELT(dlist, ch-1L)), contrast, P, nr, nc, nco, tmp);
+ }
+ }
}
@@ -700,52 +635,6 @@ SEXP extractScale(SEXP CH, SEXP W, SEXP G, SEXP NR, SEXP NC, SEXP NTIPS){
}
-// old version
-void moveLL(double *LL, double *child, double *P, int *nr, int *nc, double *tmp){
- int j;
- F77_CALL(dgemm)(transa, transb, nr, nc, nc, &one, child, nr, P, nc, &zero, tmp, nr);
- for(j=0; j<(*nc * *nr); j++) LL[j]/=tmp[j];
- F77_CALL(dgemm)(transa, transb, nr, nc, nc, &one, LL, nr, P, nc, &zero, tmp, nr);
- for(j=0; j<(*nc * *nr); j++) tmp[j] *= child[j];
-}
-
-
-void helpDAD3(double *dad, double *child, double *P, int *nr, int *nc, double *res){
- F77_CALL(dgemm)(transa, transb, nr, nc, nc, &one, child, nr, P, nc, &zero, res, nr);
- for(int j=0; j<(*nc * *nr); j++) dad[j]/=res[j];
-}
-
-void getDAD3(int *dad, int *child, double *eva, double *eve, double *evi, double *el, double *g, int *nr, int *nc, int *k, int *ntips){
- double *tmp, *P;
- int j;
- tmp = (double *) R_alloc(*nr * *nc, sizeof(double));
- P = (double *) R_alloc(*nc * *nc, sizeof(double));
- for(j = 0; j < *k; j++){
- getP(eva, eve, evi, *nc, el[0], g[j], P);
- helpDAD3(&LL[LINDEX2(*dad, j)], &LL[LINDEX2(*child, j)], P, nr, nc, tmp);
- }
-}
-
-
-// children
-void helpDAD4(double *dad, int *child, double *contrast, double *P, int *nr, int *nc, int *nco, double *res){
- matp(child, contrast, P, nr, nc, nco, res);
- for(int j=0; j<(*nc * *nr); j++) res[j]=dad[j]/res[j];
-}
-
-
-void getDAD4(int *dad, int *child, double *contrast, double *eva, double *eve, double *evi, double *el, double *g, int *nr, int *nc, int *nco, int *k, int *ntips){
- double *tmp, *P;
- int j;
- tmp = (double *) R_alloc(*nr * *nc, sizeof(double));
- P = (double *) R_alloc(*nc * *nc, sizeof(double));
- for(j = 0; j < *k; j++){
- getP(eva, eve, evi, *nc, el[0], g[j], P);
- helpDAD4(&LL[LINDEX2(*dad, j)], child, contrast, P, nr, nc, nco, tmp);
- }
-}
-
-
// dad / child * P
void helpDAD(double *dad, double *child, double *P, int nr, int nc, double *res){
F77_CALL(dgemm)(transa, transb, &nr, &nc, &nc, &one, child, &nr, P, &nc, &zero, res, &nr);
@@ -768,10 +657,8 @@ SEXP getDAD(SEXP dad, SEXP child, SEXP P, SEXP nr, SEXP nc){
return(RESULT);
}
-
-
+/*
// SEXP mixture of prepFS, prepFS & getPrep
-// while loop for moveLL mit LINDEX 1050 auf 900 ?
void prepFS(double *XX, int ch, int pa, double *eva, double *eve, double *evi, double el, double *g, int nr, int nc, int ntips, int k){
int i;
double *P, *tmp;
@@ -783,8 +670,7 @@ void prepFS(double *XX, int ch, int pa, double *eva, double *eve, double *evi, d
helpPrep(&LL[LINDEX(ch, i)], &LL[LINDEX(pa, i)], eve, evi, nr, nc, tmp, &XX[i*nr*nc]);
}
}
-
-
+*/
SEXP getPrep(SEXP dad, SEXP child, SEXP eve, SEXP evi, SEXP nr, SEXP nc){
R_len_t i, n=length(dad);
@@ -803,23 +689,18 @@ SEXP getPrep(SEXP dad, SEXP child, SEXP eve, SEXP evi, SEXP nr, SEXP nc){
return(RESULT);
}
-
-
-
+/*
void prepFSE(double *XX, int *ch, int pa, double *eva, double *eve, double *evi, double el, double *g, int nr, int nc, int ntips, int k, double *contrast, double *contrast2, int ncs){
int i;
double *P; //, *tmp
-// tmp = (double *) R_alloc(nr * nc, sizeof(double));
P = (double *) R_alloc(nc * nc, sizeof(double));
for(i=0; i<k; i++){
getP(eva, eve, evi, nc, el, g[i], P);
helpDAD2(&ROOT[i * nr * nc], ch, contrast, P, nr, nc, ncs, &LL[LINDEX(pa, i)]);
-// helpDAD(&ROOT[i * nr * nc], &LL[LINDEX(ch, i)], P, nr, nc, &LL[LINDEX(pa, i)]);
-// helpPrep(&LL[LINDEX(ch, i)], &LL[LINDEX(pa, i)], eve, evi, nr, nc, tmp, &XX[i*nr*nc]);
helpPrep2(&LL[LINDEX(pa, i)], ch, contrast2, evi, nr, nc, ncs, &XX[i*nr*nc]);
}
}
-
+*/
SEXP getSCM(SEXP kk, SEXP nrx, SEXP nTips){
int j, nr = INTEGER(nrx)[0], ntips = INTEGER(nTips)[0], k = INTEGER(kk)[0]-1L;
@@ -840,52 +721,14 @@ SEXP getLL(SEXP ax, SEXP bx, SEXP nrx, SEXP ncx, SEXP nTips){
return(RES);
}
-SEXP getROOT(SEXP ax, SEXP nrx, SEXP ncx){
- int j, nc = INTEGER(ncx)[0], nr = INTEGER(nrx)[0], a = INTEGER(ax)[0];
- SEXP RES;
- PROTECT(RES = allocMatrix(REALSXP, nr, nc));
- for(j=0; j<(nr*nc); j++) REAL(RES)[j] = ROOT[j + a * nr * nc];
- UNPROTECT(1);
- return(RES);
-}
-
-
-// ch * (pa %*% P)
-void getMI(int child, int parent, double el, double *eva, double *eve, double *evi, double *g, int nr, int nc, int k, int ntips){
- int i, a, j;
- double *P;
- P = (double *) R_alloc(nc * nc, sizeof(double));
- for(i=0; i<k; i++){
- getP(eva, eve, evi, nc, el, g[i], P);
- a = i * nr * nc;
- F77_CALL(dgemm)(transa, transb, &nr, &nc, &nc, &one, &LL[LINDEX(parent, i)], &nr, P, &nc, &zero, &ROOT[a], &nr);
- for(j=0; j<(nc * nr); j++) ROOT[a + j]*=LL[LINDEX(child, i) + j];
- }
-}
-
-// (ch %*% P) * pa
-void getME(int *child, int parent, double el, double *eva, double *eve, double *evi, double *g, int nr,
- int nc, int k, int ntips, double *contrast, int nrs){
- int i, a, j;
- double *P;
- P = (double *) R_alloc(nc * nc, sizeof(double));
- for(i=0; i<k; i++){
- getP(eva, eve, evi, nc, el, g[i], P);
- a = i * nr * nc;
- matp(child, contrast, P, &nr, &nc, &nrs, &ROOT[a]);
- for(j=0; j<(nc * nr); j++) ROOT[a + j]*=LL[LINDEX(parent, i) + j];
- }
-}
-
void NR55(double *eva, int nc, double el, double *w, double *g, SEXP X, int ld, int nr, double *f, double *res){
int i, j, k;
double *tmp;
tmp = (double *) R_alloc(nc, sizeof(double));
-
- for(k=0; k<nr; k++) res[k] = 0.0;
-
+ for(k=0; k<nr; k++) res[k] = 0.0;
for(j=0;j<ld;j++){
+// for(i=0; i<nc ;i++) tmp[i] = w[j] * (eva[i] * g[j] * el) * exp(eva[i] * g[j] * el);
for(i=0; i<nc ;i++) tmp[i] = (eva[i] * g[j] * el) * exp(eva[i] * g[j] * el);
F77_CALL(dgemv)(transa, &nr, &nc, &w[j], REAL(VECTOR_ELT(X, j)), &nr, tmp, &ONE, &one, res, &ONE);
}
@@ -918,9 +761,37 @@ void NR66(double *eva, int nc, double el, double *w, double *g, SEXP X, int ld,
for(i=0; i<nc ;i++) tmp[i] = exp(eva[i] * g[j] * el);
// alpha = w[j]
F77_CALL(dgemv)(transa, &nr, &nc, &w[j], REAL(VECTOR_ELT(X, j)), &nr, tmp, &ONE, &one, res, &ONE);
- }
+ }
+}
+
+
+
+//void NR55(double *eva, int nc, double el, double *w, double *g, SEXP X, int ld, int nr, double *f, double *res)
+ void NR77(double *eva, int nc, double el, double *w, double *g, double *X, int ld, int nr, double *f, double *res){
+ int i, j, k;
+ double *tmp;
+ tmp = (double *) R_alloc(nc, sizeof(double));
+ for(k=0; k<nr; k++) res[k] = 0.0;
+ int nrc = (nc+1L) * nr;
+ for(j=0;j<ld;j++){
+ for(i=0; i<nc ;i++) tmp[i] = (eva[i] * g[j] * el) * exp(eva[i] * g[j] * el);
+ F77_CALL(dgemv)(transa, &nr, &nc, &w[j], &X[j*nrc], &nr, tmp, &ONE, &one, res, &ONE);
+ }
+ for(i=0; i<nr ;i++) res[i]/=f[i];
+
}
+//void NR66(double *eva, int nc, double el, double *w, double *g, SEXP X, int ld, int nr, double *res)
+void NR88(double *eva, int nc, double el, double *w, double *g, double *X, int ld, int nr, double *res){
+ int i, j;
+ double *tmp; //*res, *dF,
+ tmp = (double *) R_alloc(nc, sizeof(double));
+ for(j=0;j<ld;j++){
+ for(i=0; i<nc ;i++) tmp[i] = exp(eva[i] * g[j] * el);
+ // alpha = w[j]
+ F77_CALL(dgemv)(transa, &nr, &nc, &w[j], &X[j*nr*nc], &nr, tmp, &ONE, &one, res, &ONE);
+ }
+}
// in ancestral.pml
@@ -1014,9 +885,6 @@ SEXP FS4(SEXP eig, SEXP nc, SEXP el, SEXP w, SEXP g, SEXP X, SEXP dad, SEXP chil
NR55(eva, ncx-1L, edle, ws, gs, X, INTEGER(ld)[0], nrx, f, tmp);
ll=0.0;
lll=0.0;
-// for(i=0; i<nrx ;i++) ll+=wgt[i]*tmp[i];
-// for(i=0; i<nrx ;i++) lll+=wgt[i]*tmp[i]*tmp[i];
-
for(i=0; i<nrx ;i++){
y = wgt[i]*tmp[i];
ll+=y;
@@ -1075,11 +943,11 @@ SEXP FS5(SEXP eig, SEXP nc, SEXP el, SEXP w, SEXP g, SEXP X, SEXP ld, SEXP nr, S
PROTECT(RESULT = allocVector(REALSXP, 3));
edle = REAL(el)[0];
- for(i=0; i<nrx; i++)f[i] = REAL(f0)[i];
+ for(i=0; i<nrx; i++)f[i] = REAL(f0)[i]; //memcpy
NR66(eva, ncx, edle, ws, gs, X, INTEGER(ld)[0], nrx, f); // ncx-1L !!!
for(i=0; i<nrx ;i++) l0 += wgt[i] * log(f[i]);
- while ( (eps > 1e-05) && (k < 5) ) {
+ while ( (eps > 1e-05) && (k < 10) ) {
if(scalep>0.6){
NR55(eva, ncx-1L, edle, ws, gs, X, INTEGER(ld)[0], nrx, f, tmp);
ll=0.0;
@@ -1123,12 +991,229 @@ SEXP FS5(SEXP eig, SEXP nc, SEXP el, SEXP w, SEXP g, SEXP X, SEXP ld, SEXP nr, S
lll=0.0;
for(i=0; i<nrx ;i++) lll+=wgt[i]*tmp[i]*tmp[i];
REAL(RESULT)[0] = edle;
- REAL(RESULT)[1] = 1.0/ lll; // variance
- REAL(RESULT)[2] = l0;
+ REAL(RESULT)[1] = 1.0 / lll; // variance
+ REAL(RESULT)[2] = l0; //l0
UNPROTECT(1);
return (RESULT);
}
+void fs3(double *eva, int nc, double el, double *w, double *g, double *X, int ld, int nr, double *weight,
+ double *f0, double *res)
+{
+ double *tmp, *f, edle, ledle, newedle, eps=10;
+ double ll=0.0, lll, delta=0.0, scalep = 1.0, l1=0.0, l0=0.0;
+ double y;
+ int i, k=0;
+ tmp = (double *) R_alloc(nr, sizeof(double));
+ f = (double *) R_alloc(nr, sizeof(double));
+ edle = el; //REAL(el)[0];
+
+ for(i=0; i<nr; i++)f[i] = f0[i];
+ NR88(eva, nc, edle, w, g, X, ld, nr, f); // nc-1L !!!
+ for(i=0; i<nr ;i++) l0 += weight[i] * log(f[i]);
+
+ while ( (eps > 1e-05) && (k < 10) ) {
+ if(scalep>0.6){
+ NR77(eva, nc-1L, edle, w, g, X, ld, nr, f, tmp);
+ ll=0.0;
+ lll=0.0;
+ for(i=0; i<nr ;i++){
+ y = weight[i]*tmp[i];
+ ll+=y;
+ lll+=y*tmp[i];
+ }
+ delta = ((ll/lll) < 3) ? (ll/lll) : 3;
+ } // end if
+ ledle = log(edle) + scalep * delta;
+ newedle = exp(ledle);
+// some error handling avoid too big small edges & too big steps
+ if (newedle > 10.0) newedle = 10.0;
+ if (newedle < 1e-8) newedle = 1e-8; // 1e-8 phyML
+
+ for(i=0; i<nr; i++)f[i] = f0[i];
+ NR88(eva, nc, newedle, w, g, X, ld, nr, f);
+ l1 = 0.0;
+ for(i=0; i<nr ;i++) l1 += weight[i] * log(f[i]);
+ eps = l1 - l0;
+// some error handling
+ if (eps < 0 || ISNAN(eps)) {
+ if (ISNAN(eps))eps = 0;
+ else {
+ scalep = scalep/2.0;
+ eps = 1.0;
+ }
+ newedle = edle;
+ l1 = l0;
+ }
+ else scalep = 1.0;
+ edle=newedle;
+ l0 = l1;
+ k ++;
+ }
+// variance
+ // NR555(eva, ncx-1L, edle, ws, gs, X, INTEGER(ld)[0], nrx, f, tmp);
+ // lll=0.0;
+ // for(i=0; i<nrx ;i++) lll+=wgt[i]*tmp[i]*tmp[i];
+ res[0] = edle;
+ res[1] = ll;
+ res[2] = l0; //l0
+}
+
+
+SEXP optE(SEXP PARENT, SEXP CHILD, SEXP ANC, SEXP eig, SEXP EVI, SEXP EL,
+ SEXP W, SEXP G, SEXP NR, SEXP NC, SEXP NTIPS, SEXP CONTRAST,
+ SEXP CONTRAST2, SEXP NCO,
+ SEXP BLUB, SEXP dlist, SEXP WEIGHT, SEXP F0){
+ int i, k=length(W), h, j, n=length(PARENT), m, lEL=length(EL);
+ int nc=INTEGER(NC)[0], nr=INTEGER(NR)[0], ntips=INTEGER(NTIPS)[0];
+ int *parent=INTEGER(PARENT), *child=INTEGER(CHILD), *anc=INTEGER(ANC); //
+ int loli, nco =INTEGER(NCO)[0]; //=INTEGER(LOLI)[0]
+ double *weight=REAL(WEIGHT), *f0=REAL(F0), *w=REAL(W);
+ double *g=REAL(G), *evi=REAL(EVI), *contrast=REAL(CONTRAST), *contrast2=REAL(CONTRAST2);
+ double *el; //=REAL(EL);
+ double *eva, *eve, *evei, *tmp, *P;
+ double *blub=REAL(BLUB), *X;
+ double oldel; //=el[ch-1L]
+ int ancloli, pa, ch; //=anc[loli]
+ double *res = (double *) R_alloc(3L, sizeof(double));
+ tmp = (double *) R_alloc(nr * nc, sizeof(double));
+ P = (double *) R_alloc(nc * nc, sizeof(double));
+ X = (double *) R_alloc(k * nr * nc, sizeof(double));
+
+ SEXP RESULT;
+ PROTECT(RESULT = allocVector(REALSXP, lEL));
+ el=REAL(RESULT);
+ for(i = 0; i < lEL; i++) el[i] = REAL(EL)[i];
+ eva = REAL(VECTOR_ELT(eig, 0));
+ eve = REAL(VECTOR_ELT(eig, 1));
+ evei = REAL(VECTOR_ELT(eig, 2));
+
+ loli = parent[0];
+
+ for(m = 0; m < n; m++){
+ pa = parent[m];
+ ch = child[m];
+ oldel=el[ch-1L];
+
+ while(loli != pa){
+ ancloli=anc[loli];
+ for(i = 0; i < k; i++){
+ getP(eva, eve, evei, nc, el[loli-1L], g[i], P);
+ moveLL5(&LL[LINDEX(loli, i)], &LL[LINDEX(ancloli, i)], P, &nr, &nc, tmp);
+ }
+ loli = ancloli;
+ }
+ // moveDad
+ if(ch>ntips){
+ for(i = 0; i < k; i++){
+ getP(eva, eve, evei, nc, oldel, g[i], P);
+ helpDADI(&LL[LINDEX(pa, i)], &LL[LINDEX(ch, i)], P, nr, nc, tmp);
+ helpPrep(&LL[LINDEX(pa, i)], &LL[LINDEX(ch, i)], eve, evi, nr, nc, tmp, &X[i*nr*nc]);
+ for(h = 0; h < nc; h++){
+ for(j = 0; j < nr; j++){
+ X[j+h*nr + i*nr*nc] *= blub[j+i*nr];
+ }
+ }
+ }
+ }
+ else{
+ for(i = 0; i < k; i++){
+ getP(eva, eve, evei, nc, oldel, g[i], P);
+ helpDAD5(&LL[LINDEX(pa, i)], INTEGER(VECTOR_ELT(dlist, ch-1L)), contrast, P, nr, nc, nco, tmp);
+ helpPrep2(&LL[LINDEX(pa, i)], INTEGER(VECTOR_ELT(dlist, ch-1L)), contrast2, evi, nr, nc, nco, &X[i*nr*nc]); //;
+ for(h = 0; h < nc; h++){
+ for(j = 0; j < nr; j++){
+ X[j+h*nr + i*nr*nc] *= blub[j+i*nr];
+ }
+ }
+ }
+ }
+ fs3(eva, nc, oldel, w, g, X, k, nr, weight, f0, res);
+ updateLL2(dlist, pa, ch, eva, eve, evei, res[0], w, g, nr,
+ nc, ntips, contrast, nco, k, tmp, P);
+ el[ch-1L] = res[0];
+ if (ch > ntips) loli = ch;
+ else loli = pa;
+ }
+ UNPROTECT(1); //RESULT
+ return(RESULT);
+}
+
+
+
+void rowMinScale(int *dat, int n, int k, int *res){
+ int i, h;
+ int tmp;
+ for(i = 0; i < n; i++){
+ tmp = dat[i];
+ for(h = 1; h< k; h++) {if(dat[i + h*n] < tmp) tmp=dat[i + h*n];}
+ if(tmp>0L){for(h = 0; h< k; h++) dat[i + h*n] -= tmp;}
+ res[i] = tmp;
+ }
+}
+
+
+SEXP PML4(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 root, SEXP nco, SEXP contrast, SEXP N){
+ int nr=INTEGER(NR)[0], nc=INTEGER(NC)[0], k=INTEGER(K)[0], i, j, indLL;
+ int nTips = INTEGER(NTips)[0], *SC, *sc;
+ double *g=REAL(G), *w=REAL(W), *tmp, *res;
+ SEXP TMP;
+ double *eva, *eve, *evei;
+ eva = REAL(VECTOR_ELT(eig, 0));
+ eve = REAL(VECTOR_ELT(eig, 1));
+ evei = REAL(VECTOR_ELT(eig, 2));
+ SC = (int *) R_alloc(nr * k, sizeof(int));
+ sc = (int *) R_alloc(nr, sizeof(int));
+ tmp = (double *) R_alloc(nr * k, sizeof(double));
+ PROTECT(TMP = allocVector(REALSXP, nr));
+
+ res=REAL(TMP);
+ for(i=0; i<(k*nr); i++)tmp[i]=0.0;
+ indLL = nr * nc * nTips;
+ for(i=0; i<k; i++){
+ lll3(dlist, eva, eve, evei, REAL(EL), g[i], &nr, &nc, INTEGER(node), INTEGER(edge), nTips, REAL(contrast), INTEGER(nco)[0], INTEGER(N)[0], &SC[nr * i], REAL(bf), &tmp[i*nr], &LL[indLL *i], &SCM[nr * nTips * i]);
+ }
+ rowMinScale(SC, nr, k, sc);
+ for(i=0; i<nr; i++){
+ res[i]=0.0;
+ for(j=0;j<k;j++)res[i] += w[j] * exp(LOG_SCALE_EPS * SC[i+j*nr]) * tmp[i+j*nr];
+ }
+ for(i=0; i<nr; i++) res[i] = log(res[i]) + LOG_SCALE_EPS * sc[i];
+ UNPROTECT(1);
+ return TMP;
+}
+
+
+
+SEXP PML5(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 root, SEXP nco, SEXP contrast, SEXP N){
+ int nr=INTEGER(NR)[0], nc=INTEGER(NC)[0], k=INTEGER(K)[0], i, j, indLL;
+ int nTips = INTEGER(NTips)[0], *SC, *sc;
+ double *g=REAL(G), *w=REAL(W), *tmp, *res;
+ SEXP TMP;
+ double *eva, *eve, *evei;
+ eva = REAL(VECTOR_ELT(eig, 0));
+ eve = REAL(VECTOR_ELT(eig, 1));
+ evei = REAL(VECTOR_ELT(eig, 2));
+ SC = (int *) R_alloc(nr * k, sizeof(int));
+ sc = (int *) R_alloc(nr, sizeof(int));
+ tmp = (double *) R_alloc(nr * k, sizeof(double));
+ PROTECT(TMP = allocVector(REALSXP, nr));
+
+ res=REAL(TMP);
+ for(i=0; i<(k*nr); i++)tmp[i]=0.0;
+ indLL = nr * nc * nTips;
+ for(i=0; i<k; i++){
+ lll(dlist, eva, eve, evei, REAL(EL), g[i], &nr, &nc, INTEGER(node), INTEGER(edge), nTips, REAL(contrast), INTEGER(nco)[0], INTEGER(N)[0], &SC[nr * i], REAL(bf), &tmp[i*nr], &LL[indLL *i]);
+ }
+ rowMinScale(SC, nr, k, sc);
+ for(i=0; i<nr; i++){
+ res[i]=0.0;
+ for(j=0;j<k;j++)res[i] += w[j] * exp(LOG_SCALE_EPS * SC[i+j*nr]) * tmp[i+j*nr];
+ }
+ UNPROTECT(1);
+ return TMP;
+}
+
diff --git a/src/sankoff.c b/src/sankoff.c
index 3a1ff8c..0a22f23 100644
--- a/src/sankoff.c
+++ b/src/sankoff.c
@@ -49,7 +49,7 @@ void rowMin2(double *dat, int n, int k, double *res){
}
-void rowMinInt(int *dat, int n, int k, double *res){
+void rowMinInt(int *dat, int n, int k, int *res){
int i, h;
int x;
for(i = 0; i < n; i++){
diff --git a/tests/testthat.R b/tests/testthat.R
new file mode 100644
index 0000000..8861f50
--- /dev/null
+++ b/tests/testthat.R
@@ -0,0 +1,4 @@
+library(testthat)
+suppressPackageStartupMessages(library(phangorn))
+
+test_check("phangorn")
diff --git a/tests/testthat/test_parsimony.R b/tests/testthat/test_parsimony.R
new file mode 100644
index 0000000..6b30c56
--- /dev/null
+++ b/tests/testthat/test_parsimony.R
@@ -0,0 +1,21 @@
+context("parsimony")
+
+tree1 = read.tree(text = "((t1,t2),t3,t4);")
+tree2 = read.tree(text = "((t1,t3),t2,t4);")
+dat <- phyDat(c(t1="a", t2="a",t3="t",t4="t"), type="USER", levels=c("a","c","g","t"))
+#tr_acctran = acctran(tree1, dat)
+#tr_ratchet = pratchet(dat, trace=0)
+#bab(dat)
+test_that("parsimony works properly", {
+ skip_on_cran()
+ expect_that(fitch(tree1, dat), equals(1))
+ expect_that(fitch(tree2, dat), equals(2))
+ expect_that(sankoff(tree1, dat), equals(1))
+ expect_that(sankoff(tree2, dat), equals(2))
+ expect_that(parsimony(tree1, dat), equals(1))
+})
+
+
+
+
+
diff --git a/tests/testthat/test_phyDat.R b/tests/testthat/test_phyDat.R
new file mode 100644
index 0000000..820165a
--- /dev/null
+++ b/tests/testthat/test_phyDat.R
@@ -0,0 +1,19 @@
+context("conversion_and_subsetting")
+
+data(Laurasiatherian)
+phy_matrix <- as.character(Laurasiatherian)
+phy_df <- as.data.frame(Laurasiatherian)
+phy_dnabin <- as.DNAbin(Laurasiatherian)
+phy_align <- phyDat2alignment(Laurasiatherian)
+
+test_that("conversion 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
diff --git a/tests/testthat/test_treedist.R b/tests/testthat/test_treedist.R
new file mode 100644
index 0000000..62d1b88
--- /dev/null
+++ b/tests/testthat/test_treedist.R
@@ -0,0 +1,14 @@
+context("treedist")
+
+test_that("Robinson-Foulds distance", {
+ skip_on_cran()
+ tree1 = read.tree(text="(t5:1.0,(t4:1.0,t3:1.0):1.0,(t1:1.0,t2:1.0):1.0);")
+ tree2 = read.tree(text="(t4:1.0,(t5:1.0,t3:1.0):1.0,(t1:1.0,t2:1.0):1.0);")
+ tree3 = read.tree(text="(t5:1.0,t4:1.0,t3:1.0,(t1:1.0,t2:1.0):1.0);")
+ expect_that(RF.dist(tree1, tree1), is_equivalent_to(0))
+ expect_that(RF.dist(tree1, tree2), is_equivalent_to(2))
+ expect_that(RF.dist(tree1, tree3), is_equivalent_to(1))
+ expect_that(treedist(tree1, tree1)[1], is_equivalent_to(0))
+ expect_that(treedist(tree1, tree2)[1], is_equivalent_to(2))
+ expect_that(treedist(tree1, tree3)[1], is_equivalent_to(1))
+})
\ No newline at end of file
diff --git a/vignettes/Networx.Rmd b/vignettes/Networx.Rmd
index 8acfaf9..4086c5a 100644
--- a/vignettes/Networx.Rmd
+++ b/vignettes/Networx.Rmd
@@ -11,7 +11,7 @@ vignette: >
---
-This tutorial gives a basic introduction on constructing phylogenetic networks and to add parameter to trees or networx using [phangorn](http://cran.r-project.org/web/packages/phangorn/) [@Schliep2011] in R.
+This tutorial gives a basic introduction on constructing phylogenetic networks and to add parameter to trees or networx using [phangorn](http://cran.r-project.org/package=phangorn) [@Schliep2011] in R.
Splits graph or phylogenetic networks are a nice way to display conflict data or summarize different trees. Here we present to popular networks, consensus networks [@Holland2004]
and neighborNet [@Bryant2004].
Often trees or networks are missing either edge weights or support values about the edges. We show how to improve a tree/networx by adding support values or estimating edge weights using non-negative Least-Squares (nnls).
diff --git a/vignettes/Trees.Rnw b/vignettes/Trees.Rnw
index 6af3b5a..8fba8fb 100644
--- a/vignettes/Trees.Rnw
+++ b/vignettes/Trees.Rnw
@@ -73,7 +73,7 @@ After reading in the alignment we can build a first tree with distance based met
After constructing a distance matrix we reconstruct a rooted tree with UPGMA and alternatively an unrooted tree using Neighbor Joining \cite{Saitou1987,Studier1988}.
<<echo=TRUE>>=
-dm = dist.dna(as.DNAbin(primates))
+dm = dist.ml(primates)
treeUPGMA = upgma(dm)
treeNJ = NJ(dm)
@
diff --git a/vignettes/phangorn.bib b/vignettes/phangorn.bib
index 2c958ef..c6167e3 100644
--- a/vignettes/phangorn.bib
+++ b/vignettes/phangorn.bib
@@ -1,3 +1,142 @@
+ at Manual{CRAN,
+ title = {R: A Language and Environment for Statistical Computing},
+ author = {{R Core Team}},
+ organization = {R Foundation for Statistical Computing},
+ address = {Vienna, Austria},
+ year = {2015},
+ url = {http://www.R-project.org/},
+}
+
+ at Manual{Matrix,
+ title = {Matrix: Sparse and Dense Matrix Classes and Methods},
+ author = {Douglas Bates and Martin Maechler},
+ year = {2015},
+ note = {R package version 1.2-0},
+ url = {http://CRAN.R-project.org/package=Matrix},
+}
+
+ at Manual{seqLogo,
+ title = {seqLogo: Sequence logos for DNA sequence alignments},
+ author = {Oliver Bembom},
+ note = {R package version 1.34.0},
+ year = {2015}
+}
+
+ at Manual{apex,
+ title = {apex: Phylogenetic Methods for Multiple Gene Data},
+ author = {Thibaut Jombart and Zhian Namir Kamvar and Klaus Schliep and Rebecca Harris},
+ note = {R package version 1.0.1},
+ year = {2015}
+}
+
+ at Manual{Turlach2013,
+ title = {quadprog: Functions to solve Quadratic Programming Problems.},
+ author = {Berwin A. Turlach and Andreas Weingessel},
+ year = {2013},
+ note = {R package version 1.5-5},
+ url = {http://CRAN.R-project.org/package=quadprog},
+}
+
+
+ at Manual{rgl,
+ title = {rgl: 3D Visualization Using OpenGL},
+ author = {Daniel Adler and Duncan Murdoch and {others}},
+ year = {2015},
+ note = {R package version 0.95.1247},
+ url = {http://CRAN.R-project.org/package=rgl},
+}
+
+
+ at article{Bouckaert2010,
+ author = {Bouckaert, Remco R.},
+ title = {{DensiTree}: making sense of sets of phylogenetic trees},
+ volume = {26},
+ number = {10},
+ pages = {1372-1373},
+ year = {2010},
+ doi = {10.1093/bioinformatics/btq110},
+ URL = {http://bioinformatics.oxfordjournals.org/content/26/10/1372.abstract},
+ eprint = {http://bioinformatics.oxfordjournals.org/content/26/10/1372.full.pdf+html},
+ journal = {Bioinformatics}
+}
+
+ at article{Cavalli1967,
+ author = {Cavalli-Sforza, L.L. and Edwards, A.W.F.},
+ year = {1967},
+ title = {Phylogenetic analysis: models and estimation procedures},
+ journal = {American Journal of Human Genetics},
+ volume = {19},
+ pages = {233--257}
+}
+
+ at Article{Huson2006,
+ Author = "D.H. Huson and D. Bryant",
+ Title = "Application of Phylogenetic Networks in Evolutionary Studies",
+ Journal = "Molecular Biology and Evolution",
+ Volume = "23",
+ Number = "2",
+ Pages = "254--267",
+ Year = "2006"
+}
+
+ at InCollection{Buneman1971,
+ author = "Peter Buneman",
+ title = "The recovery of trees from measures of dissimilarity",
+ booktitle = "Mathematics in the Archaeological and Historical Sciences",
+ editor = "Hodson, F. R. and Kendall, D. G. and Tautu, P. T.",
+ publisher = "Edinburgh University Press",
+ pages = "387--395",
+ year = "1971"
+}
+
+ at Article{Buneman1974,
+ title = {A Note on the Metric Properties of Trees},
+ author = {Peter Buneman},
+ journal = {Journal of combinatorial theory (B)},
+ year = {1974},
+ volume = {17},
+ pages = {48--50}
+}
+
+ at article{Fitch1971,
+ author = {Fitch, Walter M.},
+ title = {Toward Defining the Course of Evolution: Minimum Change for a Specific Tree Topology},
+ volume = {20},
+ number = {4},
+ pages = {406-416},
+ year = {1971},
+ doi = {10.1093/sysbio/20.4.406},
+ URL = {http://sysbio.oxfordjournals.org/content/20/4/406.abstract},
+ eprint = {http://sysbio.oxfordjournals.org/content/20/4/406.full.pdf+html},
+ journal = {Systematic Biology}
+}
+
+ at article{Lento1995,
+ author = {Lento, G M and Hickson, R E and Chambers, G K and Penny, D},
+ title = {Use of spectral analysis to test hypotheses on the origin of pinnipeds.},
+ volume = {12},
+ number = {1},
+ pages = {28-52},
+ year = {1995},
+ doi = {10.1093/oxfordjournals.molbev.a040189},
+ URL = {http://mbe.oxfordjournals.org/content/12/1/28.abstract},
+ eprint = {http://mbe.oxfordjournals.org/content/12/1/28.full.pdf+html},
+ journal = {Molecular Biology and Evolution}
+}
+
+ at article{Lewis2001,
+ author = {Lewis, Paul O.},
+ title = {A Likelihood Approach to Estimating Phylogeny from Discrete Morphological Character Data},
+ volume = {50},
+ number = {6},
+ pages = {913-925},
+ year = {2001},
+ doi = {10.1080/106351501753462876},
+ URL = {http://sysbio.oxfordjournals.org/content/50/6/913.abstract},
+ eprint = {http://sysbio.oxfordjournals.org/content/50/6/913.full.pdf+html},
+ journal = {Systematic Biology}
+}
+
@Article{ Rambaut1997,
Author = "A. Rambaut and N.C. Grassly",
Title = "Seq-Gen: an application for the Monte Carlo simulation of DNA sequence evolution along phylogenetic trees",
@@ -19,7 +158,7 @@
@Article{Revell2014,
- Title = {Rphylip: an R interface for PHYLIP},
+ Title = {Rphylip: an {R} interface for {PHYLIP}},
Author = {Revell, Liam J. and Chamberlain, Scott A.},
Journal = {Methods in Ecology and Evolution},
Year = {2014},
@@ -47,7 +186,7 @@
@article{Bryant2004,
author = {Bryant, David and Moulton, Vincent},
- title = {Neighbor-Net: An Agglomerative Method for the Construction of Phylogenetic Networks},
+ title = {{Neighbor-Net}: An Agglomerative Method for the Construction of Phylogenetic Networks},
volume = {21},
number = {2},
pages = {255-265},
@@ -214,6 +353,33 @@
Year = "1981"
}
+ at Article{ Felsenstein1985,
+ Author = "Joseph Felsenstein",
+ Title = "Confidence limits on phylogenies. An approach using the bootstrap",
+ Journal = "Evolution",
+ Volume = "39",
+ Pages = "783--791",
+ Year = "1985"
+}
+
+ at Article{Penny1985,
+ Author = "Penny, D. and Hendy, M.D.",
+ Title = "Testing methods evolutionary tree construction",
+ Journal = "Cladistics",
+ Volume = "1",
+ Pages = "266--278",
+ Year = "1985"
+}
+
+ at Article{Penny1986,
+ Author = "Penny, D. and Hendy, M.D.",
+ Title = "Estimating the reliability of evolutionary trees",
+ Journal = "Molecular Biology and Evolution",
+ Volume = "3",
+ Pages = "403--417",
+ Year = "1986"
+}
+
@Book{ Yang2006,
Author = "Ziheng Yang",
Title = "Computational Molecular evolution",
@@ -224,7 +390,7 @@
@Article{ Paradis2004,
Author = "E. Paradis and J. Claude and K. Strimmer",
- Title = "APE: Analyses of Phylogenetics and Evolution in R language",
+ Title = "{APE}: Analyses of Phylogenetics and Evolution in {R} language",
Journal = "Bioinformatics",
Volume = "20",
Number = "2",
@@ -244,7 +410,7 @@
@Book{ Paradis2012,
Author = "Emmanuel Paradis",
Title = "Analysis of Phylogenetics and Evolution with R",
- Edition = "Second",
+ Edition = "Second",
Publisher = "Springer",
Address = "New York",
Year = "2012"
@@ -260,7 +426,7 @@
series = "Biological and Medical Physics, Biomedical Engineering",
pages = "207--232",
address = "New York",
- publisher = "Springer Verlag",
+ publisher = "Springer",
note = "{ISBN :} 978-3-540-35305-8"
}
@@ -296,7 +462,7 @@
@article{Lanfear2012,
author = {Lanfear, Robert and Calcott, Brett and Ho, Simon Y. W. and Guindon, Stephane},
- title = {PartitionFinder: Combined Selection of Partitioning Schemes and Substitution Models for Phylogenetic Analyses},
+ title = {Partition{F}inder: Combined Selection of Partitioning Schemes and Substitution Models for Phylogenetic Analyses},
volume = {29},
number = {6},
pages = {1695-1701},
@@ -307,6 +473,16 @@
journal = {Molecular Biology and Evolution}
}
+ at Article{ Rokas2003,
+ title = "Genome-scale approaches to resolving incongruence in molecular phylogenies",
+ author = "Rokas, A. and Williams, B. L. and King, N., and Carroll, S. B.",
+ journal = "Nature",
+ year = "2011",
+ volume = "425",
+ number = "6960",
+ pages = "798--804"
+}
+
@Article{ Schliep2011b,
title = "Harvesting Evolutionary Signals in a Forest of Prokaryotic Gene Trees",
author = "Klaus Schliep and Philippe Lopez and Fran\c{c}ois-Joseph Lapointe and Eric Bapteste",
@@ -328,7 +504,7 @@
@article{Posada1998,
author = {Posada, D. and Crandall, K.A.},
- title = {MODELTEST: testing the model of DNA substitution.},
+ title = {{MODELTEST}: testing the model of {DNA} substitution.},
volume = {14},
number = {9},
pages = {817-818},
@@ -338,7 +514,7 @@
@article{Posada2008,
author = {Posada, David},
- title = {jModelTest: Phylogenetic Model Averaging},
+ title = {{jModelTest}: Phylogenetic Model Averaging},
volume = {25},
number = {7},
pages = {1253-1256},
@@ -364,7 +540,14 @@
Year = "1973"
}
-
+ at Book{ Burnham2002,
+ Author = "Burnham, K. P. and Anderson, D. R",
+ Title = "Model selection and multimodel inference: a practical information-theoretic approach",
+ Edition = "Second",
+ Publisher = "Springer",
+ Address = "New York",
+ Year = "2002"
+}
@Article{ Lapointe2010,
title = "Clanistics: a multi-level perspective for harvesting unrooted gene trees",
@@ -377,7 +560,6 @@
author = "Fran\c{c}ois-Joseph Lapointe and Philippe Lopez and Yan Boucher and Jeremy Koenig and Eric Bapteste"
}
-
@article{Wilkinson2007,
title = "Of clades and clans: terms for phylogenetic relationships in unrooted trees",
journal = "Trends in Ecology and Evolution",
@@ -389,7 +571,6 @@
author = "Mark Wilkinson and James O. McInerney and Robert P. Hirt and Peter G. Foster and T. Martin Embley"
}
-
@article{Warnow2012,
title = "Standard maximum likelihood analyses of alignments with gaps can be statistically inconsistent",
journal = "PLOS Currents Tree of Life",
--
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