Strange behavior of save_at=false

[erlebach]

I am getting strange behavior with the different versions of save_at and save_start and save_end. Here is a MWE:

using DifferentialEquations

function exponen!(du, u, param, t)
	du[1] = -u[1]
end

tspan = (0., 3.)
prob = ODEProblem(exponen!, [2.], tspan)
solve(prob, Tsit5(), saveat=[0.2, 1.3], save_start=false, save_end=false)

I expect the solution at the two times 0.2 and 1.3 to be printed. This is what happens. Next, I explicitly add the start and end times to saveat, and specify that save_start and save_end are false. I expect the solution at four times values to be printed. However, only three are printed.

sol = solve(prob, Tsit5(), saveat=[0., 0.2, 1.3, 3.], save_start=false, save_end=false)
println("sol.t: ", sol.t)

This might be a missing edge case.

[ArneBouillon]

I also encountered this issue today. When solving an ODEProblem and passing a saveat that contains tspan[1] and/or tspan[2], the following behaviour is observed:

• When tspan[1] is included in saveat and save_start = false, that time point is not contained in solve's solution (and so the solution has one time point less than does saveat.
• When tspan[2] is included in saveat and save_end = false, the terminal time point is still included in the solution.

An example is given by the following code.

function f(du, u, p, t)
    du[1] = -cos(u[1]) * u[1]
end
prob = ODEProblem(f, [10], (0.0, 0.4))

vcat(solve(prob, Tsit5(); saveat = 0:.1:.4).t...)
# =>  [0.0; 0.1; 0.2; 0.3; 0.4]
vcat(solve(prob, Tsit5(); saveat = 0:.1:.4, save_start = true, save_end = true).t...)
# =>  [0.0; 0.1; 0.2; 0.3; 0.4]
vcat(solve(prob, Tsit5(); saveat = 0:.1:.4, save_start = false, save_end = false).t...)
# =>       [0.1; 0.2; 0.3; 0.4]

This apparent inconsistency between how save_start and save_end are treated is what convinced me that this behaviour is probably a bug, instead of intended, although the differences in the documentation of the two options make it so that I'm not 100% sure. @ChrisRackauckas, could you potentially shed some light on that?

If it is concluded that this is indeed a bug, I would propose the following solution (and can take a stab at implementing it).

• Change the behaviour so that the last example also returns [0.0; 0.1; 0.2; 0.3; 0.4].
• Potentially tweak the documentation such that the docs of the save_start option more closely resemble those of save_end, which read "Denotes whether the final timepoint is forced to be saved, regardless of the other saving settings."

If the example above is the intended behaviour, it might be worthwhile to document this very explicitly (unless it already is and I missed it).

Thanks!

[oscardssmith]

this seems relatively bugish. IMO, the behavior should probably be that saveat is respected. if I were writing the API from scratch, I would probably delete save_start and save_end and just make the default saveat be the tspan, but that's obviously breaking to do now

[ChrisRackauckas]

That would be insufficient for many applications. It's running under the assumption that there's a an easy static representation for any array which just isn't true. If you try to represent the saveat=1e-1, save_first =false, save_end = true for tspan = (0.0,10000.0) you wouldn't be able to do it in what you're suggesting in a GPU kernel, not without creating a new array type at least.

Next, I explicitly add the start and end times to saveat, and specify that save_start and save_end are false. I expect the solution at four times values to be printed.

The booleans are explicit overrides. It's assumed that if you're going to set those booleans then it's going to be respected. IIRC that is the documented behavior, "regardless of the other saving settings", precisely for the reasons mentioned above.

When tspan[2] is included in saveat and save_end = false, the terminal time point is still included in the solution.

That one is a bug.

[oscardssmith]

wouldn't that example just be saveat=1e-1:1e-1:10000?

[ChrisRackauckas]

Nowadays mostly yes, but since ranges have interesting higher precision corrections to stepping points they are not compatible with most GPU kernels. On CUDA that now has a fallback to the inexact form, and that won't work on many platforms like Metal still.

[oscardssmith]

can't you collect the saveat points in advance? you need all of them in the solution object anyway.

[ChrisRackauckas]

Not if you want it to be both static and different between solves (can be different sized)

[ArneBouillon]

That one is a bug.

Sure, that's fair too, as long as it's consistent!

The booleans are explicit overrides. It's assumed that if you're going to set those booleans then it's going to be respected. IIRC that is the documented behavior, "regardless of the other saving settings", precisely for the reasons mentioned above.

Interesting, this is the exact opposite of how I interpreted that sentence 😅 When I read "Denotes whether the final timepoint is forced to be saved, regardless of the other saving settings.", I think "When this option is true, the final timepoint will definitely be saved; otherwise, it will depend on the other saving settings." Hence, when the boolean is false, I assume that the explicit override is off and the other saving settings are respected.

