SavingCallback used together with event handling callbacks

[ivborissov]

Event handling callbacks (DiscreteCallback, ContinuousCallback) allow us to save u before and after the affect! with save_positions=(true,true). It would be a nice option for SavingCallback to save values before and after the affect! of event handling callback when those callbacks are used together in a CallbackSet.

[ChrisRackauckas]

Right now it's before the affect! of a DiscreteCallback if you place it before, and after if you place it after. I am not sure if there's a good way to handle this.

[ivborissov]

Then a simple workaround is to place SavingCallback after the DiscreteCallback and add something like tstops-1e-6/tstops+1e-6 to saveat of the SavingCallback

[ChrisRackauckas]

Yeah, you can add that to the affect! of the first callback.

[ivborissov]

I have the following code to use Saving Callback together with ContinuousCallback. The example is simple but it illustrates the issue I still haven't solved. The output order of timepoints in out.t is t=12 - event time before t=10.0 - saveat time. So the SavingCallback saving is done before saveat intepolation... Is there something wrong with my code?

using DiffEqBase, OrdinaryDiffEq, DiffEqCallbacks

ode_func = function(du, u, p, t)
    # static
    k1, comp1 = p
    # dynamic
    A, B = u ./ [comp1, comp1] # factor
    # rules
    v1 = k1 * A * comp1

    du .= [-v1, v1]
end

function saving_func(u, t, integrator)
    # static
    k1, comp1 = integrator.p
    # dynamic
    A, B = u ./ [comp1, comp1] # factor
    # rules
    v1 = k1 * A * comp1

    # return vector
    [A, B, k1, v1] # outputIds
end

function condition(u,t,integrator) # Event when event_f(u,t) == 0
  t-12.0
end

function affect!(integrator)
  saveAddComponents(integrator)
  #k1, comp1 = integrator.p
  #t = integrator.t
  u0 = integrator.sol.prob.u0
  A__, B__ = integrator.u

  A__, B__ = u0
  
  integrator.u .= A__, B__
  saveAddComponents(integrator)
end


event1 = ContinuousCallback(condition,affect!,save_positions=(false,false))

function saveAddComponents(integrator::OrdinaryDiffEq.ODEIntegrator)
    affect! = integrator.opts.callback.discrete_callbacks[1].affect!

    affect!.saveiter += 1
    copyat_or_push!(affect!.saved_values.t, affect!.saveiter, integrator.t)
    copyat_or_push!(affect!.saved_values.saveval, affect!.saveiter, affect!.save_func(integrator.u, integrator.t, integrator),Val{false})

    return nothing
end

out = SavedValues(Float64, Vector{Float64})

#saveat
saveat = collect(0.0:10.0:100.0)

scb = SavingCallback(saving_func, out, saveat = saveat)

# ode problem
prob = ODEProblem(ode_func, [11.0,0.0], (0.0,100.0), [0.01, 1.1], callback = scb)

# solution is stored in out variable
solve(
    prob,
    Tsit5(),
    save_start=true,
    save_end=false,
    save_everystep=false,
    callback = event1
    )

[ivborissov]

And the output out is

SavedValues{tType=Float64, savevalType=Array{Float64,1}}
t:
[0.0, 12.0, 12.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0]
saveval:
Array{Float64,1}[[10.0, 0.0, 0.01, 0.11], [8.8692, 1.1308, 0.01, 0.0975612], [10.0, 0.0, 0.01, 0.11], [9.04753, 0.952467, 0.01, 0.0995229], [9.23116, 0.768838, 0.01, 0.101543], [8.3527, 1.6473, 0.01, 0.0918797], [7.55784, 2.44216, 0.01, 0.0831362], [6.83861, 3.16139, 0.01, 0.0752247], [6.18783, 3.81217, 0.01, 0.0680662], [5.59898, 4.40102, 0.01, 0.0615887], [5.06616, 4.93384, 0.01, 0.0557278], [4.58406, 5.41594, 0.01, 0.0504246], [4.14783, 5.85217, 0.01, 0.0456261]]

