Nothing
# Simulates the within-host coalescent model
# @param times times at which N samples are taken(counted forward in time from infection time)
# @param neg is the product of the within-host effective population size and the generation duration in days
# @return array of size(2N)*3 where each row is a node,the first column indicate the date of the node and the last two columns indicate the two children. This array has internal nodes sorted in order of most recent to most ancient node(and remains so during the algorithm). The last node corresponds to infection time and only has one child
withinhost = function(times,neg) {
prob <- 0
ind=order(times,decreasing=T);tim=times[ind]
n <- length(tim)
nodes <- cbind(0,ind[1],0);#Start with one node at time 0 and with the first isolate connected to it
i <- 2
while (i <= n) {#Graft branches one by one
curt <- tim[i]
fi <- which( nodes[ ,1] < curt );fi<-fi[1]
trunc=0
for (j in (fi:nrow(nodes))) {
trunc=trunc+(curt-nodes[j,1]) * (i-j)
curt <- nodes[j,1]
}
r=rexpT(1/neg,trunc)
prob=prob+r$prob
r=r$x
curt <- tim[i]#Current time:start with date of isolate and go back in time until coalescence happens
fi <- which( nodes[ ,1] < curt );fi<-fi[1]
for (j in (fi:nrow(nodes))) {
if (r > (curt-nodes[j,1]) * (i-j)) {
r <- r-(curt-nodes[j,1]) * (i-j)
curt <- nodes[j,1]
} else {
curt <- curt-r/(i-j)#Found the time for grafting
r <- 0
break
}
}
if (r>0) warning('Warning: r>0 should not happen')
#Create new node
a <- nodes[ ,2:3];a[a >= j + n] <- a[a >= j + n] + 1;nodes[ ,2:3] <- a;#Renumbering according to table insertion in next line
nodes <- rbind(nodes[seqML(1,j-1), ],c(curt,ind[i],0),nodes[seqML(j,nrow(nodes)),])
#Now choose on which branch to regraft amongst the branches alive at time curt
no <- j
side <- 2
#prob <- prob + log(1/(nrow(nodes)-j))
w <- 1 + floor(runif(1) * (nrow(nodes)-j))
while (w > 0) {
no <- no + side-1
side <- 3-side
if (nodes[no,side + 1] <= n ||(nodes[no,side + 1] > n && nodes[nodes[no,side + 1]-n,1] > curt)) {
w <- w-1
}
}
nodes[j,3] <- nodes[no,side + 1]
nodes[no,side + 1] <- n + j
i <- i + 1
}
nodes <- rbind(matrix(0, nrow = n, ncol = 3),nodes)
nodes[1:n,1] <- times
return(list(nodes = nodes,prob = prob))
}
# Simulates the linear within-host coalescent model
# @param times times at which N samples are taken(counted forward in time from infection time)
# @param neg is the reciprocal of the linear growth rate
# @return array of size(2N)*3 where each row is a node,the first column indicate the date of the node and the last two columns indicate the two children. This array has internal nodes sorted in order of most recent to most ancient node(and remains so during the algorithm). The last node corresponds to infection time and only has one child
withinhostLinear = function(times,neg) {
rate=1/neg
prob <- 0
ind=order(times,decreasing=T);tim=times[ind]
n <- length(tim)
nodes <- cbind(0,ind[1],0);#Start with one node at time 0 and with the first isolate connected to it
i <- 2
while (i <= n) {#Graft branches one by one
curt <- tim[i]
fi <- which( nodes[ ,1] < curt );fi<-fi[1]
r=-log(runif(1))
curt <- tim[i]#Current time:start with date of isolate and go back in time until coalescence happens
fi <- which( nodes[ ,1] < curt );fi<-fi[1]
for (j in (fi:nrow(nodes))) {
ma=(i-j)*(log(curt)-log(nodes[j,1]))/rate
if (r > ma) {
r <- r-ma
curt <- nodes[j,1]
} else {
curt2 <- curt*exp(-rate*r/(i-j))#Found the time for grafting
r <- r- (i-j)*(log(curt)-log(curt2))/rate
curt <- curt2
break
}
}
if (r>1e-6) warning(sprintf('Warning: r>1e-6 should not happen, r=%f',r))
#Create new node
a <- nodes[ ,2:3];a[a >= j + n] <- a[a >= j + n] + 1;nodes[ ,2:3] <- a;#Renumbering according to table insertion in next line
nodes <- rbind(nodes[seqML(1,j-1), ],c(curt,ind[i],0),nodes[seqML(j,nrow(nodes)),])
#Now choose on which branch to regraft amongst the branches alive at time curt
no <- j
side <- 2
#prob <- prob + log(1/(nrow(nodes)-j))
w <- 1 + floor(runif(1) * (nrow(nodes)-j))
while (w > 0) {
no <- no + side-1
side <- 3-side
if (nodes[no,side + 1] <= n ||(nodes[no,side + 1] > n && nodes[nodes[no,side + 1]-n,1] > curt)) {
w <- w-1
}
}
nodes[j,3] <- nodes[no,side + 1]
nodes[no,side + 1] <- n + j
i <- i + 1
}
nodes <- rbind(matrix(0, nrow = n, ncol = 3),nodes)
nodes[1:n,1] <- times
return(list(nodes = nodes,prob = prob))
}
seqML <- function(from, to, by=1) {if (from > to) integer(0) else seq.int(from, to, by)}
#truncated exponential distribution
rexpT = function(rate,trunc) {
x=-log(1-runif(1)*(1-exp(-rate*trunc)))/rate
prob=log(rate)-rate*x-log(1-exp(-rate*trunc))
return(list(x=x,prob=prob))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.