[Mattice-commits] r29 - pkg/R

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

Author: andrew_hipp
Date: 2008-11-14 20:50:49 +0100 (Fri, 14 Nov 2008)
New Revision: 29

Modified:
   pkg/R/batchHansen.R
   pkg/R/informationCriterion.R
   pkg/R/summarizingAnalyses.R
   pkg/R/treeTraversal.R
Log:
various revisions focused on summarizing analysis from a batchHansen run

Modified: pkg/R/batchHansen.R
===================================================================
--- pkg/R/batchHansen.R	2008-11-14 08:05:41 UTC (rev 28)
+++ pkg/R/batchHansen.R	2008-11-14 19:50:49 UTC (rev 29)
@@ -72,9 +72,9 @@
   treeCounter = 0
   for (i in 1:length(ouchTrees)) {
     tree <- ouchTrees[[i]]
-    regimesList = regimeVectors(tree, cladeMembersList, maxNodes)
+    rl = regimeVectors(tree, cladeMembersList, maxNodes)
     if(identical(regimeTitles, NULL)) {
-      regimeTitles <- as.character(1:length(regimesList))
+      regimeTitles <- as.character(1:length(rl$regimeList))
       if(brown) regimeTitles <- c(regimeTitles, 'brown')
       }
     
@@ -95,12 +95,13 @@
     ## send it off to batchHansen and just stick the results in hansenBatch... this won't work as the number of regimes gets large, 
     ##   so there should be some option here to just hang onto the coefficients for each run (i.e., hang onto 'coef(hansen(...))' rather than 'hansen(...)')
     ##   there could also be an option to save the entire object as a series of files in addition to hanging onto 
-    hansenBatch[[i]] = batchHansen(tree, dataIn, regimesList, regimeTitles, brown, ...)
+    hansenBatch[[i]] = batchHansen(tree, dataIn, rl$regimeList, regimeTitles, brown, ...)
     message(paste("Tree",i,"of",length(ouchTrees),"complete"))}
     
     ## right now no summary is returned; one is needed, summarizing over trees what is summarized for each tree in batchHansen
-  
-  return(list(hansens = hansenBatch, regimes = regimesList)) }
+  outdata <- list(hansens = hansenBatch, regimeList = rl$regimeList, regimeMatrix = rl$regimeMatrix)
+  class(outdata) <- 'hansenBatch'
+  return(outdata)}
 
 batchHansen <-
 # Runs hansen.fit and brown.fit on a tree over a batch of selective regimes

Modified: pkg/R/informationCriterion.R
===================================================================
--- pkg/R/informationCriterion.R	2008-11-14 08:05:41 UTC (rev 28)
+++ pkg/R/informationCriterion.R	2008-11-14 19:50:49 UTC (rev 29)
@@ -17,3 +17,7 @@
   BICwi = as.vector(lapply(deltaBIC, function(x, allDelta) {exp(-0.5 * x) / sum(exp(-0.5 * allDelta))}, allDelta = deltaBIC), mode = "numeric")
   return(list(names = names, u = u, K = K, AIC = AIC, AICc = AICc, BIC = BIC, AICwi = AICwi, AICcwi = AICcwi, BICwi = BICwi)) }
 
+informationCriterion.hansenBatch <- function(hansenBatch)
+## call informationCriterion for a 'hansen.batch' object
+## recall that there are a bunch of trees in here... summarize over these
+

Modified: pkg/R/summarizingAnalyses.R
===================================================================
--- pkg/R/summarizingAnalyses.R	2008-11-14 08:05:41 UTC (rev 28)
+++ pkg/R/summarizingAnalyses.R	2008-11-14 19:50:49 UTC (rev 29)
@@ -21,3 +21,9 @@
     hansenStatsTemp[counter, "maxAICweight"] = max(i[0:16, "AICweight"])}
   return(hansenStatsTemp)}
 
+
+summary.hansenBatch <- function(hansenBatch){
+## items in output: hansens, regimeList, regimeMatrix
+## the summary will sum weights over all nodes over all trees
+
+  }
\ No newline at end of file

Modified: pkg/R/treeTraversal.R
===================================================================
--- pkg/R/treeTraversal.R	2008-11-14 08:05:41 UTC (rev 28)
+++ pkg/R/treeTraversal.R	2008-11-14 19:50:49 UTC (rev 29)
@@ -168,9 +168,9 @@
       regime[[i]] = sort(n[!is.na(n)]) 
     }
     regime[[numberOfRegimes]] = rep("0", times = as.integer(log2(i)) + 1) 
-    if(nodeMatrix == T) {
-      #n <- ifelse(length(changeNodes) == 1, as.numeric(changeNodes), length(changeNodes))
-      regimesNameMatrix = matrix(
+
+    ## make node matrix
+    regimesNameMatrix = matrix(
         data = NA, nrow = numberOfRegimes, ncol = length(changeNodes), dimnames = list(
           as.character(seq(numberOfRegimes)), as.character(seq(length(changeNodes)))
           )
@@ -181,19 +181,20 @@
           else regimesNameMatrix[i,j] = 1 
         }
       }
-      outdata <- regimesNameMatrix
+      outmatrix <- regimesNameMatrix
       if(!identical(maxNodes, NULL)) {
-        outdata <- outdata[apply(outdata,1,sum) <= maxNodes, ]
-        dimnames(outdata)[[1]] = as.character(seq(dim(outdata)[1]))
+        outmatrix <- outmatrix[apply(outmatrix,1,sum) <= maxNodes, ]
+        dimnames(outmatrix)[[1]] = as.character(seq(dim(outmatrix)[1]))
         }
-    }
-    else {
-      outdata <- regime 
+
+    ## prune list
+    outlist <- regime 
       if(!identical(maxNodes, NULL)) {
-        outdata <- outdata[sapply(outdata, length) <= maxNodes]
-        outdata[[length(outdata) + 1]] <- rep("0", length(changeNodes))
+        outlist <- outlist[sapply(outlist, length) <= maxNodes]
+        outlist[[length(outlist) + 1]] <- rep("0", length(changeNodes))
         }
       }
+  outdata <- list(regimeList = outlist, regimeMatrix = outMatrix)
   return(outdata) }
 
 regimeVectors <-
@@ -221,11 +222,14 @@
   #changeNodesVector = vector("character", length(changeNodesList))
   #for (i in 1:length(changeNodesList)) # Changing cladeMemberList to a 1-d vector
   #  {changeNodesVector[i] = changeNodesList[[i]]}
-  allRegimes = allPossibleRegimes(changeNodesVector, maxNodes)
+  apr = allPossibleRegimes(changeNodesVector, maxNodes)
+  allRegimes <- apr$regimeList
+  regimeMatrix <- apr$regimeMatrix
   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))}
-  return(regimePaintings) }
+  outdata <- list(regimeList = regimePaintings, regimeMatrix = regimeMatrix)
+  return(outdata) }