Different solutions for in-place and not in-place functions with unconstrained algorithm

[devmotion]

Hi!

As already mentioned in #14, the unconstrained algorithm computes different solutions for problems with in-place and not in-place functions. A small example:

julia> using DelayDiffEq, DiffEqBase, OrdinaryDiffEq, DiffEqProblemLibrary
julia> alg = MethodOfSteps(BS3(); constrained=false)
julia> u₀ = 1.0
julia> prob = prob_dde_1delay_scalar_notinplace(u₀)
julia> sol = solve(prob, alg)
julia> prob2 = prob_dde_1delay(u₀)
julia> sol2 = solve(prob2, alg)
julia> sol.t == sol2.t
false

julia> sol.u == sol2[1, :]
false

julia> sol.errors
Dict{Symbol,Float64} with 3 entries:
  :l∞    => 3.33289e-5
  :final => 1.86351e-5
  :l2    => 1.40652e-5

julia> sol2.errors
Dict{Symbol,Float64} with 3 entries:
  :l∞    => 3.31446e-5
  :final => 1.86005e-5
  :l2    => 1.401e-5

I assume the problem originates from the DDE integrator and the ODE integrator.integrator pointing, among others, both to the same arrays u, uprev, and k in the in-place case. So for some assignments such as https://github.com/JuliaDiffEq/DelayDiffEq.jl/blob/master/src/integrator_interface.jl#L96, which is executed in the first iteration, the value is only updated once for not in-place integrators but immediately in every step for in-place integrators. Moreover, even during calculations in https://github.com/JuliaDiffEq/DelayDiffEq.jl/blob/master/src/integrator_interface.jl#L69 any change of u or k is immediately propagated to integrator.integrator, and thus might also change the interpolation which is not desired, I guess. I was able to reduce the difference of the errors in the example above by completely separating integrator.u and integrator.k also for in-place functions, i.e. by creating a copy of both arrays during initialization and always using recursivecopy! instead of assignments such as integrator.integrator.u = integrator.u, but I got still a tiny difference.

However, I'm still confused by the current implementation of perform_step!. For step sizes below the minimal delay the next step can be directly calculated by the correspond method in OrdinaryDiffEq.jl (https://github.com/JuliaDiffEq/DelayDiffEq.jl/blob/master/src/integrator_interface.jl#L114), whereas for larger step sizes one calculates a first approximation of the next step (also via OrdinaryDiffEq.jl, https://github.com/JuliaDiffEq/DelayDiffEq.jl/blob/master/src/integrator_interface.jl#L69), which is then used for interpolation during the next iteration since integrator.f depends on a HistoryFunction that contains the initial history function, the current solution of integrator.integrator, and the current state of integrator.integrator (which is relevant for interpolation). If this is correct, I have some questions/remarks:

• Can the line https://github.com/JuliaDiffEq/DelayDiffEq.jl/blob/master/src/integrator_interface.jl#L41 removed and https://github.com/JuliaDiffEq/DelayDiffEq.jl/blob/master/src/integrator_interface.jl#L45 changed to integrator.integrator.tprev = integrator.t? To enable extrapolation in the first iteration only integrator.integrator is relevant, so this should be sufficient? And we can remove the redundant line https://github.com/JuliaDiffEq/DelayDiffEq.jl/blob/master/src/integrator_interface.jl#L94, can't we?
• Should integrator.uprev and integrator.integrator.uprev not always equal the same value since also in subsequent iterations the value at time t of the integrator, that was reached in the last step, stays the same? So why do we need https://github.com/JuliaDiffEq/DelayDiffEq.jl/blob/master/src/integrator_interface.jl#L96?
• Why do we need to cache integrator.t, integrator.tprev, and integrator.uprev, even though they are not changed by perform_step! in OrdinaryDiffEq.jl? Can't we remove lines https://github.com/JuliaDiffEq/DelayDiffEq.jl/blob/master/src/integrator_interface.jl#L54-L56 and https://github.com/JuliaDiffEq/DelayDiffEq.jl/blob/master/src/integrator_interface.jl#L107-L109?
• Why are the values such as u, k, fsallast, and EEst that are changed by perform_step! in integrator not copied to integrator.integrator at the end of the method?

[ChrisRackauckas]

