dhsl.R

next_poisson <- function(lambdas) {
    lambda <- sum(lambdas)
    t <- rexp(1, rate = lambda)
    i <- sample(length(lambdas), 1, prob = lambdas / lambda)
    return(c(t, i))
}

duplication <- function(s) {
    i <- sample(length(s), 1, prob = s / sum(s))
    s[i] <- s[i] + 1
    return(s)
}

host_switch <- function(s) {
    i <- sample(length(s), 1, prob = s / sum(s))
    js <- 1:length(s)
    js <- js[js != i]
    if (length(js) > 1)
        j <- sample(js, 1)
    else
        j <- js
    s[j] <- s[j] + 1
    return(s)
}

loss <- function(s) {
    i <- sample(length(s), 1, prob = s / sum(s))
    s[i] <- s[i] - 1
    return(s)
}

events = c(duplication, host_switch, loss)

simulate <- function(s, t_f) {
    t <- 0
    repeat {
        np <- next_poisson(rates)
        t <- t + np[1]
        if (t >= t_f || s == rep(0, length(s)))
            break
        s <- events[[np[2]]](s)
    }
    return(s)
}

lambda <- 0
tau <- 1
mu <- 0

rates <- c(lambda, tau, mu)

state <- c(1, 0)

t_f <- 2

N <- 10000

# states <-vector(mode = "list", length = N)
writeLines(paste("state", "one", "two", sep = "\t"))
for (i in 0:N) {
#    states[i] <- simulate(state, t_f)
    s <- simulate(state, t_f)
    writeLines(paste(i, s[1], s[2], sep = "\t"))
}

# for (i in 1:length(state)) {
#     s <- states[1:length(states)][i]
#     print(s)
#     hist(s, 0:max(s))
# }