Fix for SurvIT_ode() in case mcmc_size=1

[nkehrein]

Currently, SurvIT_ode() fails if mcmc_site==1 because the apply-function expects a matrix. Thanks.

[pveber]

cc @virgile-baudrot

[virgile-baudrot]

Indeed it's super fast using this approach (10x faster on the small example).

I did a quick check of speed and results.

Speed

Original version:

library(deSolve)
model_SD <- function(t, State, parms, input)  {
  with(as.list(c(parms, State)), {
    TK <- State[1:mcmc_size]
    TD <- State[(mcmc_size+1):(2*mcmc_size)]
    dTK <- kd * (input(t) - TK)     # internal concentration
    dTD <- kk * pmax(TK - z, 0) + hb # risk function
    res <- c(dTK, dTD)
    list(res, conc = input(t))
  })
}

kd = 0.5
hb = 0.2
z = 16
kk = 1.9
mcmc_size = length(kd)
time = c(0,1,2,3)
Cw = c(50,21,3,2)
signal = data.frame(times=time, import=Cw)
sigimp <- stats::approxfun(
  signal$times, signal$import, method = "linear", rule = 2
)
parms = c(mcmc_size,kd,hb,z,kk)
## model
desolve_tktd_solve <- function (i){
  xstart = c(rep(c(D=0),mcmc_size),rep(c(H=0),mcmc_size))
  parms = c(mcmc_size,kd,hb,z,kk)
  deSolve::ode(
    y = xstart,
    times = time,
    func = model_SD,
    parms = parms,
    input = sigimp
  )
}
system.time({ lapply(1:100,desolve_tktd_solve) })

I obtained 0.54 second compared to 0.05 with yours 👍

I isolated your version like this:

library(deSolve)
# with RTools installed:
# system("R CMD SHLIB gutsRedSD.c")
dyn.load("gutsRedSD.dll")

kd = 0.5
hb = 0.2
z = 16
kk = 1.9
mcmc_size = length(kd)

time = c(0,1,2,3)
Cw = c(50,21,3,2)
signal = data.frame(times=time, import=Cw)
times <- signal$times
xstart = c(rep(c(D=0),mcmc_size), rep(c(H=0),mcmc_size))
# ordering of parameters required by compiled function
parms = c(mcmc_size,kd,hb,z,kk)
interpolate_method = "linear" # or "constant"
# solve model
desolve_tktd_solve_C <- function (i){
  xstart = c(rep(c(D=0),mcmc_size),rep(c(H=0),mcmc_size))
  parms = c(mcmc_size,kd,hb,z,kk)
  deSolve::ode(y=xstart,
            times=times,
            parms=parms,
            method="lsoda",
            dllname="gutsRedSD", # specify gutsRedSD
            initfunc="gutsredsd_init",
            func="gutsredsd_func",
            initforc="gutsredsd_forc",
            forcings=signal,
            fcontrol=list(method=interpolate_method,rule=2,ties="ordered"),
            nout=1
  )
}
system.time({ lapply(1:100,desolve_tktd_solve_C) })

Results

However from my script I obtained 2 different results... I've to explore it a little deeper. I probably made something wrong. I keep track here, I'll check tomorrow.

out <- ode(
  y = xstart,
  times = time,
  func = model_SD,
  parms = parms,
  input = sigimp)

out_C <- deSolve::ode(y=xstart,
            times=times,
            parms=parms,
            method="lsoda",
            dllname="gutsRedSD", # specify gutsRedSD
            initfunc="gutsredsd_init",
            func="gutsredsd_func",
            initforc="gutsredsd_forc",
            forcings=signal,
            fcontrol=list(method=interpolate_method,rule=2,ties="ordered"),
            nout=1
)

R> out_C
  time        D         H  3
1    0  0.00000  0.000000 50
2    1 19.67346  1.013234 21
3    2 20.19541  8.729816  3
4    3 13.42958 11.055586  3

R> out
  time         D   H conc
1    0  0.000000 0.0   50
2    1 13.494688 0.2   21
3    2 12.612718 0.4    3
4    3  8.617356 0.6    2