Can the line https://github.com/JuliaDiffEq/DelayDiffEq.jl/blob/master/src/integrator_interface.jl#L41 removed and https://github.com/JuliaDiffEq/DelayDiffEq.jl/blob/master/src/integrator_interface.jl#L45 changed to integrator.integrator.tprev = integrator.t? To enable extrapolation in the first iteration only integrator.integrator is relevant, so this should be sufficient? And we can remove the redundant line https://github.com/JuliaDiffEq/DelayDiffEq.jl/blob/master/src/integrator_interface.jl#L94, can't we?

I think that shoiuld work? Of course, if the tests fail after doing that then the answer is no. Likely when first implementing it I added t changes along with each u change to be safe.

Should integrator.uprev and integrator.integrator.uprev not always equal the same value since also in subsequent iterations the value at time t of the integrator, that was reached in the last step, stays the same? So why do we need https://github.com/JuliaDiffEq/DelayDiffEq.jl/blob/master/src/integrator_interface.jl#L96?

It needs to update the current interpolation to use the new values of u and cover the new interval after the first iteration in order to correct it. Otherwise the "history values in the current interval" will not take into account the newest estimates of u.

Why do we need to cache integrator.t, integrator.tprev, and integrator.uprev, even though they are not changed by perform_step! in OrdinaryDiffEq.jl? Can't we remove lines https://github.com/JuliaDiffEq/DelayDiffEq.jl/blob/master/src/integrator_interface.jl#L54-L56 and https://github.com/JuliaDiffEq/DelayDiffEq.jl/blob/master/src/integrator_interface.jl#L107-L109?

Looks like it can be removed and was part of a bad attempt at getting this to work.

Why are the values such as u, k, fsallast, and EEst that are changed by perform_step! in integrator not copied to integrator.integrator at the end of the method?

That looks like a mistake and could be the root of the issue.

[ChrisRackauckas]

For ODEs with BS3: does inplace and out of place give the exact same answer on a size-1 array? I believe the answer should be yes but that's the first thing to check.

[devmotion]

It needs to update the current interpolation to use the new values of u and cover the new interval after the first iteration in order to correct it. Otherwise the "history values in the current interval" will not take into account the newest estimates of u

So for the interpolation to work correctly, it is not sufficient to update u (and not uprev) after every iteration? I thought we want to interpolate in the time interval between t and t+dt and only calculate a new estimate of u at time t+dt in every iteration (which of course needs to be taken into account in the next iteration), but the function value at time t is fixed to uprev.

For ODEs with BS3: does inplace and out of place give the exact same answer on a size-1 array? I believe the answer should be yes but that's the first thing to check.

In the example above the functions of problems prob_dde_1delay_scalar_notinplace and prob_dde_1delay consist only of one component, so u is a size-1 array in case of the in-place function.

[ChrisRackauckas]

So for the interpolation to work correctly, it is not sufficient to update u (and not uprev) after every iteration? I thought we want to interpolate in the time interval between t and t+dt and only calculate a new estimate of u at time t+dt in every iteration (which of course needs to be taken into account in the next iteration), but the function value at time t is fixed to uprev.

This probably needs an extra comment. It's more like this. When dt > smallest delay, the history function is implicit and requires values in [t,t+dt]. But the interpolation integrator(t) will do current interpolation and then will use extrapolation for t beyond integrator.t. This is used here for the startup, and thus looks like this:

0. Previous step ends, the current interval for the interpolation is [t-dt_old,t], and we look at to step from t to t+dt.
1. In the first iteration, we solve using extrapolation given by the interpolation on [t-dt_old,t] to give projected "history values" in [t,t+dt]. With these history values, we solve a perform_step to get u.
2. Since this equation was implicit, we didn't necessarily get the correct values for the "projected history". So now we update t and u. By doing so, the "current interval" in the integrator changes to [t,t+dt], i.e it now uses the value of u, and now the adjacent history of [t-dt_old,t] gives exactly the same values but is handled by the interpolation on sol (through the definition of the history function).
3. We iteratively update u (Picard iteration, or soon to be Anderson) until we solve the implicit equation (implicit only because of those history values).

So it's moreso that the domain of interpolation needs to change after the first iteration. That is the part that is tricky and is why all of those iter >1 checks are there. I hope that makes sense.

