[Mattice-commits] r245 - pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Mar 14 17:35:48 CET 2011
Author: andrew_hipp
Date: 2011-03-14 17:35:48 +0100 (Mon, 14 Mar 2011)
New Revision: 245
Removed:
pkg/R/OLDregimes.R
Log:
Moved to MISC
Deleted: pkg/R/OLDregimes.R
===================================================================
--- pkg/R/OLDregimes.R 2011-03-14 16:35:28 UTC (rev 244)
+++ pkg/R/OLDregimes.R 2011-03-14 16:35:48 UTC (rev 245)
@@ -1,167 +0,0 @@
-# ---------------------------------------------------
-# FUNCTIONS FOR PAINTING REGIMES ON AN S4 OUCH TREE #
-# ---------------------------------------------------
-# Adapted from functions used in Hipp 2007 Evolution paper
-# Initially written for ouch v 1.2-4
-# updated to ouch >= 2.4-2 Nov 2008
-# updated to accommodate multiple trees Nov 2008
-
-regimeVectors <-
-# This is the basic call to get the full range of regimes over a set of trees
-# Generates the list of painted branches representing all possible selective regimes for OU analyses, taking as argument
-# species vectors that describe the clades at the bases of which regimes are specified to change.
-# Arguments:
-# "tree" = the standard tree specification vectors of the OUCH-style tree
-# "cladeMembersList" = list of vectors containing names of the members of each clade (except for the root of the tree)
-# "maxNodes" = maximum number of nodes at which regime is permitted to change
-# Value:
-# "regList" = list of vectors that can each be plugged directly into OU analysis as the "regimes" argument
-# "nodeMatrix" = matrix of trees (rows) by nodes (columns) indicating what nodes are present on which trees
-# 19 nov 08: changing to accept a list of trees and trimmed down greatly
-function(ouchTrees, cladeMembersList, maxNodes = NULL) {
- nnode <- length(cladeMembersList)
- regMatrix <- regimeMatrix(n = nnode, maxNodes = maxNodes)
- apr = regimeMaker(ouchTrees, regMatrix, cladeMembersList)
- outdata <- list(regList = apr$regList, regMatrix = apr$regMatrix, nodeMatrix = apr$nodeMatrix)
- return(outdata)
-}
-
-paintBranches <-
-# Paints branches with regimes changing at nodes specified
-# arguments
-# "tree" = OUCH-style (S4) tree
-# "regimeShiftNodes" = a vector of nodes or a list of taxa defining nodes at which selective regimes shift:
-# root may be included--if not present, will be added--but tips are meaningless in this context
-# "regimeTitles" = a vector of titles for the regimes that begin at the root and at the nodes indicated in "regimeShiftNodes",
-# in order of description in "regimeShiftNodes", except that the root is listed first in "regimeTitles"
-# but not at all in "regimeShiftNodes"... defaults to "node[x]regime
-# Value: a vector of regimes that can be handed to hansen
-# this is an especially ugly function, and I'm sure there are prettier ways of doing it. The 'paint' function in
-# ouch could be adapted to this purpose, if one makes a series of consecutive calls going down the tree,
-# and down the road it would probably make sense to turn this into just such a function.
-
-function(regimeShiftNodes, tree, regimeTitles = NULL) {
- ## ------------------ begin ouchtree block -----------------
- ## check to see if tree inherits 'ouchtree'
- if (!is(tree,'ouchtree'))
- stop(paste('This function has been rewritten to use the new S4 ', sQuote('ouchtree'), ' class.',
- '\nYou can generate a tree of this class by calling ', sQuote('ouchtree()'), '.', sep = ""))
- ## get the vectors we need:
- ancestor <- tree at ancestors # class = "character"
- node <- tree at nodes # class = "character"
- species <- tree at nodelabels # class = "character" -- note that nodelabels is more general than this indicates and the name should be changed throughout at some point
- times <- tree at times # class = "numeric"
- ## ------------------ end ouchtree block -------------------
-
- if(class(regimeShiftNodes) == "list") regimeShiftNodes <- unlist(lapply(regimeShiftNodes, mrcaOUCH, tree = tree))
- regimeShiftNodes <- unique(c(as.character(tree at root), regimeShiftNodes))
- if(identical(regimeTitles, NULL)) regimeTitles <- as.character(regimeShiftNodes)
- names(regimeTitles) = as.character(regimeShiftNodes)
- colorsVector = character(length(node))
- for (i in 1:length(ancestor)) {
- # First three lines fill up the vector for nodes that are hit in order
- if (is.na(ancestor[i])) {
- colorsVector[i] = regimeTitles["1"]
- next }
- if (as.character(ancestor[i]) %in% as.character(regimeShiftNodes)) {
- colorsVector[i] = regimeTitles[as.character(ancestor[i])]
- next }
- if (colorsVector[as.integer(ancestor[i])] != "") {
- colorsVector[i] = colorsVector[as.integer(ancestor[i])]
- next }
- # These lines fill up the vector for nodes run reached before their immediate ancestor
- nodeQ = integer(length(node))
- ii = i
- repeat {
- nodeQ = c(ii, nodeQ)
- ii = as.numeric(ancestor[ii])
- if (as.character(ancestor[ii]) %in% as.character(regimeShiftNodes)) {
- colorsVector[ii] = colorsVector[as.integer(ancestor[ii])]
- break}
- if (colorsVector[as.integer(ancestor[ii])] != "") {
- colorsVector[ii] = colorsVector[as.integer(ancestor[ii])]
- break} }
-
- for(j in nodeQ) {
- colorsVector[j] = colorsVector[as.integer(ancestor[j])] }
-
- } # closes for(i in 1:length(ancestor)) loop
-
- # a hack to fix a problem I don't understand... with the undesired side effect that it colors the stem of some subtrees rather than the crown as originally written
- for(i in 1:length(colorsVector)) if(colorsVector[i] == "") colorsVector[i] <- as.character(i)
-
- # colors terminal branches if any terminal branches are in the regimeShiftNodes
- for(i in regimeShiftNodes) if(i %in% tree at term) colorsVector[as.numeric(i)] <- as.character(i)
- colorsVector <- as.factor(colorsVector)
- names(colorsVector) <- tree at nodes
- return(colorsVector) }
-
-regimeMaker <- function(ouchTrees, regMatrix, nodeMembers) {
-## supplants the old 'allPossibleRegimes'
-## takes a list of ouchtree objects, a regimeMatrix ouput, and a list of nodeMembers (the taxa definining each node of interest)
-## Value:
-## regList = a list of nodes defining the change points for each tree (i.e., a list of lists)
-## nodeMatrix = a matrix of trees (rows) by nodes (columns) indicating whether the node is present in each tree
-
- # set up variables
- numTrees <- length(ouchTrees)
- numNodes <- length(nodeMembers)
- if(numNodes != dim(regMatrix)[2]) stop('Number of nodes (columns) in regMatrix must equal number of items in nodeMembers list')
- nodeMatrix <- matrix(NA, nrow = numTrees, ncol = numNodes, dimnames = list(seq(numTrees), dimnames(regMatrix)[[2]]))
- changeNodes <- list(numTrees)
- regList <- list(numTrees)
- regMatrices <- list(numTrees)
-
- # fill outdata
- for(i in seq(numNodes)) nodeMatrix[, i] <- unlist(lapply(ouchTrees, isMonophyletic, taxa = nodeMembers[[i]]))
- for(i in seq(numTrees)) {
- tree <- ouchTrees[[i]]
- regMatrices[[i]] <- regMatrix * as.numeric(matrix(nodeMatrix[i, ], dim(regMatrix)[1], dim(regMatrix)[2], byrow = TRUE)) # multiplies regMatrix by nodes present
- regMatrices[[i]][1:(dim(regMatrices[[i]])[1] - 1), ][which(apply(regMatrices[[i]][1:(dim(regMatrices[[i]])[1] - 1), ], 1, sum) == 0), ] <- rep(NA, numNodes) # set to NA regimes that have no nodes, except for OU1 model
- regMatrices[[i]][duplicated(apply(regMatrices[[i]], 1, as.decimal)), ] <- rep(NA, numNodes) ## set to NA non-unique regimes
- dimnames(regMatrices[[i]]) <- list(seq(dim(regMatrices[[i]])[1]), dimnames(regMatrices[[i]])[[2]])
- numTreeRegs <- dim(regMatrices[[i]])[1]
- treeRegs <- list(numTreeRegs) # this will be assigned to regList[[i]]
- nodesVector <- unlist(lapply(nodeMembers, mrcaOUCH, tree = ouchTrees[[i]])) # as written, gets the MRCA for even invalid nodes just so indexing stays right
- for(j in seq(numTreeRegs)) {
- if(any(is.na(regMatrices[[i]][j, ]))) treeRegs[[j]] <- NA
- else {
- treeRegs[[j]] <- as.factor(paintBranches(c("1", nodesVector[as.logical(regMatrices[[i]][j, ])]), tree))
- names(treeRegs[[j]]) <- tree at nodes
- }
- }
- regList[[i]] <- treeRegs
- }
- regMatrices$overall <- regMatrix # this is the matrix that includes all regimes without regard to any tree
- outdata <- list(regList = regList, nodeMatrix = nodeMatrix, regMatrix = regMatrices)
- return(outdata)
-}
-
-regimeMatrix <- function(n, maxNodes) {
-## recursive function that returns the same thing as oldRegimeMatrix, but much more efficient, at least for small maxNodes
-## actually, it appears to be more efficient even at n = maxNodes
- if(n == 1) return(matrix(1:0, nrow = 2, ncol = 1))
- outmat <- matrix(NA, nrow = 0, ncol = n)
- for (i in 1:(n-1)) {
- temp <- c(rep(0, (i-1)), 1)
- remainder <- n - i
- if (maxNodes > 1 && remainder > 0) {
- nextMat <- regimeMatrix(remainder, maxNodes - 1)
- temp <- cbind(matrix(temp, dim(nextMat)[1], length(temp), byrow = TRUE), nextMat)
- }
- else temp[(i+1):n] <- rep(0, length((i+1):n))
- outmat <- rbind(outmat, temp)
- }
- outmat <- rbind(outmat, c(rep(0, n-1), 1))
- outmat <- rbind(outmat, rep(0,n))
- dimnames(outmat) = list(seq(dim(outmat)[1]), seq(dim(outmat)[2]))
- return(outmat)
-}
-
-as.decimal <- function(n) {
-# takes a binary vector and makes it a decimal
- digits <- length(n)
- result <- 0
- for(i in digits:1) result <- result + n[i] * 2 ^ (digits - i)
- result
-}
More information about the Mattice-commits
mailing list