DiffEqCallbacks problem with DDE

zamk · February 9, 2021, 12:08pm

I have created a SavingCallback which prints out and saves the du, here is my code:

prob = DDEProblem(model,u0,h,(0.0,120.0),p)
alg = MethodOfSteps(Tsit5())
saved_values = SavedValues(Float64,Array{Float64})
saving = SavingCallback((u,t,integrator)->begin
                            du = get_du(integrator)
                            println(du)
                            du
                         end,saved_values,saveat=0.0:1.0:120.0)
positive_domain = PositiveDomain(zeros(6))
sol = solve(prob,alg,reltol=1e-5,abstol=1e-5,verbose=false,
            callback=CallbackSet(positive_domain,saving))
saved_values.saveval

However, I get only copies of the same array in saved_values.saveval:

121-element Array{Array{Float64,N} where N,1}:
 [-0.10130122735023514, -0.12747378214637814, -0.9356362119311361, -0.14477642620179343, -3.214542262837749e-9, 0.0]
 [-0.10130122735023514, -0.12747378214637814, -0.9356362119311361, -0.14477642620179343, -3.214542262837749e-9, 0.0]
 [-0.10130122735023514, -0.12747378214637814, -0.9356362119311361, -0.14477642620179343, -3.214542262837749e-9, 0.0]
 [-0.10130122735023514, -0.12747378214637814, -0.9356362119311361, -0.14477642620179343, -3.214542262837749e-9, 0.0]
 ⋮

whereas the values that are printed out are:

[-9.758880078641988, 9.73926587110036, 0.0, 0.0010190401331306772, 1.731440765537098e-11, 0.0]
[-9.2710620972309, 8.437221518401948, 0.0, 0.02908298475881324, 7.167640978098638e-10, 0.0]
[-9.001511975595097, 7.750094269936375, 0.0, 0.03564680683011713, 1.0608697222391056e-9, 0.0]
[-8.68452543732975, 6.971780891579974, 0.0, 0.03913534215472937, 1.4286156914998178e-9, 0.0]
[-8.313055858601432, 6.100969249011655, 0.0, 0.039895854945804096, 1.8100876888516305e-9, 0.0]
[-7.890979385401242, 5.166178043953941, 0.0, 0.038391043227802596, 2.180181419839846e-9, 0.0]
[-7.419848644971299, 4.192318009030844, 0.0, 0.03509693089910976, 2.5152511392120506e-9, 0.0]

Many thanks in advance for any input!

ChrisRackauckas · February 9, 2021, 12:21pm

can you open an issue in the DelayDiffEq repo? Seems like a bug but I don’t know what it would be off the top of my head.

zamk · February 9, 2021, 12:34pm

I opened the issue, thanks @ChrisRackauckas!

david-pl · February 9, 2021, 1:32pm

@zamk You’re not copying the du, so it is still in use by the solver. So effectively you end up with the derivative of the final step stored in your saved_values. Try doing

saving = SavingCallback((u,t,integrator)->begin
                            du = get_du(integrator)
                            copy(du)
                         end,saved_values,saveat=0.0:1.0:120.0)

zamk · February 9, 2021, 2:14pm

Thanks a lot @david-pl! Since the saving is supposed to be done during solving, I find it strange that it only saves the last step. I was trying to find a workaround by writing du into IOBuffer and collecting the array from there instead of printing it out.

Is there a particular reason why it works like this, is it specific to delay differential equations? Because the following example (for ODE) seems to be working:

saved_values = SavedValues(Float64, Float64)
cb = SavingCallback((u,t,integrator)->get_du(integrator), saved_values,saveat=0.0:0.1:1.0)
f(u,p,t) = 1.01*u
u0 = 1/2
tspan = (0.0,1.0)
prob = ODEProblem(f,u0,tspan)
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8,callback=cb)

david-pl · February 9, 2021, 2:37pm

You’re welcome!

Since the saving is supposed to be done during solving, I find it strange that it only saves the last step.

To be clear: saving is done at every step, however the object that you save is then modified by the solver. To avoid allocations, get_du does not create a new copy, but simply points directly to the du in use by the solver. This is then overwritten at each step so finally you end up with only the derivative of the final step.
Essentially what happens is this:

du = [1.0]
saved = Vector{Float64}[]
for i=1:10
    du[1] = i
    push!(saved, du)
end

i.e. at each iteration du is pushed into the array saved, but since du is modified until the iteration is complete, you end up with a list containing the final du 10 times.

Is there a particular reason why it works like this, is it specific to delay differential equations?

