Skip to content

[R-sig-dyn-mod] Modify times parameter in lsoda at runtime

3 messages · Faelens Ruben, Thomas Petzoldt

Original
Faelens Ruben
# 
Hello,

I have a dynamic model where I do not know my observation times
beforehand. I only have a single event declared at t=10, and this event
determines the next observation time for the model.

I currently have implemented this as follows:

IV <- c(15, 20, 30)
t <- 0
next_t <- 15
while( ! is.na(next_t) ) {
    out <- lsoda(
      y=IV,
      model=MyModel,
      times=c(t, next_t),
      parms=NA
    )
    IV <- out[2, 2:5]
    t <- next_t
    next_t <- if( IV[1] < 10 ) { t + 1 * IV[1] } else { t + 30 * IV[1] }
}

Unfortunately, this is very slow, because LSODA needs to initialize the
stepsize each time.
I can fire an event using the 'events' functionality of LSODA, but I
cannot modify the times parameter or the events array once the ODE
integration is running. Preferably, I would do something like this:

FireEvent <- function(t, y, parms) {
    events$func <<- FireEvent
    events$time <<- if( IV[1] < 10 ) {
	t + 1 * IV[1]
    } else {
      t + 30 * IV[1]
    }
}
IV <- c(15, 20, 30)
events <- as.data.frame(list(func=FireEvent, time=10))
out <- lsoda(
  y=IV,
  model=MyModel,
  times=c(0, Inf),
  parms=NA,
  events=events
)    

Is it at all possible to modify the 'times' or 'events' parameter for
LSODA while it is running?

/ Ruben Faelens
Thomas Petzoldt
# 
Dear Ruben,

changing 'times' on the fly is not possible for efficiency reasons.

In contrast to this it *is* possible to trigger events or to stop a 
simulation at any point in time by using a root function. Examples 3-5 
of the deSolve ?events help page may help you.

Also, re-thinking a model approach and some restructuring can often 
simplify things so that 'tricks' can be avoided.

Hope it helps

Thomas
Faelens Ruben
# 
Thank you for your help. I have succeeded into creating this "event
simulator" by using the root finding functionality of LSODA. Thank you
for your help!

For future reference, here is the code:

rootfunc <- function(t,y,p) {
	if( length(EventQueue) == 0)
		stop()
	else
		return min(t-EventQueue$t)
}

eventfunc <- function(t,y,p) {
	# Execute event
	[...]
	# Remove event from event queue
	[...]

	return(y)
}

lsoda(
	y=y,
	times=c(0, 1e8),
	func=Model,
	parms=parms,
	events=list(func=eventfunc, root=TRUE),
	rootfunc=rootfunc
)