Iteratively simulate a parasite population

I'm trying to build a simulation in R to model the effects of different parasite egg laying behaviors on how many offspring the parasite has, given a limited amount of resources. In brief, this simulation iterates through successive generations of the parasite and follows the number of parasites in each category of traits (how many eggs they lay per berry). I'm looking to see if any particular population goes to extinction more often than others.

Question: Is it possible to speed up this simulation? I'm interested in running sensitivity analyses with it under varying parameter values, but to get to a reasonable number of generations (n >= 100) and simulations per parameter set, it takes quite a while to run. I'm fairly comfortable with the apply/map family of functions, but I don't have enough experience with R to know if there are places in this script where I could replace any loops or un-vectorized functions with vectorized versions.

One other note: I'm currently only returning the data frame object generation_tracker, but I will eventually be returning the berry_mat and oviposition_location objects as well.

I've included the code for the simulation below. The code does require a few packages to run and the simulation has two helper functions which are provided at the top of those code chunk. Thanks!

library(tidyverse)
library(magrittr)
initBerry <- function(n_berries) {
berry_mat <- matrix(0, nrow = 4, ncol = n_berries)
return(berry_mat)
}

initWasps <- function(n_wasps,
n_eggs,
eggs_per_berry,
pr_double) {

wasps_df <- data.frame(matrix(NA, nrow = n_wasps, ncol = 5))
colnames(wasps_df) <- c("starting_eggs",
"tot_eggs",
"eggs_per_berry",
"pr_double",
"larvae_survive")
wasps_df$starting_eggs <- n_eggs wasps_df$tot_eggs <- n_eggs
wasps_df$eggs_per_berry <- eggs_per_berry wasps_df$pr_double <- pr_double
return(wasps_df)
}
simMultiGenerations <- function(n_berries,
n_wasps,
n_eggs,
pr_double,
double_rate,
n_gens,
n_sims)
{

# Initialize simulation ---------------------------------------------------

cl <- makeCluster(detectCores() - 1, type = "SOCK")
registerDoSNOW(cl)
output <- foreach(
j = 1:n_sims,
.export = c("initBerry",
"initWasps"),
.packages = c("magrittr",
"tidyverse",
) %dopar%
{
# NOTE: need to assign to something
berry_mat <- initBerry(n_berries = n_berries)
wasp_df <- initWasps(
n_wasps = n_wasps,
n_eggs = rtpois(n_wasps, n_eggs,a = 4, b = 20),
eggs_per_berry = sample(1:4, n_wasps, replace = T),
pr_double = pr_double
)

# Track time with t, as t increases, increase pr_double
t <- 0
# get a parameter for logistic equation for probability of laying
a <- (1 - pr_double) / pr_double

# nested list to store which wasps oviposit in which seeds
oviposition_location <- vector("list", n_berries)
oviposition_location %<>% map(function(.x) {
.x <- vector("list", 4)
})
generation_tracker <- data.frame(
generation = 1:n_gens,
eggs_per_berry1 = NA,
eggs_per_berry2 = NA,
eggs_per_berry3 = NA,
eggs_per_berry4 = NA
)
for (g in 1:n_gens) {
while (any(wasp_df$tot_eggs > 0)) { # Repeat loop while any wasps have eggs, stop when no eggs remain for (i in 1:nrow(wasp_df)) { # Oviposition happens discretely, one wasp at a time if (wasp_df$tot_eggs[i] == 0) {
# If the current wasp has no eggs, jump to next wasp
next
} else {
# Pick which berry to oviposit in. According to Etsuro, this is random
berry <-
sample(ncol(berry_mat),
size = 1,
replace = TRUE)
# Pick which seed to oviposit in.
if (wasp_df$tot_eggs[i] >= wasp_df$eggs_per_berry[i]) {
seed <- sample(
nrow(berry_mat),
size = wasp_df$eggs_per_berry[i], replace = FALSE ) } else { seed <- sample(nrow(berry_mat), size = 1, replace = FALSE) } # If the seed is empty, oviposit and subtract an egg for (j in 1:length(seed)) { temp_seed <- seed[j] if (berry_mat[temp_seed, berry] == 0) { berry_mat[temp_seed, berry] <- berry_mat[temp_seed, berry] + 1 wasp_df$tot_eggs[i] <-
wasp_df$tot_eggs[i] - 1 oviposition_location[[berry]][[temp_seed]] <- c(i, oviposition_location[[berry]][[temp_seed]]) } else { # Wasp decides whether to oviposit in already oviposited in egg lay_egg <- rbinom(1, 1, prob = wasp_df$pr_double[i])
# If 0, wasp does not oviposit
# If 1, wasp oviposits
if (lay_egg == 1) {
berry_mat[temp_seed, berry] <- berry_mat[temp_seed, berry] + 1
wasp_df$tot_eggs[i] <- wasp_df$tot_eggs[i] - 1
oviposition_location[[berry]][[temp_seed]] <-
c(i, oviposition_location[[berry]][[temp_seed]])
}
}
}
}
# increment time forward for the purpose of increasing oviposition rate
t <- t + 1

if (t %% nrow(wasp_df) == 0) {
# Every time all the wasps have a chance to oviposit, increase
# their likelihood to oviposit in an already occupied seed
wasp_df$pr_double <- 1 / (1 + a * exp(-double_rate * (t / nrow(wasp_df)))) } } } # End while --------------------------------------------------------------- oviposition_location <- oviposition_location %>% modify_depth(2, function(.x) { if (is.null(.x)) { .x <- 0 } else { # Randomly select which parasitoid develops .x <- .x[sample(length(.x), 1)] } }) id <- unlist(oviposition_location) wasp_larvae <- as.data.frame(table(id)) if (any(wasp_larvae$id == 0)) {
wasp_larvae <- wasp_larvae[-1, ]
}
wasp_larvae$id <- as.integer(as.character(wasp_larvae$id))
wasp_df$larvae_survive[wasp_larvae$id] <- wasp_larvae$Freq wasp_df$larvae_survive <- ifelse(is.na(wasp_df$larvae_survive), 0, wasp_df$larvae_survive)
wasp_df %<>% mutate(prop_surv = larvae_survive / starting_eggs)
wasp_df$eggs_per_berry <- factor(as.character(wasp_df$eggs_per_berry),
levels = c("1", "2", "3", "4"))
wasp_df %<>% complete(eggs_per_berry, fill = list(larvae_survive = 0))
n_offspring <- wasp_df %>%
group_by(eggs_per_berry) %>%
summarise(offspring = sum(larvae_survive))
generation_tracker[g, 2:5] <- t(n_offspring[, 2])
frequency_eggs_p_berry <- c(rep(1, n_offspring$offspring[1]), rep(2, n_offspring$offspring[2]),
rep(3, n_offspring$offspring[3]), rep(4, n_offspring$offspring[4]))
n_eggs <- rtpois(sum(n_offspring$offspring), n_eggs,a = 4, b = 20) wasp_df <- data.frame( n_wasps = 1:sum(n_offspring$offspring),
starting_eggs = n_eggs,
tot_eggs = n_eggs,
eggs_per_berry = sample(frequency_eggs_p_berry,
sum(n_offspring$offspring), replace = FALSE), pr_double = rep(pr_double, sum(n_offspring$offspring)),
larvae_survive = NA
)
oviposition_location <- vector("list", n_berries)
oviposition_location %<>% map(function(.x) {
.x <- vector("list", 4)
})
}