In the example above the functions of problems prob_dde_1delay_scalar_notinplace and prob_dde_1delay consist only of one component, so u is a size-1 array in case of the in-place function.

I would still check the ODE solver works out here like you think it does. I am 90% certain it does, but I don't think I CI test every solver to make sure it does.

[devmotion]

Thanks for the explanation! It makes sense, that's the way I understood it mathematically, as well. However, I'm not sure I got it how exactly an integrator uses uprev, u, tprev, and t to either extrapolate or intrapolate.

I would still check the ODE solver works out here like you think it does. I am 90% certain it does, but I don't think I CI test every solver to make sure it does.

Sorry, I misread it in the first place, I didn't get that you were talking about ODE solvers. Yes, this was one of the first things I checked and the ODE solver seems to work in the same way for both in-place and not in-place functions:

julia> using DiffEqProblemLibrary, OrdinaryDiffEq
julia> alg = BS3()
julia> prob_notinplace = ODEProblem(DiffEqProblemLibrary.linear, 1.0, (0.0, 5.0))
julia> prob_inplace = ODEProblem(DiffEqProblemLibrary.f_2dlinear, [1.0], (0.0, 5.0))
julia> sol_notinplace = solve(prob_notinplace, alg)
julia> sol_inplace = solve(prob_inplace, alg)
julia> sol_notinplace.t == sol_inplace.t
true

julia> sol_notinplace.u == sol_inplace[1, :]
true

julia> sol_notinplace.errors
Dict{Symbol,Float64} with 3 entries:
  :l∞    => 0.311186
  :final => 0.311186
  :l2    => 0.0920546

julia> sol_inplace.errors
Dict{Symbol,Float64} with 3 entries:
  :l∞    => 0.311186
  :final => 0.311186
  :l2    => 0.0920546

[devmotion]

Coming back to the question and remarks from above:

Can the line https://github.com/JuliaDiffEq/DelayDiffEq.jl/blob/master/src/integrator_interface.jl#L41 removed and https://github.com/JuliaDiffEq/DelayDiffEq.jl/blob/master/src/integrator_interface.jl#L45 changed to integrator.integrator.tprev = integrator.t? To enable extrapolation in the first iteration only integrator.integrator is relevant, so this should be sufficient? And we can remove the redundant line https://github.com/JuliaDiffEq/DelayDiffEq.jl/blob/master/src/integrator_interface.jl#L94, can't we?

I think that shoiuld work? Of course, if the tests fail after doing that then the answer is no. Likely when first implementing it I added t changes along with each u change to be safe.

All tests still pass after applying those changes, and also my example still fails (with the same errors).

Why do we need to cache integrator.t, integrator.tprev, and integrator.uprev, even though they are not changed by perform_step! in OrdinaryDiffEq.jl? Can't we remove lines https://github.com/JuliaDiffEq/DelayDiffEq.jl/blob/master/src/integrator_interface.jl#L54-L56 and https://github.com/JuliaDiffEq/DelayDiffEq.jl/blob/master/src/integrator_interface.jl#L107-L109?

Looks like it can be removed and was part of a bad attempt at getting this to work.

Same here, changed code passes tests, but my example still fails.

Why are the values such as u, k, fsallast, and EEst that are changed by perform_step! in integrator not copied to integrator.integrator at the end of the method?

That looks like a mistake and could be the root of the issue.

Tried to update those values in integrator.integrator as well. Unfortunately this doesn't make any difference. When thinking about it, this result seems reasonable - we usually work with integrator and its values, so values of integrator.integrator doesn't matter, besides in the methods savevalues!, perform_step!, and postamble!, where the necessary fields of integrator.integrator are updated.

So I think the root of the problem actually might be

I assume the problem originates from the DDE integrator and the ODE integrator.integrator pointing, among others, both to the same arrays u, uprev, and k in the in-place case. So for some assignments such as https://github.com/JuliaDiffEq/DelayDiffEq.jl/blob/master/src/integrator_interface.jl#L96, which is executed in the first iteration, the value is only updated once for not in-place integrators but immediately in every step for in-place integrators. Moreover, even during calculations in https://github.com/JuliaDiffEq/DelayDiffEq.jl/blob/master/src/integrator_interface.jl#L69 any change of u or k is immediately propagated to integrator.integrator, and thus might also change the interpolation which is not desired, I guess.

I see (at least) the following problems:

