0

I'm writing a agent-based model for phytoplankton. Right now it is fairly basic and includes functions for growth, uptake of a resource ("S" in the code), division, mortality, and a function to keep track of the total extracellular resource concentration outside of the cells (q represents the amount of resource in an individual, then I sum the q of all the individuals and subtract it from the total S). It has become quite a time hog from all the for loops taking more than 24 hours to run when it should be much faster. Are there ways that I can optimize my code to help it run faster?

uptake <- function(x){
  vmax <- x[1]
  km <- x[2]
  S <- x[3]
  result <- vmax*(S/(km+S))
  return(result)
}

growth <- function(x){
  mumax <- x[1]
  qnaught <- x[2]
  q <- x[3]
  result <- mumax*(1-(qnaught/q))
  return(result)
}

division <- function(x){
  mu <-x[1]
  m <- x[2]
  result <- mu*m
  return(result)
}

extracellular <- function(x){
  S <- x[1]
  flow_rate <- x[2]
  nutrient_inflow <- x[3]
  quota_per_time <- x[4]
  result <- ((flow_rate*nutrient_inflow)-(flow_rate*S)-quota_per_time) 
  return(result)
}


## Initial agents and time step parameters
num_agents <- 200
num_time_steps <- 400

## Set up data frame for agents
agents <- data.frame("ID" = 1:num_agents, "vmax" = 0.0000014, "km" = 17, 
                     "mumax" = 1.1, "qnaught" = 0.00000036, "cell_size" = 0.00075, 
                     "alive" = 1)
## Set up data frame for state variables
state <- data.frame("S" = 1.2)
agents <- cbind(agents, state)
state$flow_rate <- 0.86
state$nutrient_inflow <- 1.4
## Set up other parameter values
## q represents the amount of resource in an individual
q <- rep(0, num_agents)
size <- rep(0, num_agents)
m0 <- 0.00075
output <- data.frame()
agents_time <- data.frame()

## For loop for model
for (i in 1:num_time_steps) {
  ## mortality of agents in each time step
  for (k in agents$ID){
    random <- runif(1, min = 1, max = 100)
    if (random < 10.0){
      agents[k,7] = 0
    }
  }
  agents_alive <- filter(agents, alive == 1)
  
  ## update S from state data frame
  agents_alive$S <- state$S
  q <- q[1:nrow(agents_alive)]
  ## sum uptake per time step to use for extracellular concentration

  quota_per_agent <- apply(agents_alive[, c(2, 3, 8)], MARGIN = 1, FUN = function(x) uptake(x))
  quota_per_time <- sum(quota_per_agent)
  q <- q + quota_per_agent
  tmp <- cbind(agents_alive, q)
  mu <- apply(tmp[, c(4, 5, 9)], MARGIN = 1, FUN = function(x) growth(x))
  tmp <- cbind(tmp, mu)
  size <- size[1:nrow(agents_alive)]
  size <- size + apply(tmp[, c(6, 10)], MARGIN = 1, FUN = function(x) division(x))

  ## if size calculated by division function is below minimum cell size (0.75), 
  ## replace with minimum cell size
  size <- ifelse(size < 0.00075, 0.00075, size)
  index = 1
  ## division by cell size
  for (j in size){
    if (j > (2*tmp[index,6])){
      size[index] = m0
      q[index] = (q[index]/2)
      last <- tail(agents$ID, n = 1)
      new <- c(last + 1, 0.0000014, 17,
               1.1, 0.00000036, m0, 1, state$S)
      agents <- rbind(agents, new)
      agents_alive <- rbind(agents_alive, new)
      q <- append(q, q[index])
      mu <- append(mu, mu[index])
      size <- append(size, size[index])
    }
    index = index + 1
  }

  ## update extracellular nutrient concentration state variable
  state$quota_per_time <- quota_per_time
  new_S <- state$S + apply(state, MARGIN = 1, FUN = function(x) extracellular(x))
  state <- replace(state, 1, new_S)
  if (state$S < 0){
    state$S <- 0
  }
  #dataframe of outputs for each time step
  out_per_time <- data.frame("ID" = agents_alive$ID, 
                             "vmax" = agents_alive$vmax, 
                             "km" = agents_alive$km,
                             "S" = agents_alive$S, 
                             "mumax" = agents_alive$mumax,
                             "qnaught" = agents_alive$qnaught,
                             "time" = rep(i, nrow(agents_alive)), 
                             "q" = q,
                             "mu" = mu,
                             "cell_size" = size,
                             "alive" = agents_alive$alive)
  
  #row-append to output dataframe that stores outputs for all agents and all time steps
  output <- rbind(output, out_per_time)
}
Melissa Mattwig
  • 21
  • 1
  • 2
    I suspect those `append` and `rbind` calls are a bottleneck, since it is often much slower in R to grow objects than it would be to pre-allocate the spaces and then populate it as you go. For more, see Circle 2 of the R Inferno: https://www.burns-stat.com/pages/Tutor/R_inferno.pdf – Jon Spring Jul 26 '22 at 18:38

2 Answers2

1

