[Adephylo-commits] r27 - pkg/R

[noreply at r-forge.r-project.org]

Author: jombart
Date: 2008-11-25 17:45:31 +0100 (Tue, 25 Nov 2008)
New Revision: 27

Modified:
   pkg/R/distances.R
   pkg/R/utils.R
Log:
Added a new argument to sp.tips, include.mrca (faster when FALSE).


Modified: pkg/R/distances.R
===================================================================
--- pkg/R/distances.R	2008-11-25 15:23:37 UTC (rev 26)
+++ pkg/R/distances.R	2008-11-25 16:45:31 UTC (rev 27)
@@ -1,50 +1,70 @@
 ############
 # distNodes
 ############
-distNodes <- function(x, node1, node2,
+distTips <- function(x, tips="all",
                       method=c("brlength","nNodes","Abouheif","sumDD")){
 
     if(!require(phylobase)) stop("phylobase package is not installed")
 
-    ## conversion from phylo, phylo4 and phylo4d
+    ## handle arguments
     x <- as(x, "phylo4")
     method <- match.arg(method)
+    tips <- getnodes(x, tip1)
+    N <- nTips(x)
+    if(tips="all") { tips <- 1:N }
 
     ## some checks
     if (is.character(checkval <- check_phylo4(x))) stop(checkval)
-    node1 <- getnodes(x, node1)
-    node2 <- getnodes(x, node2)
-    if(any(is.na(c(node1,node2)))) stop("wrong node specified")
-    if(node1==node2) return(0)
+    if(any(is.na(tips))) stop("wrong tips specified")
 
-    ## get the path between node1 and node2
-    path <- shortestPath(x, node1, node2)
+    ## create all couples of observations
+    findAllPairs <- function(vec){
+        res <- list(i=NULL,j=NULL)
+        k <- 0
+        for(i in 1:(length(vec)-1)){
+            for(j in 2:length(vec)){
+                k <- k+1
+                res[[1]][k] <- i
+                res[[2]][k] <- j
+            }
+        }
+        res <- data.frame(res)
+        return(res)
+    }
 
+    allPairs <- findAllPairs(tips) # this contains all possible pairs of tips
+
+    ## get the shortest path between all pairs of tips
+    allPath <- sp.tips(x, allPairs$i, allPairs$j, useTipNames=TRUE, quiet=TRUE)
+
     ## compute distances
     if(method=="brlength"){
         if(!hasEdgeLength(x)) stop("x does not have branch length")
-        path <- c(node1, node2, path)
-        path <- path[path != MRCA(x, node1, node2)]
-        edge.idx <- getedges(x, path)
+        ## add tip1 and tip2 to the paths, so that these edges are counted
+        allPath <- as.data.frame
+        
+        allPath <- lapply(allPath, c, (node1, node2, allPath)
+        allPath <- allPath[allPath != MRCA(x, node1, node2)]
+        edge.idx <- getedges(x, allPath)
         res <- sum(edgeLength(x)[edge.idx], na.rm=TRUE)
         return(res)
     } # end brlength
 
     if(method=="nNodes"){
-        res <- length(path)
+        res <- length(allPath)
         return(res)
     } # end nNodes
 
     if(method=="Abouheif"){
         E <- x at edge
-        temp <- table(E[,1])[as.character(path)] # number of dd per node
+        temp <- table(E[,1])[as.character(allPath)] # number of dd per node
         res <- prod(temp)
         return(res)
     } # end Abouheif
 
     if(method=="sumDD"){
         E <- x at edge
-        temp <- table(E[,1])[as.character(path)] # number of dd per node
+        temp <- table(E[,1])[as.character(allPath)] # number of dd per node
         res <- sum(temp)
         return(res)
  } # end sumDD

Modified: pkg/R/utils.R
===================================================================
--- pkg/R/utils.R	2008-11-25 15:23:37 UTC (rev 26)
+++ pkg/R/utils.R	2008-11-25 16:45:31 UTC (rev 27)
@@ -1,7 +1,7 @@
 ##########
 # sp.tips
 ##########
-sp.tips <- function(x, tip1, tip2, useTipNames=FALSE, quiet=FALSE){
+sp.tips <- function(x, tip1, tip2, useTipNames=FALSE, quiet=FALSE, include.mrca=TRUE){
     if(!require(phylobase)) stop("phylobase package is not installed")
 
     ## conversion from phylo, phylo4 and phylo4d
@@ -62,6 +62,15 @@
     } # end pathTwoTips
 
 
+    pathTwoTips.no.mrca <- function(path1, path2){
+        cpath <- c(path1, path2)
+        temp <- intersect(path1, path2)
+        res <- setdiff(cpath, temp)
+        return(res)
+    } # end pathTwoTips
+
+
+
     ## main computations
     allPathToRoot <- lapply(allTips, function(i) tipToRoot(E, i))
     names(allPathToRoot) <- allTips
@@ -69,7 +78,15 @@
     allPath1 <- allPathToRoot[as.character(t1)]
     allPath2 <- allPathToRoot[as.character(t2)]
 
-    res <- lapply(1:length(allPath1), function(i) pathTwoTips(allPath1[[i]], allPath2[[i]]) )
+    if(include.mrca) {
+        res <- lapply(1:length(allPath1), function(i) pathTwoTips(allPath1[[i]], allPath2[[i]]) )
+    } else {
+        res <- lapply(1:length(allPath1), function(i) pathTwoTips.no.mrca(allPath1[[i]], allPath2[[i]]) )
+        temp.names <- names(res)
+        temp <- sapply(res, function(vec) length(vec)>0)
+        res[temp] <- lapply(res[temp], function(vec) getnodes(x, vec) ) # name the nodes
+        names(res) <- temp.names
+    }
 
     if(useTipNames) {
         names(res) <- paste(names(t1), names(t2), sep="-")
@@ -84,15 +101,15 @@
 
 # examples
 # source("/home/master/dev/adephylo/pkg/R/utils.R")
-## phy <- as(rtree(15),"phylo4")
-## plot(phy,show.n=T)
-## tip1 <- "t1"
-## tip2 <- "t2"
+#phy <- as(rtree(15),"phylo4")
+plot(phy,show.n=T)
+tip1 <- "t1"
+tip2 <- "t2"
 
 
-## sp.tips(phy, "t1", "t2")
-## sp.tips(phy, rep(1,15), 1:15)
-## sp.tips(phy, rep(1, 15), 1:15, TRUE)
+sp.tips(phy, "t1", "t2")
+sp.tips(phy, rep(1,15), 1:15)
+sp.tips(phy, rep(1, 15), 1:15, TRUE)
 
 ## heavier tree
 # x <- as(rtree(1000), "phylo4")