[Adephylo-commits] r157 - in pkg: R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Mar 14 20:48:40 CET 2010
Author: jombart
Date: 2010-03-14 20:48:40 +0100 (Sun, 14 Mar 2010)
New Revision: 157
Modified:
pkg/R/distances.R
pkg/src/distPhylo.c
Log:
Distances are now computed in C.
It works !!!
Modified: pkg/R/distances.R
===================================================================
--- pkg/R/distances.R 2010-03-12 17:01:08 UTC (rev 156)
+++ pkg/R/distances.R 2010-03-14 19:48:40 UTC (rev 157)
@@ -2,93 +2,113 @@
# distTips
###########
distTips <- function(x, tips="all",
- method=c("patristic","nNodes","Abouheif","sumDD")){
+ method=c("patristic","nNodes","Abouheif","sumDD"), useC=TRUE){
if(!require(phylobase)) stop("phylobase package is not installed")
- ## handle arguments
- x <- as(x, "phylo4")
- method <- match.arg(method)
- N <- nTips(x)
- if(tips[1]=="all") { tips <- 1:N }
- tips <- getNode(x, tips)
- tips.names <- names(tips)
+ if(useC){
+ tre <- as(x, "phylo")
+ n <- as.integer(nTips(tre))
+ resSize <- as.integer(n*(n-1)/2)
+ res <- double(resSize)
+ method <- match.arg(method)
+ method <- match(method, c("patristic","nNodes","Abouheif","sumDD"))
- ## some checks
- if (is.character(checkval <- checkPhylo4(x))) stop(checkval)
- if(any(is.na(tips))) stop("wrong tips specified")
+ temp <- .C("distalltips", as.integer(tre$edge[,1]), as.integer(tre$edge[,2]), as.double(tre$edge.length), nrow(tre$edge), n, res, resSize, as.integer(method), PACKAGE="adephylo")
+ res <- temp[[6]]
- ## 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 (i+1):length(vec)){
- k <- k+1
- res[[1]][k] <- i
- res[[2]][k] <- j
+ class(res) <- "dist"
+ attr(res, "Size") <- nTips(tre)
+ attr(res, "Diag") <- FALSE
+ attr(res, "Upper") <- FALSE
+ attr(res, "method") <- paste("Phylogenetic: ",method,sep="")
+ attr(res, "call") <- match.call()
+ attr(res, "Labels") <- tre$tip.label
+ } else {
+
+ ## handle arguments
+ x <- as(x, "phylo4")
+ method <- match.arg(method)
+ N <- nTips(x)
+ if(tips[1]=="all") { tips <- 1:N }
+ tips <- getNode(x, tips)
+ tips.names <- names(tips)
+
+ ## some checks
+ if (is.character(checkval <- checkPhylo4(x))) stop(checkval)
+ if(any(is.na(tips))) stop("wrong tips specified")
+
+ ## 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 (i+1):length(vec)){
+ k <- k+1
+ res[[1]][k] <- i
+ res[[2]][k] <- j
+ }
}
+ res <- data.frame(res)
+ return(res)
}
- res <- data.frame(res)
- return(res)
- }
- allPairs <- findAllPairs(tips) # this contains all possible pairs of tips
+ allPairs <- findAllPairs(tips) # this contains all possible pairs of tips
- ## get the shortest path between all pairs of tips
- if(method != "patristic") {
- allPath <- sp.tips(x, allPairs$i, allPairs$j, useTipNames=TRUE, quiet=TRUE)
- } else {
- allPath <- sp.tips(x, allPairs$i, allPairs$j, useTipNames=TRUE, quiet=TRUE,
- include.mrca=FALSE)
- }
+ ## get the shortest path between all pairs of tips
+ if(method != "patristic") {
+ allPath <- sp.tips(x, allPairs$i, allPairs$j, useTipNames=TRUE, quiet=TRUE)
+ } else {
+ allPath <- sp.tips(x, allPairs$i, allPairs$j, useTipNames=TRUE, quiet=TRUE,
+ include.mrca=FALSE)
+ }
- ## compute distances
- if(method=="patristic"){
- if(!hasEdgeLength(x)) stop("x does not have branch length")
- ## add tip1 and tip2 to the paths, so that these edges are counted
- allPath.names <- names(allPath)
- allPath <- lapply(1:length(allPath), function(i)
- c(allPath[[i]], allPairs[i,1], allPairs[i,2]) )
- names(allPath) <- allPath.names
+ ## compute distances
+ if(method=="patristic"){
+ if(!hasEdgeLength(x)) stop("x does not have branch length")
+ ## add tip1 and tip2 to the paths, so that these edges are counted
+ allPath.names <- names(allPath)
+ allPath <- lapply(1:length(allPath), function(i)
+ c(allPath[[i]], allPairs[i,1], allPairs[i,2]) )
+ names(allPath) <- allPath.names
- edge.idx <- lapply(allPath, function(e) getEdge(x, e) ) # list of indices of edges
- allEdgeLength <- edgeLength(x)
- res <- lapply(edge.idx, function(idx) sum(allEdgeLength[idx], na.rm=TRUE) )
- } # end patristic
+ edge.idx <- lapply(allPath, function(e) getEdge(x, e) ) # list of indices of edges
+ allEdgeLength <- edgeLength(x)
+ res <- lapply(edge.idx, function(idx) sum(allEdgeLength[idx], na.rm=TRUE) )
+ } # end patristic
- if(method=="nNodes"){
- res <- lapply(allPath, length)
- } # end nNodes
+ if(method=="nNodes"){
+ res <- lapply(allPath, length)
+ } # end nNodes
- if(method=="Abouheif"){
- E <- x at edge
- f1 <- function(onePath){ # computes product of dd for one path
- temp <- table(E[,1])[as.character(onePath)] # number of dd per node
- return(prod(temp))
- }
- res <- lapply(allPath, f1)
- } # end Abouheif
+ if(method=="Abouheif"){
+ E <- x at edge
+ f1 <- function(onePath){ # computes product of dd for one path
+ temp <- table(E[,1])[as.character(onePath)] # number of dd per node
+ return(prod(temp))
+ }
+ res <- lapply(allPath, f1)
+ } # end Abouheif
- if(method=="sumDD"){
- E <- x at edge
- f1 <- function(onePath){ # computes sum of dd for one path
- temp <- table(E[,1])[as.character(onePath)] # number of dd per node
- return(sum(temp))
- }
- res <- lapply(allPath, f1)
- } # end sumDD
+ if(method=="sumDD"){
+ E <- x at edge
+ f1 <- function(onePath){ # computes sum of dd for one path
+ temp <- table(E[,1])[as.character(onePath)] # number of dd per node
+ return(sum(temp))
+ }
+ res <- lapply(allPath, f1)
+ } # end sumDD
- ## convert res to a dist object
- res <- unlist(res)
- class(res) <- "dist"
- attr(res, "Size") <- length(tips)
- attr(res, "Diag") <- FALSE
- attr(res, "Upper") <- FALSE
- attr(res, "method") <- paste("Phylogenetic: ",method,sep="")
- attr(res, "call") <- match.call()
- attr(res, "Labels") <- tips.names
-
+ ## convert res to a dist object
+ res <- unlist(res)
+ class(res) <- "dist"
+ attr(res, "Size") <- length(tips)
+ attr(res, "Diag") <- FALSE
+ attr(res, "Upper") <- FALSE
+ attr(res, "method") <- paste("Phylogenetic: ",method,sep="")
+ attr(res, "call") <- match.call()
+ attr(res, "Labels") <- tips.names
+ }
return(res)
} # end distTips
Modified: pkg/src/distPhylo.c
===================================================================
--- pkg/src/distPhylo.c 2010-03-12 17:01:08 UTC (rev 156)
+++ pkg/src/distPhylo.c 2010-03-14 19:48:40 UTC (rev 157)
@@ -94,20 +94,40 @@
- for patristic distances, the set of edge used is: {output of sp2tips} U {tipA, tipB} \ {output of mrca2tips}
- for all others: {output of sp2tips}
*/
-int dist2tips(int *ances, int *desc, double *brlength, int N, int tipA, int tipB, int method){
+double dist2tips(int *ances, int *desc, double *brlength, int N, int tipA, int tipB, int method){
/* declarations */
int *path, *lengthPath, *myMrca;
- int i, res;
+ int i;
+ double res;
/* allocate memory */
vecintalloc(&path, N);
lengthPath = (int *) calloc(1, sizeof(int));
- myMrca = (int *) calloc(1, sizeof(int));
+ vecintalloc(&myMrca, 1); /* has to be this way, not a simple pointer, to be used in intANotInB */
+ /* /\* debugging *\/ */
+ /* printf("\n-- Input to dist2tips --\n"); */
+ /* printf("\nances:"); */
+ /* for(i=1;i<=N; i++){ */
+ /* printf(" %d", ances[i]); */
+ /* } */
+ /* printf("\ndesc:"); */
+ /* for(i=1;i<=N; i++){ */
+ /* printf(" %d", desc[i]); */
+ /* } */
+ /* printf("\nedge length:"); */
+ /* for(i=1;i<=N; i++){ */
+ /* printf(" %f", brlength[i]); */
+ /* } */
+
/* find the shortest path between the two tips */
sp2tips(ances, desc, N, tipA, tipB, path, lengthPath);
+ /* printf("\nShortest path found in dist2tips:"); */
+ /* for(i=1;i<=*lengthPath; i++){ */
+ /* printf(" %d", path[i]); */
+ /* } */
/* compute the distance */
@@ -115,10 +135,11 @@
{
case 1: /* patristic distance */
/* find the mrca */
- *myMrca = 0;
- *myMrca = mrca2tips(ances, desc, tipA, tipB, N);
+ myMrca[1] = 0;
+ myMrca[1] = mrca2tips(ances, desc, tipA, tipB, N);
/* remove mrca from the path */
+ /* printf("\nMRCA: %d", myMrca[1]); */
intANotInB(path, myMrca, *lengthPath, 1, path, lengthPath);
/* add tips to the path */
@@ -127,10 +148,17 @@
*lengthPath = *lengthPath + 1;
path[*lengthPath] = tipB;
+ /* printf("\nPath used in patristic distance:"); */
+ /* for(i=1;i<=*lengthPath; i++){ */
+ /* printf(" %d", path[i]); */
+ /* } */
+
/* compute length */
res=0.0;
for(i=1; i<=*lengthPath; i++){
- res += findedgelength(desc, brlength, N, path[i]);
+ res = res + findedgelength(desc, brlength, N, path[i]);
+ /* printf("\nlength of edge terminating by %d: %f", path[i], findedgelength(desc, brlength, N, path[i])); */
+ /* printf("\n value of res: %f", res); */
}
break;
@@ -141,14 +169,14 @@
case 3: /* prod DD (Abouheif) */
res=1.0;
for(i=1; i<=*lengthPath; i++){
- res *= findNbDD(ances, desc, N, path[i]);
+ res = res * findNbDD(ances, desc, N, path[i]);
}
break;
case 4: /* sum DD */
res=0.0;
for(i=1; i<=*lengthPath; i++){
- res += findNbDD(ances, desc, N, path[i]);
+ res = res + findNbDD(ances, desc, N, path[i]);
}
break;
@@ -161,8 +189,9 @@
/* free memory */
freeintvec(path);
free(lengthPath);
- free(myMrca);
+ freeintvec(myMrca);
+ /* printf("\nDistance between %d and %d: %f", tipA, tipB, res); */
return(res);
} /* end dist2tips */
@@ -209,7 +238,7 @@
/* create local vectors for ancestors, descendents and branch lengths */
ancesLoc[0] = *N;
descLoc[0] = *N;
- brlengthLoc[0] = *N ; /* conversion (casting) int->double*/
+ brlengthLoc[0] = *N ; /* implicit casting int->double */
for(i=0; i< *N; i++){
ancesLoc[i+1] = ances[i];
descLoc[i+1] = desc[i];
@@ -218,11 +247,12 @@
/* perform computations for all pairs of tips (indexed 'i,j') */
- k = 0; /* used to browse 'res' and 'resId' */
+ k = 0; /* used to browse 'res' */
for(i=1; i<=(*nTips-1); i++){
for(j=(i+1); j<=(*nTips); j++){
res[k] = dist2tips(ancesLoc, descLoc, brlengthLoc, *N, i, j, *method);
+ /*printf("\nDistance between tip %d and %d in main function: %f", i, j, res[k]);*/
k++;
}
}
@@ -241,25 +271,29 @@
/*
library(adephylo)
-tre=rtree(10)
+tre=rtree(5)
+#tre$edge.length <- round(tre$edge.length*10)
+#tre$edge.length[tre$edge.length<1] <- 1
+
plot(tre)
nodelabels()
tiplabels()
+edgelabels(text=tre$edge.length)
n <- as.integer(nTips(tre))
resSize=as.integer(n*(n-1)/2)
-res <- integer(resSize)
+res <- numeric(resSize)
# void distalltips(int *ances, int *desc, double *brlength, int *N, int *nTips, double *res, int *resSize, int *method){
-## nb nodes
-toto <- .C("distalltips", as.integer(tre$edge[,1]), as.integer(tre$edge[,2]), as.double(tre$edge.length), nrow(tre$edge), n, res, length(res), as.integer(2))
+## patristic
+toto <- .C("distalltips", as.integer(tre$edge[,1]), as.integer(tre$edge[,2]), as.double(tre$edge.length), nrow(tre$edge), n, res, length(res), as.integer(1))
res <- toto[[6]]
res
-## patristic
-toto <- .C("distalltips", as.integer(tre$edge[,1]), as.integer(tre$edge[,2]), as.double(tre$edge.length), nrow(tre$edge), n, res, length(res), as.integer(1))
+## nNodes
+toto <- .C("distalltips", as.integer(tre$edge[,1]), as.integer(tre$edge[,2]), as.double(tre$edge.length), nrow(tre$edge), n, res, length(res), as.integer(2))
res <- toto[[6]]
res
More information about the Adephylo-commits
mailing list