[Phylobase-commits] r640 - in pkg: R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Sep 12 03:01:43 CEST 2009
Author: pdc
Date: 2009-09-12 03:01:42 +0200 (Sat, 12 Sep 2009)
New Revision: 640
Added:
pkg/src/phyloXX.c
Modified:
pkg/R/treePlot.R
Log:
Speed up plotting. Plotting of simple trees now faster than ape for most sane cases (read < 10,000 tips).
The function phyloXXYY() that calculates the X and Y positions of the line segments has been updated to use a C function (phyloXX) to calculate the X locations more speedily.
I promise this is the last unrequested optimization to the plotting code.
Modified: pkg/R/treePlot.R
===================================================================
--- pkg/R/treePlot.R 2009-09-12 00:47:39 UTC (rev 639)
+++ pkg/R/treePlot.R 2009-09-12 01:01:42 UTC (rev 640)
@@ -230,12 +230,30 @@
## Set root x value to zero and calculate x positions
xx[1] <- 0
segs$v0x[1] <- segs$v1x[1] <- segs$h0x[1] <- 0
- for (i in edge[, 2]) {
- dex <- edge[, 1] == i
- cur <- edge[, 2] == i
- xx[dex] <- phy at edge.length[dex] + xx[cur]
- segs$h1x[cur] <- segs$v0x[dex] <- segs$v1x[dex] <- segs$h0x[dex] <- xx[cur]
- }
+ edge1 <- as.integer(edge[,1])
+ edge2 <- as.integer(edge[,2])
+ edgeLen <- phy at edge.length
+ edgeLen[is.na(edgeLen)] <- 0
+ edgeLen <- as.numeric(edgeLen)
+ nedges <- as.integer(nEdges(phy))
+ segsv0x <- as.numeric(rep.int(0, Nedges))
+ xPos <- .C("phyloxx", edge1, edge2,
+ edgeLen, nedges, xx, segsv0x)
+ ## browser()
+ xx <- xPos[[5]]
+ segs$v0x <- xPos[[6]]
+ ## test1 <- function() {
+ ## for (i in edge[, 2]) {
+ ## dex <- edge[, 1] == i
+ ## cur <- edge[, 2] == i
+ ## xx[dex] <- phy at edge.length[dex] + xx[cur]
+ ## segs$v0x[dex] <- xx[cur]
+ ## }
+ ## return(list(segs=segs, xx=xx))
+ ## }
+ ## test1out <- test1()
+ ## segs <- test1out$segs
+ ## xx <- test1out$xx
## Set y positions for terminal nodes and calculate remaining y positions
if(!is.null(tip.order)) {
@@ -245,20 +263,37 @@
}
segs$h0y[tips] <- segs$h1y[tips] <- yy[tips]
segs$v1y[tips] <- segs$v0y[tips] <- yy[tips]
- for(i in rev((Ntips + 1):nEdges(phy))) {
- dex <- edge[, 1] == i
- cur <- edge[, 2] == i
- yy[cur] <- segs$h0y[cur] <- segs$h1y[cur] <- segs$v1y[cur] <- segs$v0y[dex] <- mean(yy[dex])
+ placeHolder <- function() {
+ for(i in rev((Ntips + 1):nEdges(phy))) {
+ dex <- edge[, 1] == i
+ cur <- edge[, 2] == i
+ yy[cur] <- segs$v0y[dex] <- mean(yy[dex])
+ }
+ return(list(segs=segs, yy=yy))
}
-
+ yPos <- placeHolder()
+ segs <- yPos$segs
+ yy <- yPos$yy
+
+ ## edgeLen[is.na(edgeLen)] <- 0
+ ## edgeLen <- as.numeric(edgeLen)
+ ## ntips <- as.integer(nTips(phy))
+ ## yy <- as.numeric(yy)
+ ## segsv0y <- as.numeric(yy)
+ ## browser()
+ ## yPos <- .C("phyloyy", edge1, edge2,
+ ## ntips, nedges, yy, segsv0y)
+
+ segs$h0y <- segs$h1y <- segs$v1y <- yy
+
## scale the x values so they range from 0 to 1
Xmax <- max(xx)
segs$v0x <- segs$v0x / Xmax
- segs$v1x <- segs$v1x / Xmax
- segs$h0x <- segs$h0x / Xmax
- segs$h1x <- segs$h1x / Xmax
xx <- xx / Xmax
+ segs$h1x <- xx
+ segs$v1x <- segs$h0x <- segs$v0x
+
# TODO return an index vector instead of a second phy object
list(xx = xx, yy = yy, phy = phy, phy.orig = phy.orig, segs = segs)
}
Added: pkg/src/phyloXX.c
===================================================================
--- pkg/src/phyloXX.c (rev 0)
+++ pkg/src/phyloXX.c 2009-09-12 01:01:42 UTC (rev 640)
@@ -0,0 +1,97 @@
+/*
+ descendants.c:
+ Identify all descendants of a given node. Function inputs are
+ derived from a phylo4 edge matrix, which *must* be in preorder order.
+ The isDescendant input vector should contain 1 for the immediate
+ children of the node, and 0 otherwise. The function returns this
+ vector updated to include all further descendants.
+*/
+
+// test1 <- function() {
+// for (i in edge[, 2]) {
+// dex <- edge[, 1] == i
+// cur <- edge[, 2] == i
+// xx[dex] <- phy at edge.length[dex] + xx[cur]
+// segs$v0x[dex] <- xx[cur]
+// }
+// return(list(segs=segs, xx=xx))
+// }
+// test1out <- test1()
+// segs <- test1out$segs
+// xx <- test1out$xx
+
+// test2 <- function() {
+// for(i in rev((Ntips + 1):nEdges(phy))) {
+// dex <- edge[, 1] == i
+// cur <- edge[, 2] == i
+// yy[cur] <- segs$v0y[dex] <- mean(yy[dex])
+// }
+// return(list(segs=segs, yy=yy))
+// }
+// test2out <- test2()
+// segs <- test2out$segs
+// yy <- test2out$yy
+// segs$h0y <- segs$h1y <- segs$v1y <- yy
+
+#include <R.h>
+
+//
+// void phyloyy(int *edge1, int *edge2, int *ntips,
+// int *numEdges, double *yy, double *v0y)
+// {
+// int i;
+// int k;
+// int j;
+// int cur;
+// int des;
+// int count;
+// double tmp;
+// double theMean;
+// Rprintf("test\n");
+// for (i=*numEdges; i > *ntips ; i--) {
+// for (k=0; k<*numEdges; k++) {
+// if(i == edge2[k]) {
+// cur = k;
+// }
+// }
+// tmp=0;
+// count=0;
+// for (j=0; j<*numEdges; j++) {
+// if(i == edge1[j]) {
+// des = j;
+// tmp += yy[j];
+// count += 1;
+// }
+// }
+// theMean = tmp / count;
+// yy[cur] = theMean;
+// for (j=0; j<*numEdges; j++) {
+// if(i == edge1[j]) {
+// v0y[j] = theMean;
+// }
+// }
+//
+// }
+// }
+
+void phyloxx(int *edge1, int *edge2, double *edgeLengths,
+ int *numEdges, double *xx, double *v0x)
+{
+ int j;
+ int i;
+ int k;
+ int cur;
+ for (i=0; i <*numEdges; i++) {
+ for (k=0; k<*numEdges; k++) {
+ if(edge2[i] == edge2[k]) {
+ cur = k;
+ }
+ }
+ for (j=0; j<*numEdges; j++) {
+ if(edge2[i] == edge1[j]) {
+ xx[j] = edgeLengths[j] + xx[cur];
+ v0x[j] = xx[cur];
+ }
+ }
+ }
+}
\ No newline at end of file
More information about the Phylobase-commits
mailing list