Implement automatic absolute tolerance computation as MATLAB

[ronisbr]

Hi guys!

Simulink, as can be seen here https://www.mathworks.com/help/simulink/gui/absolute-tolerance.html, has the capability of controlling the absolute tolerance automatically. It seems a very simple algorithm:

The default value (auto) initially sets the absolute tolerance for each state to 1e-6. As the simulation progresses, the absolute tolerance for each state is reset to the maximum value that the state has thus far assumed times the relative tolerance for that state.

For example, if a state goes from 0 to 1 and the Relative tolerance is 1e-3, then by the end of the simulation, the Absolute tolerance is set to 1e-3.

If the computed setting is not suitable, you can determine an appropriate setting yourself.

I am wondering if it is possible to implement this feature in this package.

[ChrisRackauckas]

Yes, the implementation would be pretty easy to do. The hard part would be the interface. We would need to think about what to use instead of Auto though. This could be done without overhead by making it be a different type, like you can set abstol=AutoTolerance (I don't like the name Auto) and the timestepping methods could use type-checking to set this behavior.

@mauro3 you have any input?

[mauro3]

This would be something each solver would have to implement, right? Or a callback? What about setting abstol=0 to trigger automatic setting (or an error if not supported)?

[ChrisRackauckas]

This would be something each solver would have to implement, right?

Yeah. It would be possible with a callback, but then we'd have to keep expanding what's in the standard callbacks.

What about setting abstol=0 to trigger automatic setting (or an error if not supported)?

abstol = 0 has a already has a meaning though. It means "don't use absolute tolerance, just relative tolerance". Maybe we could use a negative number, or Inf?

Side note: we really need that keyword arg dispatch PR to be part of v0.6.

JuliaLang/julia#16580

[ChrisRackauckas]

Thinking about this again: how much does this actually help? I am torn between "if MATLAB has it, it must be useful enough that we shouldn't require users who want it to write a very special callback" and "is it actually useful? Won't the relative tolerance kind of handle this?"

[mauro3]

Maybe a library of standard call-backs would be good?

[ChrisRackauckas]

Maybe a library of standard call-backs would be good?

I don't think so, unless it's one solver specific (OrdinaryDiffEq.jl would be the natural target since it's the most flexible to do these things). I will be creating libraries which generate event functions because events are a high enough level that they can be generically applied to many solvers. Callbacks are much too low level. This problem is actually a good example of that. Let's say for example, to accommodate this problem, we extend the common callback proposal to be:

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

I claim this is a bad idea. For OrdinaryDiffEq.jl, reltol and abstol are immutables. In the near future I'd like to get them as arrays, but for now they are numbers. This means you cannot mutate them in the callback. So okay, we can handle this by instead passing a Ref{typeof(abstol)} instead of the numbers themselves, and then in the callback mutate them with abstol[]. Good that works for OrdinaryDiffEq.jl.

But that will work for nothing else. It won't work for PR49 because it uses an immutable type for its adaptive options (and its probably unwise to change that, so maybe you make reltol and abstol as Ref inside of AdaptiveOptions and so you pass this to the callback?). If you mutate the values for Sundials (or any other interop solver), it won't actually affect the solver until you call

flag = @checkflag CVodeSStolerances(mem, reltol, abstol)

So now the code starts to get really messy. You could create a macro

function callback(alg::Alg,t,u,du,f,cache,reltol,abstol,args...)
  @adaptive_abstol
end

which generates

function callback(alg::Alg,t,u,du,f,cache,reltol,abstol,args...)
  # calculate new abstol
  if typeof(alg) <: SundialsAlgorithm
    ...
  ...
end

and condition that all out, but you're still missing the fact that you need to cache the value for the maximum u over the timeseries, and you have to deal with the fact that each callback will have a different args... (and what's worse is we have to be using those args... here since that's where mem would be). So what you'd actually have to do is have to do is have a macro which generates the full callback function

@ode_callback begin
  @adaptive_abstol
end

and have it create a separate dispatch with "standard naming" for each algorithm (what repo would this be in? It would need to depend on every solver package), and the inner macro should still have the conditionals unless you make that outer macro "super smart" which would be very time consuming to develop / maintain.

The other way to do this with Sundials is to allow it to accept abstol and reltol Ref which can be mutated, and then after the callback check for mutation and make it automatically handle the mem call. If you make every solver have to implement this check, then it's a lot of work for the potential that someone wants to do this, and probably more work than just implementing this behavior in every solver.

I think this shows how callbacks are low-level and close enough to the actual algorithm that it would be a nightmare to maintain. But callbacks are very useful. When I give you the function signature for OrdinaryDiffEq.jl, it's pretty easy to write a callback which implements this for OrdinaryDiffEq.jl. And same for Sundials.jl (if you're familiar with the Sundials API). So callbacks let you easily extend an algorithm as much as you want, but uniformity in doing so is out of the question (unlike events).

[Maybe I should just write a library of really cool callbacks for OrdinaryDiffEq.jl and not worry about the other solvers, just for handling these kinds of odd cases. It would be restricting, but for a case like this its an easy fix]

The real problem here is that it's hard to do via events as well. Sure, we could pass Refs to reltol and abstol in apply_event!, but again, for cases like Sundials you'd have to do a bunch of "check if Ref changed" calculations as mentioned earlier (and have a running counter for the maximum value of u. Forgot about that earlier...).

[ChrisRackauckas]

Scratch what I had at the end there. For interop packages like Sundials, it will need to check for changes and use API calls to change the temporary arrays in the mem for events to u, so why not have it do the same for the tolerances? These could conceivably be added to the events as mutable references then. A running counter for the maximum value of u would be a field in the event_cache (or make apply_event! be a call overloaded type with a field which holds that value!).

So that should work everywhere except PR49 (unless it finds a way to deal with the immutability issue). This could then go into events library.

That would extend it to:

apply_event!(t,u,du,event_cache,cache,abstol,reltol,funcs...)

Is it fine to add more arguments? Will there be other cases of tolerance changing?

[mauro3]

I think this would be better combined with the init stuff:

whatever_init_returns = init(p, alg; solve_kwargs...)
...
apply_event!(t,u,du,event_cache,cache,whatever_init_returns,...)

The whatever_init_returns would contain abstol and all other used kwargs.

[ChrisRackauckas]

Wow, sorry this got so long. But here's the brain dump. There's a lot of very integrated things in this idea to mention.

---

That doesn't help in all cases, and would actually exacerbate the problem in some. For any solver which cannot work directly on the init type (any interop solver), every field which needs some nontrivial change when changed (i.e. any field which needs to be passed to the Sundials/DASKR/LSODA/... mem core) would have to have extra checks.

That said, it's worth looking into. Now that we really have something together I think we can take stock and use our current results as a point to test off of. Specifically:

1. I would like to test in OrdinaryDiffEq.jl what happens if I don't @unpack the kwargs. If you change the internal ODEIntegrator to a mutable type, and directly use its fields for the options, do you take a hit in the benchmarks? We have a full set of benchmarks, and if they convince me that the mutable type is not hit often enough for this to matter performance-wise (for options), then it makes sense to start thinking of an options type.

(Do you want to run some benchmarks, or would you be willing to wait a little bit for me to get around to it?).

2. I think the lesson we learned is that the kwargs (the common solver options) should be restricted to the common set, algorithm/package-specific options should go with the algorithm, which could then be documented in the fashion of Extended algorithm documentation in constituent packages #120 . So that would fix the issue I had before of "flexible options", and would help enforcing uniformity.

One thing to note is that the options type would need to silently contain a field default_set which is used for setting the default algorithm to avoid infinite recursion. You can see why here:

#120

That could be a set DEOptions type. Its constructor would need to be DEOptions(prob,fields...) since it needs the type information from the problem. This would allow the pass in apply_event!, and then it would be up to the individual solvers to implement its modifications.

Once we have that, I think we should go a little further with init. We should then setup the solvers to mutate the solution type, so that

sol,integrator,opts = init(prob, alg; solve_kwargs...,pre_allocated_arrays...)
solve!(sol,prob,integrator,opts)

will mutate sol. We could just fold this all into one integrator type, but it's probably easier for users to just see this an know that they can mutate prob.tspan then solve! further, rather than having to setup a new interface on integrator such that integrator.prob always exists for mutation.

Of course,

sol = solve(prob,alg,pre_allocated...;kwargs...)

would just fold into the above.

Then, actually we can just put another type parameter into the solution object and just make it have the field integrator for this, making it:

sol,opts = init(prob, alg; solve_kwargs...,pre_allocated_arrays...)
solve!(sol,prob,opts)

---

There are two more questions that this opens up. The first is why keep prob and opts separate from sol here. prob actually is in the solution type, but I think it might make sense to think of this as "solve the new problem", and internally set sol.prob=prob to make sol.prob be the last problem solved. Or users could have to modify sol.prob, sol.opts. The issue I do have there is that makes it look like those options "truly were the case", where, with events, the options will have cases (like this issue) where they are changing throughout the solution, making me weary of putting it in the solution type. Then again, if it's "values at the end" for both, maybe that's a good enough convention.

Now, the last thing that probably will get brought up here is why keep integrator separate from alg. This seems like it's going back to the original suggestion of keeping alg as a type and not its instantiation, and then using Alg(prob,specific_args...). However, the reason why I am weary there is because Alg is a giant wad of not pre-computable type information. One way this is the case is for setting symbols like linear_solve in Sundials. I can see boolean arguments to turn on/off difference pieces of the solver also being an option in the near future, so the algorithm type can get pretty complicated type parameters.

I feel the next proposal would be to then just say, oh well solve will need dynamic dispatch but you can get around that by always using the init. However, it isn't true that you can always use the init as stated above. One place where that's the case is in parameter estimation algorithms which use autodifferentiation. These will actually need to change the type of the problem, and thus of the solution arrays, and thus one cannot pass the inited sol to a function which performs parameter estimation. But the biggest reason to init would be to pre-allocate, so I don't think it's a good idea to go away from pre-allocating to get rid of the integrate type. That said, this keeps the alg types simple, and lets the integrate type proliferate with the tableaus and everything.

So in the end, I think it's best to pass an instantiated Alg which creates a type with values and all of that, so that way this can be passed to solve, init, parameter estimation functions, etc. Alg is a documented way to generate some algorithm type, init would be a new initialization interface, integrator would be a new field of the solution type which is just documented and typed as "package does whatever the hell it wants with this", with the open question of putting the opts and prob in sol for solve!.

[ChrisRackauckas]

For an update on this issue, see SciML/Roadmap#10 (comment)

[ChrisRackauckas]

This has been implemented with the new callback interface in DiffEq callbacks.

https://github.com/JuliaDiffEq/DiffEqCallbacks.jl

[ronisbr]

Hi @ChrisRackauckas,

Sorry, I was completely offline due to a travel. I read your new Callback interface and the Integrator interface and I only can say what an amazing work! This solved all my problems :) I will study it better and port my code to use these new features.

