[Picante-commits] r178 - pkg/R

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

Author: skembel
Date: 2009-02-24 00:54:11 +0100 (Tue, 24 Feb 2009)
New Revision: 178

Modified:
   pkg/R/phylosor.R
Log:
Fix phylosor code to deal with tree-community pruning mismatch

Modified: pkg/R/phylosor.R
===================================================================
--- pkg/R/phylosor.R	2009-02-23 22:28:51 UTC (rev 177)
+++ pkg/R/phylosor.R	2009-02-23 23:54:11 UTC (rev 178)
@@ -112,48 +112,41 @@
 .pdshort=function(comm,tree,keep.root=TRUE)
 {
 
+    if (!is.rooted(tree) || !is.ultrametric(tree)) {
+        stop("Rooted ultrametric tree required for phylosor calculation")
+    }
+
 	nbspecies <- length(comm)
 	species <- names(comm)
-    absent <- species[comm==0]	#species not in sample
-    present <- species[comm>0]    #species in sample
-        if(length(present)==0)
+    #fixme prune species properly but keeping track of root if needed
+
+    present <- species[comm>0]  #species in sample
+    commabsent <- species[comm==0]	#species not in sample
+    treeabsent <- tree$tip.label[which(!(tree$tip.label %in% present))]
+    
+    if(length(present)==0)
     {
         #no species present
         PD<-0
     }
     else if(length(present)==1)
     {
-        #one species present - PD = length from root to that tip
-        if (!is.rooted(tree) || !keep.root) {
-            stop("Rooted tree and keep.root=TRUE required to calculate PD of single-species communities")
-        }
-        
+        #one species present - PD = length from root to that tip        
         PD <- node.age(tree)$ages[ which(tree$edge[,2] == 
                                     which(tree$tip.label==present))]
     }
-    else if(length(absent)==0)
+    else if(length(treeabsent)==0)
     {
-        #all species present
+        #all species in tree present in community
         PD <- sum(tree$edge.length)
     }
     else
     {
-        #subset of species present
-        sub.tree<-drop.tip(tree,absent) 
-        if (keep.root) {
-            if (!is.rooted(tree) || !is.ultrametric(tree)) {
-                stop("Rooted ultrametric tree required to calculate PD keeping root node")
-            }
-            #fixme - currently assumes ultrametric tree, could fix this
-            sub.tree.depth <- max(node.age(sub.tree)$ages)
-            orig.tree.depth <- max(node.age(tree)$ages)
-            PD<-sum(sub.tree$edge.length) + (orig.tree.depth-sub.tree.depth)        
-        }
-        else {
-            PD<-sum(sub.tree$edge.length)
-        }
+        #subset of tree species present in community
+        sub.tree<-drop.tip(tree,treeabsent) 
+        sub.tree.depth <- max(node.age(sub.tree)$ages)
+        orig.tree.depth <- max(node.age(tree)$ages)
+        PD<-sum(sub.tree$edge.length) + (orig.tree.depth-sub.tree.depth)
     }
     return(PD)
 }
-
-