I am trying to find the numerical solution for an ODE at time points (t_1,...,t_n) with parameters (theta_1,...,theta_n), where n is in the order of tens of thousands. In fact, this can be seen as n "independent" ODEs. Then, I need to repeat this M times, where M in the order of tens of thousands. That is, one ODE for each time point with different parameters, and only the first n steps are parallelizable.
Is there an efficient way of doing this? Should I parallelize it as a function of t's and theta's? I am currently using a for loop, but the equation I am interested in is slightly more complicated than the example below and n is in the order of tens of thousands. Thus, I would like to make an efficient implementation of this solver. I have also tried a naive parallel implementation of foreach, but the resulting code is even slower.
A reproducible example below, for M = 1, using the logistic ODE (I know this one has an analytic solution, just trying to present a reproducible toy example)
library(deSolve)
# Logistic ODE
logisode <- function(t, y, par) {
# state variables
N <- y[1]
# parameters
lambda <- par[1]
kappa <- par[2]
# model equations
dN <- lambda*N*(1 - N/kappa)
# result
return( list(dN) )
}
lambdas = seq(0.1,1,by=0.1)
kappas = seq(0.1,1,by=0.1) + 1
ttimes <- seq(0.1, 25, length.out = 10)
out <- matrix(0, ncol = 2, nrow = length(lambdas))
for(i in 1:length(lambdas)){
params <- c(lambda = lambdas[i], kappa = kappas[i], h0 = 0.1)
init <- c(N = params[3])
times <- c(0,ttimes[i])
out[i,] <- ode(init, times, logisode, params, method = "lsode")[-1,]
}
plot(out, type = "l")
Here an example for Windows with parApply. The difficulty with Windows is, that the complete data have to be passed to each core. That's why packages and "shared" data need to be defined or loaded within runModel. This is easier for Linux, with mcapply, where the cores have access to global variables.
library(parallel)
runModel <- function(params) {
library(deSolve)
logisode <- function(t, y, par) {
N <- y[1]
lambda <- par[1]
kappa <- par[2]
dN <- lambda * N * (1 - N/kappa)
return( list(dN) )
}
init <- c(N = params[3])
times <- c(0, params[4])
ret <- ode(init, times, logisode, params, method = "lsode")[-1,]
unname(ret) # discard names
}
lambdas <- seq(0.1, 1, by = 0.1)
kappas <- seq(0.1, 1, by =0.1) + 1
ttimes <- seq(0.1, 25, length.out = 10)
params <- cbind(lambda = lambdas, kappa = kappas, h0 = 0.1, times=ttimes)
## single core version
#out <- apply(params, MARGIN = 1, FUN = runModel)
## parallel on 8 cores
cl <- makeCluster(getOption("cl.cores", 8))
out <- parApply(cl, params, MARGIN = 1, FUN = runModel)
stopCluster(cl)
out <- t(out) # transpose matrix
plot(out, type = "l")
There are of also other possibilities, e.g. the future package.