some news for ape 3.0-8

[ape.git] / R / pic.R

diff --git a/R/pic.R b/R/pic.R
index b1cb20b17921668f98e4375a95093d4ad3fda44b..2735f8f24ca82ade9e8a78c6a68a80a73906b56c 100644 (file)
--- a/R/pic.R
+++ b/R/pic.R
@@ -1,26 +1,24 @@
-## pic.R (2012-09-11)
+## pic.R (2013-02-18)
 
 ##   Phylogenetically Independent Contrasts
 
-## Copyright 2002-2012 Emmanuel Paradis
+## Copyright 2002-2013 Emmanuel Paradis
 
 ## This file is part of the R-package `ape'.
 ## See the file ../COPYING for licensing issues.
 
-pic <- function(x, phy, scaled = TRUE, var.contrasts = FALSE)
+pic <- function(x, phy, scaled = TRUE, var.contrasts = FALSE, rescaled.tree = FALSE)
 {
-    if (!inherits(phy, "phylo"))
-      stop("object 'phy' is not of class \"phylo\"")
-    if (is.null(phy$edge.length))
-      stop("your tree has no branch lengths")
+    if (!inherits(phy, "phylo")) stop("object 'phy' is not of class \"phylo\"")
+    if (is.null(phy$edge.length)) stop("your tree has no branch lengths")
     nb.tip <- length(phy$tip.label)
     nb.node <- phy$Nnode
     if (nb.node != nb.tip - 1)
-      stop("'phy' is not rooted and fully dichotomous")
+        stop("'phy' is not rooted and fully dichotomous")
     if (length(x) != nb.tip)
-      stop("length of phenotypic and of phylogenetic data do not match")
+        stop("length of phenotypic and of phylogenetic data do not match")
     if (any(is.na(x)))
-      stop("missing data in 'x': you may consider removing the species with missing data from your tree with the function 'drop.tip'.")
+        stop("missing data in 'x': you may consider removing the species with missing data from your tree with the function 'drop.tip'.")
 
     phy <- reorder(phy, "postorder")
     phenotype <- numeric(nb.tip + nb.node)
@@ -35,10 +33,10 @@ pic <- function(x, phy, scaled = TRUE, var.contrasts = FALSE)
             warning("the names of argument 'x' and the tip labels of the tree did not match: the former were ignored in the analysis.")
         }
     }
+
     ## No need to copy the branch lengths: they are rescaled
     ## in the C code, so it's important to leave the default
     ## `DUP = TRUE' of .C.
-
     ans <- .C("pic", as.integer(nb.tip), as.integer(nb.node),
               as.integer(phy$edge[, 1]), as.integer(phy$edge[, 2]),
               as.double(phy$edge.length), as.double(phenotype),
@@ -46,22 +44,6 @@ pic <- function(x, phy, scaled = TRUE, var.contrasts = FALSE)
               as.integer(var.contrasts), as.integer(scaled),
               PACKAGE = "ape")
 
-    ## The "old" R code:
-    ##for (i in seq(from = 1, by = 2, length.out = nb.node)) {
-    ##    j <- i + 1
-    ##    anc <- phy$edge[i, 1]
-    ##    des1 <- phy$edge[i, 2]
-    ##    des2 <- phy$edge[j, 2]
-    ##    sumbl <- bl[i] + bl[j]
-    ##    ic <- anc - nb.tip
-    ##    contr[ic] <- phenotype[des1] - phenotype[des2]
-    ##    if (scaled) contr[ic] <- contr[ic]/sqrt(sumbl)
-    ##    if (var.contrasts) var.con[ic] <- sumbl
-    ##    phenotype[anc] <- (phenotype[des1]*bl[j] + phenotype[des2]*bl[i])/sumbl
-    ##    k <- which(phy$edge[, 2] == anc)
-    ##    bl[k] <- bl[k] + bl[i]*bl[j]/sumbl
-    ##
-    ##}
     contr <- ans[[7]]
     lbls <-
         if (is.null(phy$node.label)) as.character(1:nb.node + nb.tip)
@@ -70,6 +52,10 @@ pic <- function(x, phy, scaled = TRUE, var.contrasts = FALSE)
         contr <- cbind(contr, ans[[8]])
         dimnames(contr) <- list(lbls, c("contrasts", "variance"))
     } else names(contr) <- lbls
+    if (rescaled.tree) {
+        phy$edge.length <- ans[[5]]
+        contr <- list(contr = contr, rescaled.tree = phy)
+    }
     contr
 }
 
@@ -77,7 +63,7 @@ pic.ortho <- function(x, phy, var.contrasts = FALSE, intra = FALSE)
 {
     n <- length(x)
     m <- n - 1L # number of nodes
-    phy <- reorder(phy, "pruningwise")
+    phy <- reorder(phy, "postorder")
     xx <- unlist(lapply(x, mean)) # 'x' in Felsenstein's paper
     xx <- c(xx, numeric(m))
     delta.v <- numeric(n + m)