[ChrisRackauckas]

@MSeeker1340

[ChrisRackauckas]

Yeah, this will be a difficult issue to solve. The DiscreteCallbacks are done after the continuous callbacks which would be required for correctness in a lot of cases where they modify the value. However, in your case you'd want the SavingCallback to be ran sometimes before the event handling. I say sometimes because if your dt was large enough, then you want to save t=10 but not t=20 and then pull back to t=12 and apply the effect. That is a really difficult interaction to get correct with separate callbacks. Hmm...

[ivborissov]

Interesting... Rewriting affect! this way seems to have solved the issue:

function affect!(integrator)
  scb.affect!(integrator, true)

  u0 = integrator.sol.prob.u0
  A__, B__ = integrator.u

  A__, B__ = u0
  integrator.u .= A__, B__

  scb.affect!(integrator, true)
end

[ChrisRackauckas]

Yes, you need to do the saveat before application of other parts. That's the difficulty in solving this generally.

[ivborissov]

Back to this issue.
My current approach to handle SavingCallback together with ContinuousCallbacks was to add the following lines to affect!(integrator) function like:

function affect!(integrator)
    # check if there are saveat points in (t, t+dt)  before the event
    save_func! = integrator.opts.callback.discrete_callbacks[1].affect!
    save_func!(integrator)
        
    # force save point before event
    save_func!(integrator, true)
    
    # change integrator with event
    # ....

    # force save point after event
    save_func!(integrator, true)
end

However with some ODE problems I have encountered wrong results applying this approch. Here is MWE with SavingCallback used together with two ContinuousCallbacks in a CallbackSet . The computed time point of the second event occures to be wrong. The correct time point can be obtained by making dtmax less than saveat step. However it is a bad workaround and I still don't understand the source of this error.
https://github.com/ivborissov/JuliaODE/blob/master/SavingCallback%20%2B%20Event%20MWE.ipynb

[ivborissov]

Well, I am still thinking of a good way to handle SavingCallback interaction with other callbacks. The workaround I am using is to trigger saving from the affect! function of event handling callback:

function evt_affect!(integrator)
  save_position(integrator) # save before event
  # perform event
  save_position(integrator) # save after event
end

function save_position(integrator)
  affect! = integrator.opts.callback.discrete_callbacks[1].affect!
  affect!.saveiter += 1
  copyat_or_push!(affect!.saved_values.t, affect!.saveiter, integrator.t)
  copyat_or_push!(affect!.saved_values.saveval, affect!.saveiter, affect!.save_func(integrator.u, integrator.t, integrator),Val{false})
end

There is a number of problems with this approach :

1. It ruins the sequence of saving when (t,t+dt) contains both saveat timepoints and "event" timepoints
2. It doesn't work with Sundials

And here is a MWE from SBML cases (sorry, not really "Minimal"):

using DiffEqBase, OrdinaryDiffEq, DiffEqCallbacks, DataFrames, Sundials

function ode_(du, u, p, t)
    (C,S2,k1,k2) = p
    (S1,S3) = u 
    reaction1 = C * k1 * S1 * S2
    reaction2 = C * k2 * S3
    du .= [
      -reaction1+reaction2,  # dS1/dt
      reaction1-reaction2,  # dS3/dt
    ]
end

event1_condition_(u, t, integrator) = 0.75 - u[1]
function event1_assignment_(integrator)
    save_position(integrator) # save point before event
    integrator.p[2] = 1.0     #event
    save_position(integrator)  # save point after event
end

event2_condition_(u, t, integrator) = u[2] - 1.4
function event2_assignment_(integrator)
    save_position(integrator)      # save point before event
    integrator.u[1] = 1.0      # event
    save_position(integrator)      # save point after event
end

event1 = ContinuousCallback(event1_condition_, event1_assignment_,save_positions=(false,false))
event2 = ContinuousCallback(event2_condition_, event2_assignment_, save_positions=(false,false))

