Event and root-finding support

[noamross]

I have a system where I use events and root-finding. Specifically, this is a 3D set of ODEs (infection stages, age and sex structure).

I use root-finding to stop the simulation when it appears to have reached a stable cycle. It requires use of dede, as it compares current values to lag values and returns zero when the differences are sufficiently small, like so.

stopfun = function(t, state, parms) {
  if(t < 10) {
    return(1)
  } else if(sum((state - lagvalue(t - 1))^2) > 0.001) {
    return(1)
  } else {
    return(0)
  }
}

I use events to update values at discrete time points. In my system, a fraction of individuals are removed from the population once per year. Also, age classes advance discretely (e.g., all age-1 juveniles are removed and added to age-2 at a fixed time point.

In theory these might be able to be incorporated as additional arguments to the obj$run() command.

[richfitz]

Thanks! This is totally on the list, and I hope to start work on it Monday. I think there's a rootfinding example in the deSolve manual I can mock up.

It might be that the stopping function etc need to be specified in compiled code to match the deSolve interface or to reasonably access variables in the model.

With your cycle detection, how do you know how far back to look? Presumably the t-1 bit needs to correspond exactly to the length of the cycle for this to work, or am I being dense?

I'm having some difficulty with lsoda and the other deSolve solvers with large systems of equations with lags (see #3 for scattered notes). So I need to also pull together a better solution for delay equations, which is slowing me down at the moment.

[noamross]

Example 6 in ?deSolve::events shows a stop function written in R used a compiled code ODE. It looks like you can have either.

Yes, in my case, I'm looking for a cycle of length exactly one. This is because I'm forcing this cycle with annual events and time-dependent seasonal variables (1 = 1 year), so I'm not looking for emergent cyclical patterns.

[richfitz]

Ah, I see - make sense.

Looks like R events with compiled code were added in the latest deSolve release (1.13) which I've not actually updated to! Will see if that can be made to work, but the issue with R events is accessing and manipulating the model simply is going to be tricky.

Hopefully back onto the model Thursday or Friday...