I don’t think so. I think the difference in your examples isn’t ODE vs DDE but rather in-place vs out-of-place. Your OP is in-place, whereas the example you showed in your last post is out-of-place, which does create a copy for du (I think). If you make it in-place, you should see the same behavior.

zamk · February 9, 2021, 2:56pm

Thank you again for the detailed explanation! However, It seems like it saves the last step only if I change the code as follows:

prob = DDEProblem(model,u0,h,(0.0,120.0),p)
alg = MethodOfSteps(Tsit5())
saved_values = SavedValues(Float64,Array{Float64})
saving = SavingCallback((u,t,integrator)->get_du(integrator),
                        saved_values,saveat=0.0:1.0:120.0)
positive_domain = PositiveDomain(zeros(6),save=false)
sol = solve(prob,alg,reltol=1e-5,abstol=1e-5,verbose=false,
            callback=CallbackSet(positive_domain,saving));

I still see this output:

121-element Array{Array{Float64,N} where N,1}:
 [-0.10130122735023514, -0.12747378214637814, -0.9356362119311361, -0.14477642620179343, -3.214542262837749e-9, 0.0]
 [-0.10130122735023514, -0.12747378214637814, -0.9356362119311361, -0.14477642620179343, -3.214542262837749e-9, 0.0]
 [-0.10130122735023514, -0.12747378214637814, -0.9356362119311361, -0.14477642620179343, -3.214542262837749e-9, 0.0]
 [-0.10130122735023514, -0.12747378214637814, -0.9356362119311361, -0.14477642620179343, -3.214542262837749e-9, 0.0]
 ⋮

In addition, I have another problem if I use the copy(du), it seems like solution has the copy of u for each step:

julia> diff(Array(sol),dims=2)
6×80 Array{Float64,2}:
 -0.294771     0.0  -2.48102      0.0  -4.65944      0.0  -5.67162     0.0  -7.07948      0.0  -8.32535     …  -3.15653     0.0  -2.4754      0.0  -1.85137     0.0  -1.39768     0.0  -1.04499     0.0  -0.207552    0.0
  0.294474     0.0   2.45502      0.0   4.49682      0.0   5.2724      0.0   6.26922      0.0   6.92605        -3.50686     0.0  -2.8579      0.0  -2.20407     0.0  -1.70467     0.0  -1.29922     0.0  -0.260717    0.0
  0.0          0.0   0.0          0.0   0.0          0.0   0.0         0.0   0.0          0.0   0.0            -5.41925     0.0  -7.03821     0.0  -7.76922     0.0  -8.08079     0.0  -7.97808     0.0  -1.86149     0.0
  1.54385e-5   0.0   0.00128393   0.0   0.00717181   0.0   0.015062    0.0   0.0253684    0.0   0.0354892      -0.476882    0.0  -0.852066    0.0  -1.08137     0.0  -1.20812     0.0  -1.23096     0.0  -0.288312    0.0
  2.61539e-13  0.0   2.28651e-11  0.0   1.41968e-10  0.0   3.4494e-10  0.0   6.91318e-10  0.0   1.17613e-9      5.13234e-8  0.0   3.39101e-8  0.0   1.51379e-8  0.0  -1.84307e-9  0.0  -1.68774e-8  0.0  -6.00476e-9  0.0
  0.0          0.0   0.0          0.0   0.0          0.0   0.0         0.0   0.0          0.0   0.0         …   0.0         0.0   0.0         0.0   0.0         0.0   0.0         0.0   0.0         0.0   0.0         0.0

david-pl · February 9, 2021, 3:04pm

Yes, that is expected because you’re not copying.

I’m not sure what you mean here. Why look at sol? That will always store u. The du is saved to saved_values.saveval.

zamk · February 9, 2021, 3:12pm

I would like to have both u and du in this case I need to save both of the values by using the SavingCallback I suppose. Still is it not a problem that duplicates of u are saved in sol? Does not it mean that I replaced the du value and interfered with the integration?

david-pl · February 9, 2021, 3:46pm

Either way should work, I think. Either return both a copy of u and du in the callback or use the callback for du only and get u from the sol (in that case you might want to set the same saveat in solve as you’re using in the callback).

As long as you don’t actually modify anything in the callback you shouldn’t be able to interfere with the integration. Returning a copy should be fine. Are you saying your integration changes depending on whether you use copy or not? That really shouldn’t happen.

zamk · February 16, 2021, 10:50am

Here’s how I solved the problem:

saving = SavingCallback((u,t,integrator)->integrator(t,Val{1}),
                        saved_values,saveat=0.0:1.0:120.0)

Instead of using get_du I used integrator(t,Val{1}). I do not have the problems that I mentioned above when I save the du like this.