Callback Interface

[ChrisRackauckas]

This is to discuss the callback interface separate from the clutter of #5. In there we started discussing a little about coming up with a common callback interface, but didn't get very far.

[ChrisRackauckas]

For reference, here was @mauro3 's suggestion:

callback_step(t1, t2, u1, u2) = ... # simplest callback, no dense output needed
callback_dense(u_dense::Callable) = ... # callback when there is dense output
callback_free_form(alg::Alg,...) # anything else, dispatch on the Algorithm could be useful

I suggested a DSL for the most advanced features. Let's get a middleground that combines them.

[ChrisRackauckas]

Here's the necessary details on my callbacks. My standard callback looks like:

@inline function ODE_DEFAULT_CALLBACK(alg,f,t,u,k,tprev,uprev,kprev,ts,timeseries,ks,dtprev,dt,
  saveat,cursaveat,saveiter,iter,save_timeseries,timeseries_steps,uEltype,ksEltype,
  dense,kshortsize,issimple_dense,fsal,fsalfirst,cache,calck,T,Ts)
  @ode_savevalues
  reeval_fsal = false
  cursaveat,saveiter,dt,t,T,reeval_fsal
end

essentially, it passes a bunch of stuff to run the saving routine. The saving routine implements the full common interface, so it is quite complex:

@def ode_savevalues begin
  if !isempty(saveat) # Perform saveat
    while cursaveat <= length(saveat) && saveat[cursaveat]<= t
      if saveat[cursaveat]<t # If we already saved at the point, ignore it
        saveiter += 1
        curt = saveat[cursaveat]
        ode_addsteps!(k,tprev,uprev,dt,alg,f)
        Θ = (curt - tprev)/dt
        val = ode_interpolant(Θ,dt,uprev,u,kprev,k,alg)
        copyat_or_push!(ts,saveiter,curt)
        copyat_or_push!(timeseries,saveiter,val)
      end
      cursaveat+=1
    end
  end
  if save_timeseries && iter%timeseries_steps==0
    saveiter += 1
    copyat_or_push!(timeseries,saveiter,u)
    copyat_or_push!(ts,saveiter,t)
    if dense
      copyat_or_push!(ks,saveiter,k)
    end
  end
  if !issimple_dense
    resize!(k,kshortsize)
  end
end

The other thing to note is that the interpolations are decomposed since there's no good way to handle the case where they hold immutable values, other than putting them in a big mutable wrapper type which would slow down the routine pretty considerably. So a lot of those values are passed for a reason.

The saving is done in the callback routine because as noted before, event handling has to hijack the saving routine in order to be both correct and fast. So simply removing it is not an option.

