[R-sig-dyn-mod] DeSolve Lagvalue Error: illegal input, lag too large

[Mike Blazanin]

Hi all, I'm having a repeated issue with the delay differential equation
solver. I've included a reproducible example at the end, but here's some
background:

I'm modeling a system with three populations: susceptible hosts, infected
hosts, and pathogens. You can see the derivs() function in my code, but
their interactions are governed by a few parameters. The delay equation
comes in because infected hosts die and release new pathogens a discrete
amount of time after they were infected. The amount of time is governed by
the parameter tau.

In general, because the system includes no decay of the pathogen
population, all simulations should end with extinction of the susceptible
and infected hosts. However, I'm having issues with some of the simulations
giving an error before reaching that point. The error is always the same:
Error in lagvalue(t - parms["tau"], 1) : illegal input in lagvalue -
lag, __, too large, at time = _______

I've made sure to write the derivs() function so that the lag is only
called after we're at least tau time from the start of the simulation. I've
also checked that the time and tau values passed to lagvalue are allowable
values (e.g. by having it print them out, not included in the code here).
I've also gone to what I think is the original C code
<https://r-forge.r-project.org/scm/viewvc.php/*checkout*/pkg/deSolve/src/lags.c?revision=534&root=desolve>
to check when the error is returned, and it seems to be saying that I'm
trying to ask for a lag value before the start of the simulation. However,
best I can tell I'm not asking for lag values before the simulation, and in
some cases (e.g. example #2 in code below) it's very clear that the lag
values being asked for are quite real.

In some cases, when this error arises I seem to be able to fix it by
shortening the step size in the times sequence. However, in other cases
that doesn't seem to help. Moreover, I really can't figure out any rhyme or
reason to why it happens sometimes and not others, so I can't do anything
to try to prevent it systematically. I'd really like to understand when and
why this issue might arise, so any help would be appreciated!

Thanks!
Mike Blazanin

Code showing errors (note that code to plot successful output stored in
yout is included at very end):

#Minimal reproduction of lag too large error
library(deSolve)
library(reshape2)
library(ggplot2)

derivs <- function(t, y, parms) {
  #The derivs function must return the derivative of all the variables at a
  # given time, in a list

  #Issue warning about too small/negative yvals
  if (any(y < 0)) {
    warning(paste("pop(s)",
                  paste(which(y < 0), collapse = ","),
                  "below 0, treating as 0, returning dY = 0"))
  }

  #Save "true" y values for future reference
  ystart <- y

  #Set negative y values to 0 so they don't affect the dN's
  y[y < 0] <- 0

  #Create output vector
  dY <- c(S = 0, I = 0, P = 0)

  ##Calculate dS

  #Old (exponential dS/dt)
  #dS/dt = rS - aSP
  #dS <- parms["r"] * y["S"] - a * y["S"] * y["P"]

  #New (logistic dS/dt)
  #dS/dt = rS((K-S/K)) - aSP
  dY["S"] <- parms["r"] * y["S"] * ((parms["K"] - y["S"])/parms["K"]) -
    parms["a"] * y["S"] * y["P"]

  ##Calculate dI
  #dI/dt = aSP - aS(t-tau)P(t-tau)
  if (t < parms["tau"]) {
    dY["I"] <- parms["a"] * y["S"] * y["P"]
  } else {
    dY["I"] <- parms["a"] * y["S"]*y["P"] -
      parms["a"] * lagvalue(t - parms["tau"], 1)*lagvalue(t - parms["tau"],
3)
  }

  ##Calculate dP
  #dP/dt = baS(t-tau)P(t-tau) - aSP
  if (t < parms["tau"]) {
    dY["P"] <- -parms["a"] * y["S"] * y["P"]
  } else {
    dY["P"] <- parms["b"] * parms["a"] *
      lagvalue(t-parms["tau"], 1)*lagvalue(t-parms["tau"], 3) -
      parms["a"]*y["S"]*y["P"]
  }

  #Issue warning about too large pop
  if (any(y > 10**100)) {
    warning(paste("pop(s)",
                  paste(which(y > 10**100), collapse = ","),
                  "exceed max limit, 10^100, returning dY = 0"))
    dY[y > 10**100] <- 0
  }

  #Change dY for too small pops
  dY[ystart < 0] <- 0

  #From documentation: The return value of func should be a list, whose
first
  #element is a vector containing the derivatives of y with respect to time
  return(list(dY))
}

yinit <- c(S = 10**6,
           I = 0,
           P = 10**5)

#First example
params <- c(r = .04, a = 10**-12, b = 5, tau = 10, K = 10**7)
times <- seq(0, 400, 0.003125)
yout <- as.data.frame(
  dede(y = yinit, times = times, func = derivs, parms = params))

#Second example
params <- c(r = .04, a = 10**-12, b = 5, tau = 50, K = 10**8)
times <- seq(0, 204800, 25.6)
yout <- as.data.frame(
  dede(y = yinit, times = times, func = derivs, parms = params))

times <- seq(0, 73190, 25.6)
yout <- as.data.frame(
  dede(y = yinit, times = times, func = derivs, parms = params))
tail(yout)
#All pops are still positive at 73141, when it says the lag is too large

#Third example
params <- c(r = .04, a = 10**-8, b = 500, tau = 100, K = 10**8)
times <- seq(0, 400, 2.5*10**-2)
yout <- as.data.frame(
  dede(y = yinit, times = times, func = derivs, parms = params))

times <- seq(0, 99, 2.5*10**-2)
yout <- as.data.frame(
  dede(y = yinit, times = times, func = derivs, parms = params))


#Code for plotting population sizes over time
ymelt <- reshape2::melt(data = as.data.frame(yout),
                        id = c("time"),
                        value.name = "Density",
                        variable.name = "Pop")

ggplot(data = ymelt,
       aes(x = time, y = Density+1, color = Pop)) +
  geom_line(lwd = 1.5, alpha = 1) +
  scale_y_continuous(trans = "log10") +
  scale_x_continuous(breaks = seq(from = 0, to = max(ymelt$time),
                                  by = round(max(ymelt$time)/10))) +
  geom_hline(yintercept = 1, lty = 2) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  NULL

	[[alternative HTML version deleted]]