Thanks you very much!
Ronan

[ChrisRackauckas]

No worries. You have time because it hasn't been released yet (hence the "developer preview" and the lack of docs change for the callbacks so far). If it does get released and you haven't changed your code, just Pkg.pin("OrdinaryDiffEq"), Pkg.pin("DiffEqBase") to the version you need. But I'll have an announcement.

[ronisbr]

Good! So I will clone master and test it!

[ronisbr]

Btw, can you point me what packages should I track from master?

[ChrisRackauckas]

OrdinaryDiffEq.jl, DiffEqBase.jl, and I think you may need ParameterizedFunctions.jl.

[ronisbr]

Hi @ChrisRackauckas,

Does this package support a vector as the absolute tolerance? I think this can only be fully implemented with this feature. The code you proposed at DiffEqCallbacks seems to compute the maximum value of the state vector, and then it multiplies by the relative tolerance to obtain the new absolute tolerance.

However, it seems that Matlab adapts the absolute tolerance for each value separately. Hence, we should have a vector of absolute tolerances with the dimension equals to that of the state vector. Am I right?

[ChrisRackauckas]

Yes. I just don't support vector tolerances yet because it will be free in v0.6. What I mean is that the following code:

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

currently is a loop. Without broadcast, I would need at least two different loops: one where the tolerances are indexed and one where they're not (and every combination if you want to support more!). However, instead of looping this could be broadcast:

atmp .= ((utilde.-u)./(integrator.opts.abstol.+max.(abs.(uprev),abs.(u)).*integrator.opts.reltol))

On v0.5, that will allocate more memory and slow things down. On v0.6 that will be just as efficient as the loop, and will handle both scalar and vector tolerances, along with every combination of them. So that's why I am waiting. Then, the same could happen to the callback in DiffEqCallbacks.