I have tried to think about something like, "if the method signature is callback_step(t1, t2, u1, u2), then save outside of the callback", but then what's the callback useful for? You wouldn't want to (or shouldn't want to) do events like this, even in the form callback_dense(u_dense::Callable) (which because of the decomposed nature of the interpolations, would be tough in the first place both for me, and for ODEInterface/Sundials to conform to) you won't be able to (or shouldn't want to) do events. So what's the short-form useful for? Printing intermediate values? Maybe I'm missing something.

[ChrisRackauckas]

One user brought this up as "logging". So I guess some people are interesting in printing intermediate values. Should the common interface just be enhanced for standard logging procedures?

• verbose: Print the values at each verbose_steps steps.
• verbose_steps: How many steps between each print.

If this was included, would there be a need for anything other than the more specific callbacks (for dealing with event handling)?

[ChrisRackauckas]

@tshort might be interested in helping us figure this out.

Maybe the answer could be that callbacks have to be algorithm specific, but event interfaces and logging doesn't have to. Here's what I've been doing. If you check the progress monitoring

https://juliadiffeq.github.io/DiffEqDocs.jl/latest/basics/common_solver_opts.html#Progress-Monitoring-1

you'll see there's a new option: progress_message. The reason for this originally was because Juno added the ability to have a message with the progressbar (when you scroll over it), so I had it show dt, t, and u by default. Now with Juno what we're doing it making this have a fallback to ProgressMeter.jl when Juno isn't present, and display the progress meter and message in the console. Since this can show current times and stuff, is a simple callback callback_step(t1, t2, u1, u2) then not necessary?

If that's not needed, then the only other common thing which callbacks are used for is event handling. But we can make that simply have the options

• event_f : the event function. Positive function with event at 0.
• apply_event!: what to do when an event occurs
• rootfind_event_loc: whether to rootfind the event time
• interp_points: how many interpolated points to check the event trigger
• dt_safety : factor to knock down dt after an event occurs. Useful for discontinuities

For the functions, we should take something like:

• event_f(t,u,du,event_cache) : Is there another piece of information which people use to place events? Note that call overloading can be used to put user caches in here (like pre-allocations). Event cache would be some mutable type <: EventCache which a user can make. Explained below.

• apply_event!(t,u,du,event_cache,cache=[]) : This is used to mutate u (and du for DAEs?). event_cache is mutated in the event_f. cache is optionally an array with pointers to the caches in the ODE/SDE/DAE solver. This allows one to resize the problem on events (see the cell example from here:

https://juliadiffeq.github.io/DiffEqDocs.jl/latest/features/callback_functions.html#Example-2:-Growing-Cell-Population-1

).

With full-featured monitoring and events as part of the common interface, the vast majority of callback usage is then done in the common interface. However, there will still be the option to pass a callback function callback = ... which can be algorithm specific to do whatever crazy specific thing you want with a given algorithm.

Opinions on this?

[ChrisRackauckas]

I did find a problem which needed more than just the events and progress monitoring. However, this would easily be doable if we standardized just the beginning arguments to the callbacks.

SciML/DifferentialEquations.jl#100

what about something like:

callback_free_form(alg::Alg,t,u,du,args...)

where the extra args can be algorithm-specific?

[ChrisRackauckas]

I have found that being able to control the cache in callbacks is very necessary in many cases like hybrid dynamical systems, so I'd suggest it next:

SciML/DifferentialEquations.jl#117

callback_free_form(alg::Alg,t,u,du,cache,args...)

[mauro3]

Could you use a ParameterizedFunction-field for a cache instead?

[ChrisRackauckas]

No, let me describe a little what the cache is. Take a look at RK4:

https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/integrators/fixed_timestep_integrators.jl#L115

The cache is found here:

https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/integrators/fixed_timestep_integrators.jl#L128

cache = (u,tmp,k₁,k₂,k₃,k₄,kprev,uprev)

It's a collection which holds the pointers to every temporary array in the algorithm. Now, this is very non-standard so it should be justified why it's there. The reason why it's there is because it allows you to "fully mutate". Two examples. The first example is the ODE which changes size:

https://juliadiffeq.github.io/DiffEqDocs.jl/latest/features/callback_functions.html#Example-2:-Growing-Cell-Population-1

Notice that in order to have the right size on everything, you have to go through all of the cache arrays and resize all of them in addition to the problem, which is why it's necessary here. The other example is a hybrid dynamical system where one makes a type with the continuous variables are indexed and the discontinuous variables just as fields for the type. In order to fully update so that way the discontinuous changes are seen in all parts of the method, you have to mutate every cache variable as well.

So it is very powerful, but maybe it should be a cache object to handle even more cases nicely via dispatch. Then resize!(cache) would be easier than writing the loop. More importantly, we can upgrade it to handle resizing Jacobians. This is more difficult: if you resize a problem how do you resize the cache for the Jacobian? It's not the same as an array, and it's a difficult calculation when the Jacobian is actually for an N-dimensional array, so it would be nice if resize!(cache) handled that as well.

[mauro3]

I see. Note that this cache is essentially the state variable in PR49's iteration setup.

[ChrisRackauckas]

Something that came up in this thread was the possibility for multiple events:

SciML/OrdinaryDiffEq.jl#16

The idea is just to use an array of functions for event_f and apply_event!, and apply the events sequentially. The implementations that I have all seem to make this pretty easy, so there's no problem with implementation or efficiency.

[ChrisRackauckas]

Intermediate Summary

Since I just posted this on Discourse and want to see if I can get @tshort 's opinion (I know you need events in DASKR.jl) and @ronisbr 's opinion, here's a summary of where the proposal is at now.

Events split from callbacks. The event framework is:

• event_f : the event function. Positive function with event at 0.
• apply_event!: what to do when an event occurs
• rootfind_event_loc: whether to rootfind the event time
• interp_points: how many interpolated points to check the event trigger
• dt_safety : factor to knock down dt after an event occurs. Useful for discontinuities
• save_positions: tuple of 2 booleans for saving before/after event application

The function signatures being:

event_f(t,u,du,event_cache,f)
apply_event!(t,u,du,event_cache,f,cache=nothing)

Cached is described in a previous comment:

#10 (comment)

event_cache lets you calculate something in event_f and use it in the apply_event! (is this necessary?). Additionally, these can be set as tuples for the case where you have multiple events, and those events will be applied sequentially (it would be a tuple because a vector of functions isn't strictly typed, while a tuple of functions is).

Callbacks would have a generic part and an algorithm-specific part. The interface would be:

callback(alg::Alg,t,u,du,f,cache,args...)

This is simple enough to handle most cases where events won't apply, and if anything really specific is needed, the solvers can supply a bunch of args... (which would be documented in their own documentation) from which users can inject sufficiently arbitrary code to their hearts content. However, for most purposes users should be using events since it's at a higher level and so, once all of the solvers implement events, they can act pretty uniformly across the ecosystem whereas I'm convinced callbacks never can.

Two points of contention still:

1. What if the user wants to mutate something like t, u, or du, and those values are immutable? t is usually immutable so this presents a problem. We can just have it as part of the interface that the solvers have to accept mutation of these values, and so if these are immutables, make them wrap it in a size 0 array? Or we can have a separate setup where t is always returned and, if u is immutable, return u and du. This is more similar to how the ODE solvers work themselves, but the caveat here is that it increases user-burden. I'm not quite sure about the performance implications either: using immutables directly in the solver is much faster than size 0 arrays, but callbacks are called only at the end of each iteration and so it may not matter as much.

2. What about solvers which cannot implement parts of this interface? I know that something like Sundials.jl will have trouble properly implementing cache because you'd actually have to change the internal caches of Sundials to make this work properly (is this even possible?). But I wouldn't want to leave cache out of the interface because it's required for non-standard Julia-defined arrays and changing the problem sizes. These things will probably never be supported by Sundials.jl for this reason, but is it fine to then still have it as part of the common interface? Just make a note that cache is only usable in native Julia solvers, and even then restrict that it's only usable in the native Julia solvers which implement the argument properly? Or does it make sense as a keyword argument? Require that users do kwargs... and parse the kwargs... is not nice though.

3. If the functions are in there, is there a need for the du? The reason why I ask is because most methods will not have evaluated this derivative yet so it would incur an extra function evaluation (on any non-FSAL algorithm), and the user can explicitly calculate it themselves using f. Normally this won't be a problem because it can then be used in the next interval, but that's not the case if there's a discontinuity. If there's a discontinuity, then that first du calculation would have to be thrown away, so it would incur an extra function cost without a gain if du is not used in the event handling. However, this does hurt the FSAL case where this du is already calculated so it would be free.