saving_(u, t, integrator) = [u[1], integrator.p[2], u[2]] # S1, S2, S3
saveat = [0.0,0.06,0.12,0.18,0.24,0.78,0.84,0.9,0.96,1.02,1.08,1.14]
scb = SavingCallback(saving_, SavedValues(Float64, Vector{Float64}), saveat=saveat)

cbset = CallbackSet(scb, event1, event2)
prob = ODEProblem(ode_, [1.0, 1.0], (0.0, 3.0), [1.0, 2.0, 0.75, 0.25], callback=cbset)

sol = solve(prob, Tsit5(), save_start = false, save_end = false, save_everystep = false)

The output is:

 Row │ t         S1        S2       S3
     │ Float64   Float64   Float64  Float64
─────┼──────────────────────────────────────
   1 │ 0.0       1.0           2.0  1.0
   2 │ 0.06      0.928803      2.0  1.0712
   3 │ 0.246162  0.75          2.0  1.25
   4 │ 0.246162  0.75          1.0  1.25
   5 │ 0.12      0.875017      1.0  1.12498
   6 │ 0.18      0.854871      1.0  1.14513
   7 │ 0.24      0.834435      1.0  1.16556
   8 │ 0.78      0.646588      1.0  1.35341
   9 │ 1.16245   0.6           1.0  1.4
  10 │ 1.16245   1.0           1.0  1.4
  11 │ 0.84      0.638051      1.0  1.36195
  12 │ 0.9       0.630011      1.0  1.36999
  13 │ 0.96      0.62244       1.0  1.37756
  14 │ 1.02      0.61531       1.0  1.38469
  15 │ 1.08      0.608595      1.0  1.39141
  16 │ 1.14      0.602271      1.0  1.39773
  17 │ 2.14329   0.75          1.0  1.65
  18 │ 2.14329   0.75          1.0  1.65

The sequence of solution timepoints is wrong because integration intervals (t,t+dt) here contain both saveats and event timepoints. So, within such intervals saving allways occures after ContinuousCallback. I tried to modify my solution by checking if (t,t+dt) interval contains saveats before event timepoint. So , I modified the evt_affect! function the following way:

function evt_affect!(integrator)
  save_func! = integrator.opts.callback.discrete_callbacks[1].affect! # saving function
  save_func!(integrator) # save timepoints (if any) before event
  save_position(integrator) # save before event
  # perform event
  save_position(integrator) # save after event
end

This seems to solve my issue. The sequence of saved timepoints is now correct:

 Row │ t         S1        S2       S3
     │ Float64   Float64   Float64  Float64
─────┼──────────────────────────────────────
   1 │ 0.0       1.0           2.0  1.0
   2 │ 0.06      0.928803      2.0  1.0712
   3 │ 0.12      0.864703      2.0  1.1353
   4 │ 0.18      0.806992      2.0  1.19301
   5 │ 0.24      0.755033      2.0  1.24497
   6 │ 0.246162  0.75          2.0  1.25
   7 │ 0.246162  0.75          1.0  1.25
   8 │ 0.78      0.64136       1.0  1.35864
   9 │ 0.84      0.633128      1.0  1.36687
  10 │ 0.9       0.625376      1.0  1.37462
  11 │ 0.96      0.618076      1.0  1.38192
  12 │ 1.02      0.6112        1.0  1.3888
  13 │ 1.08      0.604725      1.0  1.39528
  14 │ 1.12616   0.6           1.0  1.4
  15 │ 1.12616   1.0           1.0  1.4
  16 │ 1.14      0.998319      1.0  1.40168
  17 │ 2.14131   0.75          1.0  1.65
  18 │ 2.14131   0.75          1.0  1.65

But... this doesn't work with Sundials.CVODE_BDF : [CVODES ERROR] CVodeGetDky Illegal value for t.t = 0.308972 is not between tcur - hu = 0.309712 and tcur = 0.316776.
@ChrisRackauckas is this workaround acceptable (in the sense it doesn't ruin anything within DiffEq.jl)? What should be changed to make it work with Sundials