• Setting integrator.integrator.uprev = integrator.u (https://github.com/JuliaDiffEq/DelayDiffEq.jl/blob/master/src/integrator_interface.jl#L96) after the first iteration leads to two fundamentally different results for in-place and not in-place functions: for not in-place functions ìntegrator.integrator.uprev is updated once after the first iteration and fixed to certain number throughout all iterations whereas for in-place functions the elements of the array (!) integrator.integrator.uprev are updated also in every other iteration when integrator.u is updated

• Calculating new steps via perform_step! in OrdinaryDiffEq.jl (https://github.com/JuliaDiffEq/DelayDiffEq.jl/blob/master/src/integrator_interface.jl#L69) also has different implications: for not in-place functions updates to u, k, EEst, and fsallast are executed after all function evaluations using the inter/extrapolation of integrator.integrator with the updated u of the last iteration step (https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/integrators/low_order_rk_integrators.jl#L38-L45) whereas for in-place functions the elements of the array integrator.u (which equals the arrays integrator.integrator.uprev (by https://github.com/JuliaDiffEq/DelayDiffEq.jl/blob/master/src/integrator_interface.jl#L96) and integrator.integrator.u (by definition)) are updated (https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/integrators/low_order_rk_integrators.jl#L98) before the last function evaluation (https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/integrators/low_order_rk_integrators.jl#L100) which hence already uses the updated value of u instead of the one of the iteration before!

[ChrisRackauckas]

Setting integrator.integrator.uprev = integrator.u (https://github.com/JuliaDiffEq/DelayDiffEq.jl/blob/master/src/integrator_interface.jl#L96) after the first iteration leads to two fundamentally different results for in-place and not in-place functions: for not in-place functions ìntegrator.integrator.uprev is updated once after the first iteration and fixed to certain number throughout all iterations whereas for in-place functions the elements of the array (!) integrator.integrator.uprev are updated also in every other iteration when integrator.u is updated

Yes, now that I think about it, that's it. That's actually really bad. A good example which would break this is the SSPRK methods. In those methods, u is actually used as a cache variable. This means that the during the perform_steps call, the history function would be wildly changing and would make absolutely no sense. In other methods this is less of a big deal. For BS3() it is accidentally helpful, which is why it works. In fact, for any truly FSAL method which doesn't use u as an intermediate cache, you have a line

https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/integrators/low_order_rk_integrators.jl#L98

that updates u before the end (the FSAL property), and then a function call:

https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/integrators/low_order_rk_integrators.jl#L100

At that function call, it will use the updated history function because we have set the pointers and not the values the same. That for sure will cause a difference. Here that's slightly beneficial (updated u is exactly the correct u), but in almost all other cases (i.e. not an FSAL RK method, but I only tested with FSAL RK methods... omg me...) it's not beneficial and will cause the algorithm to fail. This should be changed to keep the pointers separate and copy. recursivecopy! is very cheap anyways and I was silly to avoid it.

And now I see you noticed the same thing 👍 . u and k should be fixed. EEst is a number so it's fine. fsalfirst probably doesn't matter.

[devmotion]

I already tried to separate u and k of both integrators and always use recursivecopy!. The first problem I encountered was that for not in-place functions which operate on arrays, both k contain undefined elements after initialization in line

https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/integrators/low_order_rk_integrators.jl#L3,

and hence it's not possible to use recursivecopy! immediately; but one can work around that by testing for undefined entries with isassigned. The second problem originates from the initialization of k for in-place functions in lines

https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/integrators/low_order_rk_integrators.jl#L50-L54,

since both integrators point to the same cache and so even with recursivecopy! k[1] and k[2] always point to the same arrays, which are used for interpolation during the calculation of the next step. So one either needs two use two different cache for both integrators as well - or one could use the cache just for the DDE integrator and do not initialize the ODE integrator in line

https://github.com/JuliaDiffEq/DelayDiffEq.jl/blob/master/src/integrator_interface.jl#L133.

The second option seems easier, but I'm not sure whether a separate cache would be preferable for some other steps of the solver.

[ChrisRackauckas]

No, a separate cache should not be necessary. After the ODE integrator initializes, it just needs to get integrator.k, avoiding the undefined references. Or it might be better to just initialize fsallast with zeros in initialize! which would stop this problem and be safer anyways.

since both integrators point to the same cache and so even with recursivecopy! k[1] and k[2] always point to the same arrays, which are used for interpolation during the calculation of the next step.

A separate cache shouldn't be necessary anywhere. In every case it should act just like the ODEIntegrator, the only difference being that integrator.k and integrator.u need to be separated from integrator.integrator.k and integrator.integrator.u during the perform_step. So I don't think a separate cache would be a good idea. They can both point to the same cache, but with integrator.k being a recursivecopy(integrator.integrator.k) (and integrator.u broken from integrator.integrator.u as well), and then after steps copy the new values over. We would just need to be careful in event handling to make sure to update t, u, k right after any backtracking occurs.

[devmotion]

No, a separate cache should not be necessary. After the ODE integrator initializes, it just needs to get integrator.k, avoiding the undefined references. Or it might be better to just initialize fsallast with zeros in initialize! which would stop this problem and be safer anyways.

I'm not sure whether I understand you correctly, but I guess the problem I'm dealing with can't be fixed by an initialization of fsallast alone: the problem is that for this special case (not in-place function, u as array such as [1.0]) integrator.k is initialized as array of arrays in which every element is undefined (e.g. for BS3 in line https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/integrators/low_order_rk_integrators.jl#L3). Its elements are only assigned to the arrays of integrator.fsalfirst and integrator.fsallast after the first step, in lines https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/integrators/low_order_rk_integrators.jl#L43-L44. So probably a safe fix would be to copy those lines to the initialization as well and to assign integrator.fsallast = zero(integrator.fsalfirst) before. And one has to do it for every other not in-place cache, too. Should I open a PR which fixes that problem?

A separate cache shouldn't be necessary anywhere. In every case it should act just like the ODEIntegrator, the only difference being that integrator.k and integrator.u need to be separated from integrator.integrator.k and integrator.integrator.u during the perform_step. So I don't think a separate cache would be a good idea. They can both point to the same cache, but with integrator.k being a recursivecopy(integrator.integrator.k) (and integrator.u broken from integrator.integrator.u as well), and then after steps copy the new values over. We would just need to be careful in event handling to make sure to update t, u, k right after any backtracking occurs.

I tried to implement it similarly in the beginning, but somehow I still ended up with different time points and a small difference in the error estimates (same example as above):

julia> sol.errors
Dict(:l∞=>3.6147e-5,:final=>1.91515e-5,:l2=>1.49162e-5)

julia> sol2.errors
Dict(:l∞=>3.61461e-5,:final=>1.91514e-5,:l2=>1.49159e-5)

Hence I was worried there might be different variables pointing to the same array that somehow lead to a different result for in-place functions compared to not in-place functions. The only other thing I could think of was that the call to ode_addsteps! during calculation of the interpolation might already use new values of the cache for in-place functions. On the other hand, since there are no steps added for BS3 and ode_addsteps! can't use the function integrator.integrator.f of the ODE integrator (it still contains the dependency on the history function, so it has the wrong number of arguments) this problem can not occur in this case. However, it might be a problem for other algorithms, I don't know.

I tried also an implementation with two separate caches, but (as expected with the reasoning above) the example fails with the same error estimates. So the still existent discrepancy might be due to a completely different problem, but I'm out of ideas at the moment. Nevertheless, I'll have a look at my code again, and I can also upload a simple version which just copies u and k.

[ChrisRackauckas]

The only other thing I could think of was that the call to ode_addsteps! during calculation of the interpolation might already use new values of the cache for in-place functions. On the other hand, since there are no steps added for BS3 and ode_addsteps! can't use the function integrator.integrator.f of the ODE integrator (it still contains the dependency on the history function, so it has the wrong number of arguments) this problem can not occur in this case. However, it might be a problem for other algorithms, I don't know.

Integrators with lazy evaluation of extra steps just won't work right now: #16. So that can be dealt with separately in the near future.

I tried also an implementation with two separate caches, but (as expected with the reasoning above) the example fails with the same error estimates. So the still existent discrepancy might be due to a completely different problem, but I'm out of ideas at the moment. Nevertheless, I'll have a look at my code again, and I can also upload a simple version which just copies u and k.

PR what you have there since it's a good improvement already. As for the last little bit, the best way is to print everything out in the first step and find out what the first term that's calculated differently is.