[Mattice-commits] r42 - pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Nov 19 23:51:40 CET 2008
Author: andrew_hipp
Date: 2008-11-19 23:51:40 +0100 (Wed, 19 Nov 2008)
New Revision: 42
Modified:
pkg/R/regimes.R
Log:
all regime painting functions work properly with a list of S4 ouch trees, and nodes missing from a tree are disregarded for that tree
Modified: pkg/R/regimes.R
===================================================================
--- pkg/R/regimes.R 2008-11-19 22:40:19 UTC (rev 41)
+++ pkg/R/regimes.R 2008-11-19 22:51:40 UTC (rev 42)
@@ -1,23 +1,20 @@
# ---------------------------------------------------
# FUNCTIONS FOR PAINTING REGIMES ON AN S4 OUCH TREE #
# ---------------------------------------------------
-
-# Modified from functions used in Hipp 2007 Evolution paper
+# 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
-# FINISH REGIME MAKER
-
-
regimeVectors <-
# 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: list of vectors that can each be plugged directly into OU analysis as the "regimes" argument
-# 19 nov 08: changing to accept a list of trees
-
+# 19 nov 08: changing to accept a list of trees and trimmed down greatly
function(ouchTrees, cladeMembersList, maxNodes = NULL) {
## ------------------ begin ouchtree block -----------------
## check to see if tree inherits 'ouchtree'
@@ -32,28 +29,15 @@
## ------------------ end ouchtree block -------------------
nnode <- length(cladeMembersList)
- regMatrix <- regimeMatrix(n = nnode)
+ regMatrix <- regimeMatrix(n = nnode, maxNodes = maxNodes)
apr = regimeMaker(ouchTrees, regMatrix, cladeMembersList) ## HOLD IT! NOW REGIME MAKER WORKS ON ALL TREES AT ONCE... RETHINK THIS
- # apr$regList = a list of vectors, each indicating changeNodes
- # apr$nodeMatrix = a matrix of trees (rows) by nodes (columns) indicating whether the node is present in each tree
- nodeMatri <- unlist(changeNodesList)
- #changeNodesVector = vector("character", length(changeNodesList))
- #for (i in 1:length(changeNodesList)) # Changing cladeMemberList to a 1-d vector
- # {changeNodesVector[i] = changeNodesList[[i]]}
- allRegimes <- regimesList
- regimePaintings = vector("list", length(allRegimes))
- for (i in 1:length(allRegimes)) {
- allRegimes[[i]] <- c("1", allRegimes[[i]])
- regimePaintings[[i]] <- as.factor(paintBranches(tree, allRegimes[[i]], as.character(allRegimes[[i]])))
- names(regimePaintings[[i]]) <- tree at nodes
- message(paste('Created regime',i))}
- outdata <- list(regimeList = regimePaintings, regimeMatrix = regMatrix)
+ outdata <- list(regList = apr$regList, nodeMatrix = apr$nodeMatrix)
return(outdata) }
paintBranches <-
# Paints branches with regimes changing at nodes specified
# arguments
-# "node" "ancestor" "times" = the standard tree specification vectors of the OUCH-style tree
+# "tree" = OUCH-style (S4) tree
# "regimeShiftNodes" = a vector of nodes at which selective regimes shift: root must be included, 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"
@@ -108,66 +92,6 @@
for(i in 1:length(colorsVector)) if(colorsVector[i] == "") colorsVector[i] <- as.character(i)
return(colorsVector) }
-allPossibleRegimes <-
-# Generates a list of vectors of all possible 2^n regimes for a given list of n ancestor nodes, assuming that each node
-# can either be a change node or not, and that all combinations are meaningful.
-# NOTE: This routine returns all possible permutations, ignoring the fact that the root must be included as a changeNode
-# for paintBranches to work properly. Exclude the root when calling this routine.
-# Arguments:
-# "changeNodes" = vector of all change nodes in all possible scenarios, or number of regimes assumed if nodeMatrix = T
-# "maxNodes" = single number indicating the maximum number of nodes at which a regime can change
-# "nodeMatrix" = indicates whether to return a binary table for interp or a list of changeNode vectors for analysis
-# Value:
-# regimeList = a list of changeNode vectors (assumes type = "character"), one for each possible scenario
-# regimeMatrix = : a binary table indicating whether a regime node is present or absent based on allPossibleRegimes output;
-# presumes nodes are labeled 1:n
-# 10 nov 08: this function now takes over regimeNodes
-function(changeNodes, maxNodes = NULL, nodeMatrix = F) {
- if(!identical(maxNodes, NULL) && maxNodes > length(changeNodes)) warning(paste(sQuote('maxNodes'), 'cannot be larger than the number of nodes; maxNodes ignored'))
- numberOfRegimes = ifelse(length(changeNodes) == 1, 2, 2^length(changeNodes))
- regime = vector("list", numberOfRegimes)
- for (i in 1:(numberOfRegimes - 1)) {
- remainder = i
- n = NULL
- for (j in as.integer(log2(i)):0) {
- if (2^j > remainder) n[j+1] = NA
- else {
- n[j+1] = changeNodes[j+1]
- remainder = remainder %% 2^j
- }
- }
- regime[[i]] = sort(n[!is.na(n)])
- }
- regime[[numberOfRegimes]] = rep("0", times = as.integer(log2(i)) + 1)
-
- ## make node matrix
- regimesNameMatrix = matrix(
- data = NA, nrow = numberOfRegimes, ncol = length(changeNodes), dimnames = list(
- as.character(seq(numberOfRegimes)), as.character(seq(length(changeNodes)))
- )
- )
- for (i in seq(numberOfRegimes)) {
- for (j in seq(length(changeNodes))) {
- if (is.na(match(changeNodes[j],regime[[i]]))) regimesNameMatrix[i,j] = 0 # changed this so that j indexes changeNodes
- else regimesNameMatrix[i,j] = 1
- }
- }
- outmatrix <- regimesNameMatrix
- if(!identical(maxNodes, NULL)) {
- outmatrix <- outmatrix[apply(outmatrix,1,sum) <= maxNodes, ]
- dimnames(outmatrix)[[1]] = as.character(seq(dim(outmatrix)[1]))
- }
-
- ## prune list
- outlist <- regime
- if(!identical(maxNodes, NULL)) {
- outlist <- outlist[sapply(outlist, length) <= maxNodes]
- outlist[[length(outlist) + 1]] <- rep("0", length(changeNodes))
- }
- # }
- outdata <- list(regimeList = outlist, regimeMatrix = outmatrix)
- return(outdata) }
-
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)
More information about the Mattice-commits
mailing list