Using data.tables to implement the model will speed things up massively. A lot of the time here is spent on copying data within the data.frames and using a lot of memory to do so. data.tables use copy-in-place and much, much faster. I've done my best to understand the logic and replicate it... Where I may be off, hopefully an easy adjustment. Further, a number of for loops have been avoided by using data.table "set" operations and using mapply instead of apply functions.

library(data.table)

uptake2 <- function(vmax,km,S){
  vmax*(S/(km+S))
}

growth2 <- function(mumax,qnaught,q){
  mumax*(1-(qnaught/q))
}

division2 <- function(mu,cell_size,min_cell_size){
  max(mu*cell_size,min_cell_size)
}

extracellular2 <- function(S,flow_rate,nutrient_inflow,quota_per_time){
  max(((flow_rate*nutrient_inflow)-(flow_rate*S)-quota_per_time),0) 
}


## Initial agents and time step parameters
num_agents <- 10
num_time_steps <- 400
min_cell_size <- 0.00075

state <- data.table(S = 1.2)

## Set up data frame for agents
agents <- data.table(ID = 1:num_agents, 
                     vmax = 0.0000014, 
                     km = 17, 
                     mumax = 1.1, 
                     qnaught = 0.00000036, 
                     min_cell_size = min_cell_size, 
                     alive = 1,
                     S = state$S)

state[,flow_rate:= 0.86]
state[,nutrient_inflow:= 1.4]
## Set up other parameter values


output <- data.table(ID = numeric(0),
                     vmax = numeric(0),
                     km = numeric(0),
                     S = numeric(0),
                     mumax = numeric(0),
                     qnaught = numeric(0),
                     time = numeric(0),
                     q = numeric(0),
                     mu = numeric(0),
                     cell_size = numeric(0),
                     alive = numeric(0))


#agents_time <- data.frame()

## For loop for model
for (i in 1:num_time_steps) {
  ## mortality of agents in each time step
  for (k in agents$ID){
    random <- runif(1, min = 1, max = 100)
    if (random < 10.0){
      agents[k,alive:=0]
    }
  }
  agents_alive <- agents[alive==1]

  if (nrow(agents_alive)==0) next;
  
  ## sum uptake per time step to use for extracellular concentration
  agents_alive[,q:=mapply(uptake2,
                          vmax=vmax,
                          km=km,
                          S=state$S)]
  
  quota_per_time <- sum(agents_alive$q)

  agents_alive[,mu:=mapply(growth2,
                           mumax=mumax,
                           qnaught=qnaught,
                           q=q)]
  
  agents_alive[,size:=mapply(division2,
                             mu=mu,
                             cell_size=min_cell_size,
                             min_cell_size=min_cell_size)]

  ## division by cell size
  create_new_agents <- nrow(agents_alive[size > 2*min_cell_size])
  if (create_new_agents > 0) {
    new_agents <- data.table(ID = max(agents$ID)+(1:create_new_agents),
                             vmax = 0.0000014, 
                             km = 17, 
                             mumax = 1.1, 
                             qnaught = 0.00000036, 
                             cell_size = min_cell_size, 
                             alive = 1,
                             S = state$S)

    agents_alive[size > 2*cell_size,
                 `:=`(size=min_cell_size,
                      q=q/2)]
    
    agents <- rbindlist(list(agents,
                             new_agents))
    
    agents_alive <- rbindlist(list(agents_live,
                                   new_agents))
    
  }


  
  ## update extracellular nutrient concentration state variable
  state$quota_per_time <- quota_per_time
  
  state[,S:=extracellular2(S=S,
                           flow_rate=flow_rate,
                           nutrient_inflow=nutrient_inflow,
                           quota_per_time=quota_per_time)]
  
  #row-append to output dataframe that stores outputs for all agents and all time steps
  output <- rbindlist(list(output, 
                           agents_alive[,
                                        .(ID,
                                          vmax,
                                          km,
                                          S,
                                          mumax,
                                          qnaught,
                                          time=i,
                                          q,
                                          mu,
                                          cell_size=size,
                                          alive)]
                           ))
  print(i)
}
Soren
  • 1,792
  • 1
  • 13
  • 16
0

You can likely avoid the internal for loop with the uniformly distributed random variate. You could try using (perhaps with some adjustment) something like:

agents$random <- runif(length(agents$ID), min = 1, max = 100)
agents$alive[agents$random < 10] <- 0

head(agents)
  ID    vmax km mumax qnaught cell_size alive   S    random
1  1 1.4e-06 17   1.1 3.6e-07   0.00075     0 1.2  5.313743
2  2 1.4e-06 17   1.1 3.6e-07   0.00075     1 1.2 78.682158
3  3 1.4e-06 17   1.1 3.6e-07   0.00075     1 1.2 49.320952
4  4 1.4e-06 17   1.1 3.6e-07   0.00075     1 1.2 56.237847
5  5 1.4e-06 17   1.1 3.6e-07   0.00075     0 1.2  5.142588
6  6 1.4e-06 17   1.1 3.6e-07   0.00075     1 1.2 99.543501

If you still need the outer for-loop to iterate at each time step, you can essentially erase the "random" column afterwards with:

agents$random <- NULL
Jeffrey Harkness
  • 96
  • 1
  • 5