[Ecopd-commits] r47 - in pkg: . R man

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

Author: regetz
Date: 2009-10-24 02:01:55 +0200 (Sat, 24 Oct 2009)
New Revision: 47

Added:
   pkg/man/pairdist.Rd
Modified:
   pkg/NAMESPACE
   pkg/R/utilities.R
Log:
wrote pairdist function to compute pairwise node distances for a phylo4
object; added basic documentation


Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE	2009-10-23 16:34:50 UTC (rev 46)
+++ pkg/NAMESPACE	2009-10-24 00:01:55 UTC (rev 47)
@@ -1,4 +1,5 @@
 export(tipLength)
+export(pairdist)
 export(abundance)
 export(`abundance<-`)
 export(minTL)

Modified: pkg/R/utilities.R
===================================================================
--- pkg/R/utilities.R	2009-10-23 16:34:50 UTC (rev 46)
+++ pkg/R/utilities.R	2009-10-24 00:01:55 UTC (rev 47)
@@ -141,3 +141,67 @@
   gsub("_.*$", "", tipLabels(phy))
 }
 
+
+## this works as implementation of dist.nodes for phylo4 objects, albeit
+## about 1.5x slower than dist.nodes
+pairdist <- function(phy, type=c("all", "tip"), use.labels=FALSE) {
+
+  if (hasPoly(phy) || hasRetic(phy)) {
+    stop("pairdist can't currently handle polytomies or reticulation")
+  }
+
+  type <- match.arg(type)
+  edge <- edges(phy, drop.root=TRUE)
+  elen <- edgeLength(phy)
+  elen <- elen[!names(elen) %in% edgeId(phy, "root")]
+  nodes <- nodeId(phy, "all")
+  n <- length(nodes)
+  d <- matrix(NA_real_, nrow=n, ncol=n)
+  diag(d) <- 0
+  d[edge] <- d[edge[,2:1]] <- elen
+  ntips <- nTips(phy)
+  tips <- nodeId(phy, "tip")
+  tip.parents <- sapply(tips, function(n) edge[edge[,2]==n, 1])
+  nodes.todo <- tip.parents[duplicated(tip.parents)]
+  done <- tips
+  root <- if (isRooted(phy)) rootNode(phy) else ntips+1
+  descendants <- matrix(edge[order(edge[,1]),2], nrow=2)
+  ancestors <- edge[match(nodes, edge[,2]), 1]
+
+  while(length(nodes.todo)>0) {
+    nod <- nodes.todo[1]
+    if (nod==root && length(nodes.todo)>1) {
+      nod <- nodes.todo[2]
+    }
+    des1 <- descendants[1, nod-ntips]
+    des2 <- descendants[2, nod-ntips]
+    if (des1>ntips) des1 <- which(!is.na(d[des1, ]))
+    if (des2>ntips) des2 <- which(!is.na(d[des2, ]))
+    for (y in des1) {
+      d[y, des2] <- d[des2, y] <- d[nod, y] + d[nod, des2]
+    }
+    if (nod!=root) {
+      anc <- ancestors[nod]
+      l <- elen[paste(anc, nod, sep="-")]
+      d[des2, anc] <- d[anc, des2] <- d[nod, des2] + l
+      d[des1, anc] <- d[anc, des1] <- d[nod, des1] + l
+      done <- c(done, nod)
+      if (all(descendants[, anc-ntips] %in% done)) {
+        nodes.todo <- c(nodes.todo, anc)
+      }
+    }
+    nodes.todo <- nodes.todo[nod!=nodes.todo]
+  }
+
+  if (type=="tip") {
+    d <- d[1:ntips, 1:ntips]
+  }
+
+  if (use.labels) {
+    rownames(d) <- colnames(d) <- unname(labels(phy, type=type))
+  } else {
+    rownames(d) <- colnames(d) <- nodeId(phy, type=type)
+  }
+
+  return(d)
+}

Added: pkg/man/pairdist.Rd
===================================================================
--- pkg/man/pairdist.Rd	                        (rev 0)
+++ pkg/man/pairdist.Rd	2009-10-24 00:01:55 UTC (rev 47)
@@ -0,0 +1,39 @@
+\name{pairdist}
+\alias{pairdist}
+\title{Calculate pairwise distances between nodes}
+\description{
+  Returns a symmetric matrix containing the pairwise distances between
+  all nodes, or possibly the pairwise distances between tips only.
+}
+\usage{
+  pairdist(phy, type=c("all", "tip"), use.labels=FALSE)
+}
+\arguments{
+  \item{phy}{A \code{phylo4} object (or one that inherits from it).}
+  \item{type}{Should pairwise distances be returned for "all" nodes
+    (default), or only for "tip" nodes?}
+  \item{use.labels}{(logical) If FALSE, row and column names are the
+    numeric node IDs, otherwise labels will be used (possibly resulting
+    in \code{NA} names in the case of missing node labels).}
+}
+\details{
+  TODO
+}
+\value{
+  Numeric matrix of distances, with either node IDs or node labels as
+  row and column names.
+}
+\author{Jim Regetz (regetz at nceas.ucsb.edu)}
+\seealso{
+  ape package functions \code{cophenetic} and \code{dist.nodes}, on
+  which \code{pairdist} is based.
+}
+\examples{
+  data(weeds)
+
+  # all pairwise distances
+  pairdist(weeds$A)
+
+  # pairwise distances between tips, with labels and dim names
+  pairdist(weeds$A, type="tip", use.labels=TRUE)
+}


Property changes on: pkg/man/pairdist.Rd
___________________________________________________________________
Name: svn:eol-style
   + native