[Mattice-commits] r199 - in pkg: R misc
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jun 2 23:14:20 CEST 2009
Author: andrew_hipp
Date: 2009-06-02 23:14:16 +0200 (Tue, 02 Jun 2009)
New Revision: 199
Added:
pkg/misc/ouSim.ouchtree-2jun09WorkingSnap-versionFromTextfileSVN.R
pkg/misc/ouSim.ouchtree-2jun09WorkingSnap.R
Modified:
pkg/R/ouSim.ouchtree.R
Log:
something is amiss with ouSim.ouchtree or the plotting function... works locally, not working in compiled version
Modified: pkg/R/ouSim.ouchtree.R
===================================================================
--- pkg/R/ouSim.ouchtree.R 2009-06-01 07:14:00 UTC (rev 198)
+++ pkg/R/ouSim.ouchtree.R 2009-06-02 21:14:16 UTC (rev 199)
@@ -53,6 +53,7 @@
}
names(branchList) <- tree at nodes # branches are indexed by their end nodes
timesList <- lapply(lapply(branchList, length), seq) # creates a list of sequential numbers corresponding to the rnorm vectors
+
for(i in which(sapply(branchList, length) == 0)) timesList[[i]] <- numeric(0) # added b/c seq(numeric(0)) == 1:0
startTimes <- round(unlist(sapply(tree at nodes, ancestorTime, tree = tree)) * steps) ## ASSUMES ULTRAMETRICITY
for (i in 1:length(timesList)) timesList[[i]] <- timesList[[i]] + startTimes[i]
Added: pkg/misc/ouSim.ouchtree-2jun09WorkingSnap-versionFromTextfileSVN.R
===================================================================
--- pkg/misc/ouSim.ouchtree-2jun09WorkingSnap-versionFromTextfileSVN.R (rev 0)
+++ pkg/misc/ouSim.ouchtree-2jun09WorkingSnap-versionFromTextfileSVN.R 2009-06-02 21:14:16 UTC (rev 199)
@@ -0,0 +1,72 @@
+ouSim.ouchtree <- function(object, rootState = 0, alpha = 0, variance = 1, theta = rootState, steps = 1000, ...) {
+## function to plot a simulated dataset under brownian motion or Ornstein-Uhlenbeck (OU) model
+## Arguments:
+## object is an ouch-style (S4) tree
+## alpha and theta are either single values or vectors of length (length(branchList))
+tree <- object
+message(paste("running sim with root =", rootState, ", alpha =", mean(alpha), ", var =", variance, "theta =", mean(theta)))
+
+ ## embedded function---------------------
+ ## could be released to the wild, but more arguments would need to be passed around
+ preorderOU <- function(branchList, tree, startNode, startState, alpha, theta) {
+ ## Recursive function to generate the data under a Brownian motion or OU model
+ ## modified for ouchtree (s4) Dec 08
+ ## branch times back from each tip are in tree at epochs, indexed by tip number
+ ## nodes back from each node (including tips) are in tree at lineages, indexed by node number
+ ## a branch length is the time of a node - the time of its ancestor
+ startBranch <- startNode ## startNode really means start branch... it's the end node of the branch starting this process
+ workingBranch <- branchList[[startBranch]]
+ daughterBranches <- tree at nodes[which(tree at ancestors == startNode)]
+ if(length(workingBranch) == 0) endState <- startState
+ else {
+ for (brStep in 1:length(workingBranch)) {
+ workingBranch[brStep] <-
+ startState + workingBranch[brStep] + alpha[startBranch] / steps * (theta[startBranch] - startState) # denom was mult'd by steps... should be?
+ startState <- workingBranch[brStep]
+ }
+ branchList[[startBranch]] <- workingBranch
+ endState <- branchList[[startBranch]][length(branchList[[startBranch]])]
+ }
+ if(!identical(as.integer(daughterBranches), integer(0))) {
+ for(i in daughterBranches) branchList <- preorderOU(branchList, tree, i, endState, alpha, theta) }
+ return(branchList)
+ }
+ ## --------------------------------------
+
+ ## 1. initialize
+ if(length(alpha) == 1) alpha <- rep(alpha, tree at nnodes)
+ if(length(theta) == 1) theta <- rep(theta, tree at nnodes)
+ brLengths <- c(0, unlist(lapply(2:tree at nnodes, branchLength, tree = tree))) # assumes first node is root; this should be relaxed
+ names(brLengths) <- tree at nodes # branches are indexed by end node
+ names(alpha) <- tree at nodes
+ names(theta) <- tree at nodes
+
+ ## 2. The following creates a list of random draws from the normal distribution, with standard deviation scaled by total
+ ## tree length and the number of draws for each branch equal to the number of steps in that branch.
+ ## If there is a separate variance for each branch, I assume the variance is expressed in terms of total tree-length,
+ ## so the scaling is the same for all branches
+ sdStep <- sqrt(variance / steps)
+ if(length(variance) == 1) branchList <- lapply(round(brLengths*steps), rnorm, sd = sdStep)
+ else {
+ branchList <- vector("list", length(tree at nnodes))
+ for(i in seq(tree at nnodes)) branchList[[i]] <- rnorm(n = brLengths[i], mean = 0, sd = sdStep[i])
+ }
+ names(branchList) <- tree at nodes # branches are indexed by their end nodes
+ timesList <- lapply(lapply(branchList, length), seq) # creates a list of sequential numbers corresponding to the rnorm vectors
+ for(i in which(sapply(branchList, length) == 0)) timesList[[i]] <- numeric(0) # added b/c seq(numeric(0)) == 1:0
+ startTimes <- round(unlist(sapply(tree at nodes, ancestorTime, tree = tree)) * steps) ## ASSUMES ULTRAMETRICITY
+ for (i in 1:length(timesList)) timesList[[i]] <- timesList[[i]] + startTimes[i]
+
+ ## 3. traverse
+ for(i in which(tree at ancestors == tree at root)) { ## calls preorderOU for each descendent from the root.
+ branchList <- preorderOU(branchList, tree, tree at nodes[i], rootState, alpha, theta)
+ }
+
+ value <- list(branchList = branchList, timesList = timesList, steps = steps, parameters = list(rootState = rootState, alpha = alpha, variance = variance, theta = theta))
+ class(value) <- "ouSim"
+ return(value)
+}
+
+nodeTime <- function(node, tree) {
+ return(tree at times[which(tree at nodes == node)])
+}
Added: pkg/misc/ouSim.ouchtree-2jun09WorkingSnap.R
===================================================================
--- pkg/misc/ouSim.ouchtree-2jun09WorkingSnap.R (rev 0)
+++ pkg/misc/ouSim.ouchtree-2jun09WorkingSnap.R 2009-06-02 21:14:16 UTC (rev 199)
@@ -0,0 +1,70 @@
+function(object, rootState = 0, alpha = 0, variance = 1, theta = rootState, steps = 1000, ...) {
+## function to plot a simulated dataset under brownian motion or Ornstein-Uhlenbeck (OU) model
+## Arguments:
+## object is an ouch-style (S4) tree
+## alpha and theta are either single values or vectors of length (length(branchList))
+tree <- object
+message(paste("running sim with root =", rootState, ", alpha =", mean(alpha), ", var =", variance, "theta =", mean(theta)))
+
+## embedded function---------------------
+## could be released to the wild, but more arguments will need to be passed around
+preorderOU <- function(branchList, tree, startNode, startState, alpha, theta) {
+ ## Recursive function to generate the data under a Brownian motion or OU model
+ ## modified for ouchtree (s4) Dec 08
+ ## branch times back from each tip are in tree at epochs, indexed by tip number
+ ## nodes back from each node (including tips) are in tree at lineages, indexed by node number
+ ## a branch length is the time of a node - the time of its ancestor
+ startBranch <- startNode ## startNode really means start branch... it's the end node of the branch starting this process
+ workingBranch <- branchList[[startBranch]]
+ daughterBranches <- tree at nodes[which(tree at ancestors == startNode)]
+ if(length(workingBranch) == 0) endState <- startState
+ else {
+ for (brStep in 1:length(workingBranch)) {
+ workingBranch[brStep] <-
+ startState + workingBranch[brStep] + alpha[startBranch] / steps * (theta[startBranch] - startState) # denom was mult'd by steps... should be?
+ startState <- workingBranch[brStep]
+ }
+ branchList[[startBranch]] <- workingBranch
+ endState <- branchList[[startBranch]][length(branchList[[startBranch]])]
+ }
+ if(!identical(as.integer(daughterBranches), integer(0))) {
+ for(i in daughterBranches) branchList <- preorderOU(branchList, tree, i, endState, alpha, theta) }
+ return(branchList)
+}
+## --------------------------------------
+
+ #debug(preorderOU)
+
+ ## 1. initialize
+ if(length(alpha) == 1) alpha <- rep(alpha, tree at nnodes)
+ if(length(theta) == 1) theta <- rep(theta, tree at nnodes)
+ brLengths <- c(0, unlist(lapply(2:tree at nnodes, branchLength, tree = tree))) # assumes first node is root; this should be relaxed
+ names(brLengths) <- tree at nodes # branches are indexed by end node
+ names(alpha) <- tree at nodes
+ names(theta) <- tree at nodes
+
+ ## 2. The following creates a list of random draws from the normal distribution, with standard deviation scaled by total
+ ## tree length and the number of draws for each branch equal to the number of steps in that branch.
+ ## If there is a separate variance for each branch, I assume the variance is expressed in terms of total tree-length,
+ ## so the scaling is the same for all branches
+ sdStep <- sqrt(variance / steps)
+ if(length(variance) == 1) branchList <- lapply(round(brLengths*steps), rnorm, sd = sdStep)
+ else {
+ branchList <- vector("list", length(tree at nnodes))
+ for(i in seq(tree at nnodes)) branchList[[i]] <- rnorm(n = brLengths[i], mean = 0, sd = sdStep[i])
+ }
+ names(branchList) <- tree at nodes # branches are indexed by their end nodes
+ timesList <- lapply(lapply(branchList, length), seq) # creates a list of sequential numbers corresponding to the rnorm vectors
+ for(i in which(sapply(branchList, length) == 0)) timesList[[i]] <- numeric(0) # added b/c seq(numeric(0)) == 1:0
+ startTimes <- round(unlist(sapply(tree at nodes, ancestorTime, tree = tree)) * steps) ## ASSUMES ULTRAMETRICITY
+ for (i in 1:length(timesList)) timesList[[i]] <- timesList[[i]] + startTimes[i]
+
+ ## 3. traverse
+ for(i in which(tree at ancestors == tree at root)) { ## calls preorderOU for each descendent from the root.
+ branchList <- preorderOU(branchList, tree, tree at nodes[i], rootState, alpha, theta)
+ }
+
+ value <- list(branchList = branchList, timesList = timesList, steps = steps, parameters = list(rootState = rootState, alpha = alpha, variance = variance, theta = theta))
+ class(value) <- "ouSim"
+ return(value)
+}
More information about the Mattice-commits
mailing list