My confusion probably stems from the fact that it's a negative explicit override: the default is true, and a value of false forces an override. I think this is atypical enough to warrant more detailed documentation.

Note that save_start uses a different formulation: "Denotes whether the initial condition should be included in the solution type as the first timepoint." This doesn't imply an override in any direction, so expanding this documentation might also be worthwhile.

[ChrisRackauckas]

I'd be happy to take a PR to improve the docstring clarity. That might be the start here, and since I am the one who wrote it and knows what it should mean, I'm the worst person to edit it since now you know what is unclear about it 😅 . So please PR.

[ArneBouillon]

That one is a bug.

Do you want a PR that documents the current behaviour, or first to make the behaviour more consistent? And if so, in what direction?

[ArneBouillon]

I tried a few more examples in addition to the ones from my previous comment. Combined:

vcat(solve(prob, Tsit5(); saveat = 0:.1:.4).t...)
# => [0.0; 0.1; 0.2; 0.3; 0.4]
vcat(solve(prob, Tsit5(); saveat = 0:.1:.4, save_start = true, save_end = true).t...)
# => [0.0; 0.1; 0.2; 0.3; 0.4]
vcat(solve(prob, Tsit5(); saveat = 0:.1:.4, save_start = false, save_end = false).t...)
# =>      [0.1; 0.2; 0.3; 0.4]

vcat(solve(prob, Tsit5()).t...)
# => [0.0; 0.07; 0.16; 0.24; 0.34; 0.4]
vcat(solve(prob, Tsit5(); save_start = true, save_end = true).t...)
# => [0.0; 0.07; 0.16; 0.24; 0.34; 0.4]
vcat(solve(prob, Tsit5(); save_start = false, save_end = false).t...)
# =>      [0.07; 0.16; 0.24; 0.34; 0.4]

vcat(solve(prob, Tsit5(); saveat = [.2]).t...)
# =>      [0.2]
vcat(solve(prob, Tsit5(); saveat = [.2], save_start = true, save_end = true).t...)
# => [0.0; 0.2; 0.4]
vcat(solve(prob, Tsit5(); saveat = [.2], save_start = false, save_end = false).t...)
# =>      [0.2]

Replacing saveat = 0:.1:.4 by saveat = .1 yields the same results. Interestingly, due to the defaults listed below, saveat = .15 behaves differently from saveat = 0:.15:.4 when save_end is not given.

How I understand the current behaviour is as follows:

• save_start: Denotes whether the initial condition should be included in the solution type as the first timepoint. This setting overrides saveat when set to false. Defaults to save_everystep || isempty(saveat) || saveat isa Number || prob.tspan[1] in saveat.
• save_end: Denotes whether the final condition should be included in the solution type as the final timepoint. This setting is overridden by other saving settings when set to false. Defaults to save_everystep || isempty(saveat) || saveat isa Number || prob.tspan[2] in saveat.

I have created two pull requests (here and here) to add this wording to the documentation.

I do feel that these options are overly complicated in their behaviour and should currently be avoided by users...

[ArneBouillon]

Let's centralise discussion here.

This comment mentions that a bugfix is necessary, rather than a doc change. Should both arguments be handled as save_start currently is; that is, overriding saveat when explicitly given? If so, I imagine there are a lot of different solvers implementing the DiffEq.jl interface that will need to be changed (or, of course, that are currently inconsistent with OrdinaryDiffEq.jl).

[ChrisRackauckas]

We should create a central discussion on DifferentialEquations.jl for this.

Should both arguments be handled as save_start currently is; that is, overriding saveat when explicitly given?

Yes.

if so, I imagine there are a lot of different solvers implementing the DiffEq.jl interface that will need to be changed (or, of course, that are currently inconsistent with OrdinaryDiffEq.jl).

yes, we should make an issue with checkboxes so we do them all. DiffEqGPU, Sundials, ODEInterface, DASSL, LSODA, MATLABDiffEq, SciPyDiffEq, etc. off the top of my head. It's in the docs page. We can mark this as good first issue since #2095 shows how to do it (most of the packages have a similar savevalues! kind of thing going on), so it's rather straightforward to do and would probably be a good first issue for training new folks.

[ArneBouillon]

Looks good! Due to the current defaults, the following potentially unintuitive situation remains:

vcat(solve(prob, Tsit5(); saveat = 0:.15:.4).t...)
# => [0.0; 0.15; 0.3]
vcat(solve(prob, Tsit5(); saveat = .15).t...)
# => [0.0; 0.15; 0.3; 0.4]

I don't know if you consider this intended behaviour, or not.

[ChrisRackauckas]

Yeah it's a bit odd but it's intended. I would say what's odd is that 0:.15:.4 == [0.0; 0.15; 0.3], this is just a consistent reflection of that. The latter is what a lot of people actually want, which is why we cannot just rely on ranges. Maybe we can simplify this by creating some special ranges with some "end" handling behavior, but for now this system works well enough