richfitz
commented
1 year ago unfortunately this is not that surprising - if the model is dominated by a matrix multiplication, then the version that uses a linear algebra library will be much faster.
Supporting this properly has been on the back burner for a long time (#38, #134, #213 - these mostly concern multinomial distributions but the syntactic issue in #134 is shared and is the primary blocker). The actual calling convention is not that bad, though it does mean that models need to have a working copy of gfortran to compile which is quite annoying in practice, particularly for people on macs
CGMossa
Benchmarking an R+deSolve code against the equivalent odin code yielded a surprising result:
I suspect the culprit is the lack of matrix multiplication in
odin(or maybe I don't know how to invoke it). In thedeSolvepart:However, in the
odinpart I do:Here I've omitted the parts I don't think are necessary.
I'll work on a minimal testcase to check if this is indeed the problem.
Under details, I have more complete excerpts of my code:
Details
```r odin::odin({ foi[, ] <- foi_mat[i, j] * I[j] delta_transmission[] <- transmission_rate * S[i] * I[i] + transmission_rate * S[i] * sum(foi[i,]) delta_recovery[] <- recovery_rate * I[i] deriv(S[]) <- -delta_transmission[i] deriv(I[]) <- +delta_transmission[i] - delta_recovery[i] deriv(R[]) <- delta_recovery[i] source_id <- user() target_id <- user() Total[] <- S[i] + I[i] + R[i] output(source_prevalence) <- I[as.integer(source_id)] / Total[as.integer(source_id)] output(target_prevalence) <- I[as.integer(target_id)] / Total[as.integer(target_id)] transmission_rate <- user(0.05) recovery_rate <- user(0.01) S0[] <- user() I0[] <- user() foi_mat[,] <- user() initial(S[]) <- S0[i] initial(I[]) <- I0[i] initial(R[]) <- Total[i] - S0[i] - I0[i] dim(S0) <- user() dim(S) <- N dim(I) <- N dim(R) <- N dim(I0) <- N dim(foi_mat) <- c(N, N) dim(foi) <- c(N, N) dim(Total) <- N dim(delta_transmission) <- N dim(delta_recovery) <- N N <- length(S0) }, verbose = TRUE, validate = TRUE, target = "c", pretty = TRUE, skip_cache = FALSE) -> model_generator model <- model_generator$new(S0 = site_S, I0 = site_I, foi_mat = foi_matrix, source_id = as.integer(source_site_id), target_id = as.integer(target_site_id)) model$set_user(transmission_rate = 0.05, recovery_rate = 0.01) ``` Benchmarking: ```r bench::mark( odin_c = model$run(0:100), deSolve_r = { site_R <- site_S site_R[] <- 0 deSolve::ode( y = c(S = site_S, I = site_I, R = site_R), times = 0:100, func = function(time, state, parms) { with(parms, { S <- state[1:N] I <- state[(N + 1):(2 * N)] R <- state[(2 * N + 1):(3 * N)] Total <- S + I + R between_sites <- transmission_rate * S * (foi_matrix %*% I) source_prevalence <- I[[source_id]] / Total[[source_id]] target_prevalence <- I[[target_id]] / Total[[target_id]] list(c( dS = -transmission_rate * S * I - between_sites, dI = +transmission_rate * S * I + between_sites - recovery_rate * I, dR = recovery_rate * I ), source_prevalence = source_prevalence, target_prevalence = target_prevalence) }) }, parms = list( transmission_rate = 0.05, recovery_rate = 0.01, source_id = as.integer(source_site_id), target_id = as.integer(target_site_id), N = length(site_S) ) ) }, check = FALSE ) %>% print() ```