Why does the ode function in deSolve always have an event occur at time = 0?

I am writing a more complex bit of code that wasn't behaving as I expected and after some tinkering have worked out what I think is going wrong: the 'events' in the ode function of deSolve are occuring once more than I would have expected, specifically there is an event at time=0 even when there should not be. I have created a minimal example here:

library(deSolve)

q <- 0

# Define the ODE system
ode_fun <- function(t, y, parms) {
  dy_dt <- -y # y' = -y
  return(list(dy_dt))
}

# Define the event function
event_fun <- function(t, y, parms) {
  print(t)
  q <<- q +1
  y <- y + 1
}

# Set up initial conditions and time points
y0 <- 0.5 # initial value of y
times <- seq(0, 5, by = 0.1) # time points to solve ODEs

# Set up the parameters for the ODE system
parms <- NULL

# Solve the ODE system with the event function
ode_sol <- ode(y = y0, times = times, func = ode_fun, parms = parms,
               events = list(func = event_fun, times = c(2)))

# Plot the solution
plot(ode_sol, xlab = "time", ylab = "y")
print(q)

In the example, the event should only happen once, so q should be 1 at the end, but it is actually 2! Why does this happen? Is there a way to fix it?

Solution

This is nothing to be concerned about. Let's investigate what happens. First I let the event function throw an error to be able to see the call stack:

event_fun <- function(t, y, parms) {
  print(t)
  stop("test0")
  q <<- q +1
  y + 1
}

ode_sol <- ode(y = y0, times = times, func = ode_fun, parms = parms,
               events = list(func = event_fun, time = c(2)))
traceback()
#7: stop("test0") at #3
#6: events$func(time, state, parms, ...)
#5: Func(times[1], y)
#4: eval(Func(times[1], y), rho)
#3: checkEventFunc(Eventfunc, times, y, rho)
#2: lsoda(y, times, func, parms, ...)
#1: ode(y = y0, times = times, func = ode_fun, parms = parms, events = list(func = #event_fun, 
#       time = c(2)))

We see that the function was called within a call to checkEventFunc. Let's take a look at that function:

deSolve:::checkEventFunc
#function (Func, times, y, rho) 
#{
#    tmp <- eval(Func(times[1], y), rho)
#    if (length(tmp) != length(y)) 
#        stop(paste("The number of values returned by events$func() (", 
#            length(tmp), ") must equal the length of the initial conditions vector #(", 
#            length(y), ")", sep = ""))
#    if (!is.vector(tmp)) 
#        stop("The event function 'events$func' must return a vector\n")
#}
#<bytecode: 0x00000173fc02a080>
#<environment: namespace:deSolve>

It is obvious (also from its name) that the only purpose of this function is to check that the event function was defined properly.

Now let's check the other call:

event_fun <- function(t, y, parms) {
  print(t)
  if (t > 0) stop("test1")
  q <<- q +1
  y + 1
}

ode_sol <- ode(y = y0, times = times, func = ode_fun, parms = parms,
               events = list(func = event_fun, time = c(2)))
traceback()
#5: stop("test1") at #3
#4: events$func(time, state, parms, ...)
#3: (function (time, state) 
#   {
#       attr(state, "names") <- Ynames
#       events$func(time, state, parms, ...)
#   })(2, 0.0676677861933153)
#2: lsoda(y, times, func, parms, ...)
#1: ode(y = y0, times = times, func = ode_fun, parms = parms, events = list(func = event_fun, 
#       time = c(2)))

We can see that this is a call to the event function that is actually used by the solver.

You can also compare the solutions with and without event:

event_fun <- function(t, y, parms) {
  y + 1
}
ode_sol_with <- ode(y = y0, times = times, func = ode_fun, parms = parms,
               events = list(func = event_fun, time = c(2)))
ode_sol_without <- ode(y = y0, times = times, func = ode_fun, parms = parms)

plot(ode_sol_with, xlab = "time", ylab = "y", ylim = c(0, 1))
par(new = TRUE)
plot(ode_sol_without, xlab = "time", ylab = "y", ylim = c(0, 1), col = "red", lty = 2)

Clearly, both solutions are identical before the event time.