# End generations loop ----------------------------------------------------

generation_tracker

}

# End simulation loop (end foreach) ---------------------------------------

stopCluster(cl)
closeAllConnections()
return(output)
}

simMultiGenerations(n_berries = 100,
n_wasps = 20,
n_eggs = 12,
pr_double = 0.1,
double_rate = 0.1,
n_gens = 100,
n_sims = 10)


• You badly need more functions.

• apply functions will not speed up your code that much compared to for loops, but they will make it much more readable, and then you'll be able to spot inefficiencies, in this case though you will more often need purrr::accumulate, or Reduce with accumulate = TRUE, as a given iteration is using the result of the previous iteration. Chances are that all the objects that you initiate with empty values should probably be outputs from an accumulate call of some sort.

• Put all objects that belong to one iteration together inside a list, so you can do things like your_list <- ovoposite(your_list).

• I don't think generation tracker should be initiated as a big data.frame, it should be a list that you'd get from accumulate, you can use bind_rows on it afterwards to get a data.frame. for (j in 1:length(seed)) {temp_seed <- seed[j] should be for (temp_seed in seed)

• I think some of the sampling and some other operations can be done outside of loops in a vectorized way, but i'm never too sure because I can't keep track of what is modified in a loop or not, functions will help there as well.

• Data wrangling will be faster using data.table, loops will be faster with Rcpp but you really need to structure the code much more before.

• for (j in 1:length(seed)) {temp_seed <- seed[j] should be for (temp_seed in seed)

• in the if next sequence, skip the else and save a pair of brackets.

• All of this (beside the data.table / Rcpp mentions) won't speed up your code that much, but it's crucial to make the logic of your algorithm appear clearly, and have data structure and functions that follow that logic, then we might find some tricks to speedup the code by leveraging some vectorized functions.

Edit : more on accumulate

I can't build an example directly from your code because it's a bit too complex, but hopefully this will help.

This is the kind of code that you produced:

x   <- c(10,20,30)
res <- numeric(3)
a   <- 1
b   <- 5

res[[1]] <- a * b + x[1]
for(i in 2:length(res)){
a = a+1 # using the previous value
b = b^2
res[[i]] <- a * b + x[i] - res[[i-1]]
}
res
# [1]   15   55 1850


Let's simplify this, a and b could be defined as vectors outside of the loop more efficiently:

a <- 1:3
b <- accumulate(1:2, ~.^2,.init=5)
# or better b <- 5^(2^(0:2)), i needed to think deeper about the math but it's faster, and this intermediate step helped do it
library(purrr)
abx <- a * b + x
res <- accumulate(abx[-1], ~ .y - .x, .init = abx[1])
res
# [1]   15   55 1850


This will run much faster too, not because of accumulate which just hides a regular loop, but because it forced me to think more clearly about what are my inputs and outputs and vectorize what could be.

We can make it more readable by giving a name to the function that we use to get res, we'll place this function at the top of the R file with the other and we'll document it, which will be easy because they are few parameters and it does something simple. :

iterate_on_result <- function(prev, abx_i) abx_i - prev


The in our code we'll just call :

res <- accumulate(abx[-1], iterate_on_result, .init = abx[1])


Which makes it clear what is the input and what is modified, and makes it easier to simplify the rest of the code.

• Thanks for taking the time to provide these critiques - I know I've definitely got a lot of room to improve my coding abilities. I'll work on trying to incorporate more functions into this to make it easier to read, but in the mean time I'm having trouble understanding how accumulate might fit into this simulation, since I haven't used that function before. Would you be able to provide an example from something in my script that I could use accumulate for? That might be hard, given how hard to read this is... – Jake Sep 19 '18 at 15:17
• Jake it's not an easy algorithm and good job putting this together, by simplifying your code into little blocks you'll also get code that is easier to refactor, hopefully the info I added will help you understand accumulate, reduce is just the same but ouputs the last result only. apply family functions are to be applied on a vector or list when computations are independent from each other, and you might also want to check out replicate which runs a (typically random) operation several times. for loops are not a bad thing, but nested for loops usually are. – Moody_Mudskipper Sep 19 '18 at 16:43