unknown
unknown

Reputation: 883

speed problems with odesolver in R

I have a differential equation model in R that uses the odesolver from the deSolve package. However, at the moment the model is running very slowly. I think this might be something to do with the function that I feed to odesolver being poorly written, but can't figure out what exactly is slowing it down and how I might speed it up. Does anyone have any ideas?

I've made an example that works in a similar way to mine:

library(data.table)
library(deSolve)

matrix_1 <- matrix(runif(100),10,10)
matrix_1[which(matrix_1 > 0.5)] <- 1
matrix_1[which(matrix_1 < 0.5)] <- 0

matrix_2 <- matrix(runif(100),10,10)
matrix_2[which(matrix_2 > 0.5)] <- 1
matrix_2[which(matrix_2 < 0.5)] <- 0

group_ID <- rep(c(1,2), 5)
N <- runif(10, 0, 100000)

Nchange <- function(t, N, parameters) {
  with(as.list(c(N, parameters)), {
    N_per_1 <- matrix_1 * N_per_connection  
    N_per_1[is.na(N_per_1)] <- 0
    total_N_2 <- as.vector(N_per_1)
    if (nrow(as.matrix(N_per_1)) > 1) {
      total_N_2 <- colSums(N_per_1[drop = FALSE])
    }
    N_per_1_cost <- N_per_1
    for (i in possible_competition) {
      column <- as.vector(N_per_1[, i])
      if (sum(column) > 0) {
        active_groups <- unique(group_ID[column > 0])
        if (length(active_groups) > 1){
          group_ID_dets <- data.table("group_ID" = group_ID, "column"= column, "n_IDS" = 1:length(group_ID))
          group_ID_dets$portions <- ave(group_ID_dets$column, group_ID_dets$group_ID, FUN = function(x) x / sum(x))
          group_ID_dets[is.na(group_ID_dets)] <- 0
          totals <- as.vector(unlist(tapply(group_ID_dets$column, group_ID_dets$group_ID, function(x) sum(x))))
          totals[is.na(totals)] <- 0
          totals <- totals*2 - sum(totals)
          totals[totals < 0] <- 0
          group_ID_totals <- data.table("group_ID" = unique(group_ID), "totals" = as.vector(totals))
          group_ID_dets$totals <- group_ID_totals$totals[match(group_ID_dets$group_ID, group_ID_totals$group_ID)]
          N_per_1[, i] <- group_ID_dets$totals * group_ID_dets$portions
        }
      }
    }

    res_per_1 <- N_per_1 * 0.1

    N_per_2 <- matrix_2 * N_per_connection
    N_per_2[is.na(N_per_2)] <- 0 
    res_per_2 <- N_per_2 * 0.1

    dN <- rowSums(res_per_1)  - rowSums(N_per_1_cost * 0.00003) + rowSums(res_per_2)  - 
      rowSums(N_per_2 * 0.00003) - N*0.03

    list(c(dN))
  })
}  # function describing differential equations
N_per_connection <- N/(rowSums(matrix_1) + rowSums(matrix_2))
possible_competition <- which(colSums(matrix_1 != 0)>1)
times <- seq(0, 100, by = 1) 
out <- ode(y = N, times = times, func = Nchange, parms = NULL) 

Upvotes: 0

Views: 373

Answers (1)

Dan
Dan

Reputation: 12074

A good way to identify the bottle neck is with a profiler and the profvis package provides a good way of drilling down into the results. Wrapping your code in p <- profvis({YourCodeInHere}) and then viewing the results with print(p) gives the following insights:

enter image description here

The lines that are taking the most time are (in descending order of time taken):

group_ID_totals <- data.table("group_ID" = unique(group_ID), "totals" = as.vector(totals))

group_ID_dets$portions <- ave(group_ID_dets$column, group_ID_dets$group_ID, FUN = function(x) x / sum(x))

group_ID_dets <- data.table("group_ID" = group_ID, "column"= column, "n_IDS" = 1:length(group_ID))

totals <- as.vector(unlist(tapply(group_ID_dets$column, group_ID_dets$group_ID, function(x) sum(x))))

group_ID_dets$totals <- group_ID_totals$totals[match(group_ID_dets$group_ID, group_ID_totals$group_ID)]

I'm not familiar with the details of your ODE, but you should focus on optimising these tasks. I think the larger issue is that you're running these commands in a loop. Often, you'll hear that loops are slow in R, but a more nuanced discussion of this issue is found in the answers here. Some tips there might help you restructure your code/loop. Good luck!

Upvotes: 1

Related Questions