Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Change controller handling with the adaptive option #2284

Open
wants to merge 9 commits into
base: master
Choose a base branch
from
6 changes: 4 additions & 2 deletions lib/OrdinaryDiffEqDefault/src/OrdinaryDiffEqDefault.jl
Original file line number Diff line number Diff line change
@@ -1,8 +1,10 @@
module OrdinaryDiffEqDefault

using OrdinaryDiffEq: Vern7, Vern8, Vern9, Vern6, Tsit5, Rosenbrock23, Rodas5P, FBDF,
alg_stability_size, beta2_default, beta1_default, AutoSwitchCache, ODEIntegrator,
CompositeAlgorithm, OrdinaryDiffEqAlgorithm, OrdinaryDiffEqMutableCache, AutoAlgSwitch
alg_stability_size, beta2_default, beta1_default, AutoSwitchCache,
ODEIntegrator,
CompositeAlgorithm, OrdinaryDiffEqAlgorithm,
OrdinaryDiffEqMutableCache, AutoAlgSwitch
import OrdinaryDiffEq: is_mass_matrix_alg, default_autoswitch
import LinearSolve
using LinearAlgebra: I, isdiag
Expand Down
4 changes: 2 additions & 2 deletions lib/OrdinaryDiffEqDefault/src/default_alg.jl
Original file line number Diff line number Diff line change
Expand Up @@ -116,6 +116,6 @@ end

# hack for the default alg
function is_mass_matrix_alg(alg::CompositeAlgorithm{
<:Any, <:Tuple{Tsit5, Vern7, Rosenbrock23, Rodas5P, FBDF, FBDF}})
<:Any, <:Tuple{Tsit5, Vern7, Rosenbrock23, Rodas5P, FBDF, FBDF}})
true
end
end
2 changes: 1 addition & 1 deletion lib/OrdinaryDiffEqVerner/src/alg_utils.jl
Original file line number Diff line number Diff line change
Expand Up @@ -10,4 +10,4 @@ alg_order(alg::Vern9) = 9
alg_stability_size(alg::Vern6) = 4.8553
alg_stability_size(alg::Vern7) = 4.6400
alg_stability_size(alg::Vern8) = 5.8641
alg_stability_size(alg::Vern9) = 4.4762
alg_stability_size(alg::Vern9) = 4.4762
2 changes: 1 addition & 1 deletion lib/OrdinaryDiffEqVerner/src/algorithms.jl
Original file line number Diff line number Diff line change
Expand Up @@ -116,4 +116,4 @@ end
AutoVern6(alg; lazy = true, kwargs...) = AutoAlgSwitch(Vern6(lazy = lazy), alg; kwargs...)
AutoVern7(alg; lazy = true, kwargs...) = AutoAlgSwitch(Vern7(lazy = lazy), alg; kwargs...)
AutoVern8(alg; lazy = true, kwargs...) = AutoAlgSwitch(Vern8(lazy = lazy), alg; kwargs...)
AutoVern9(alg; lazy = true, kwargs...) = AutoAlgSwitch(Vern9(lazy = lazy), alg; kwargs...)
AutoVern9(alg; lazy = true, kwargs...) = AutoAlgSwitch(Vern9(lazy = lazy), alg; kwargs...)
2 changes: 1 addition & 1 deletion lib/OrdinaryDiffEqVerner/src/controllers.jl
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why does this file exist at all?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It looks like this was a mistake in the big refactors, we should delete it.

Original file line number Diff line number Diff line change
Expand Up @@ -5,4 +5,4 @@ end

@inline function accept_step_controller(integrator, controller::PIDController)
return integrator.qold >= controller.accept_safety
end
end
2 changes: 1 addition & 1 deletion lib/OrdinaryDiffEqVerner/src/interp_func.jl
Original file line number Diff line number Diff line change
Expand Up @@ -28,4 +28,4 @@ function DiffEqBase.interp_summary(::Type{cacheType},
Union{Vern9Cache, Vern9ConstantCache
}}
dense ? "specialized 9th order lazy interpolation" : "1st order linear"
end
end
14 changes: 7 additions & 7 deletions lib/OrdinaryDiffEqVerner/src/interpolants.jl
Original file line number Diff line number Diff line change
Expand Up @@ -6,15 +6,15 @@ RK_WITH_SPECIAL_INTERPOLATIONS = Union{
}

function _ode_interpolant(Θ, dt, y₀, y₁, k,
cache::RK_WITH_SPECIAL_INTERPOLATIONS,
idxs, T::Type{Val{D}}, differential_vars) where {D}
throw(DerivativeOrderNotPossibleError())
cache::RK_WITH_SPECIAL_INTERPOLATIONS,
idxs, T::Type{Val{D}}, differential_vars) where {D}
throw(DerivativeOrderNotPossibleError())
end

function _ode_interpolant!(out, Θ, dt, y₀, y₁, k,
cache::RK_WITH_SPECIAL_INTERPOLATIONS,
idxs, T::Type{Val{D}}, differential_vars) where {D}
throw(DerivativeOrderNotPossibleError())
cache::RK_WITH_SPECIAL_INTERPOLATIONS,
idxs, T::Type{Val{D}}, differential_vars) where {D}
throw(DerivativeOrderNotPossibleError())
end

## Vern6
Expand Down Expand Up @@ -1040,4 +1040,4 @@ end
k[18][idxs] * b24Θdiff + k[19][idxs] * b25Θdiff +
k[20][idxs] * b26Θdiff) * invdt3
out
end
end
33 changes: 28 additions & 5 deletions src/integrators/controllers.jl
Original file line number Diff line number Diff line change
@@ -1,6 +1,10 @@
abstract type AbstractController end
using OrdinaryDiffEq

@inline function accept_step_controller(integrator, ::AbstractController)
return integrator.EEst <= 1
end

@inline function stepsize_controller!(integrator, alg)
stepsize_controller!(integrator, integrator.opts.controller, alg)
end
Expand All @@ -23,6 +27,30 @@ reset_alg_dependent_opts!(controller::AbstractController, alg1, alg2) = nothing

DiffEqBase.reinit!(integrator::ODEIntegrator, controller::AbstractController) = nothing

"""
NonAdaptiveController()

This Controller exists to match the interface when one does not want to use a controller,
basically if you want to keep a fixed time step.
"""
struct NonAdaptiveController <: AbstractController
end

@inline function stepsize_controller!(integrator, ::NonAdaptiveController, alg)
nothing
end

@inline function accept_step_controller(integrator, ::NonAdaptiveController)
return true
end

function step_accept_controller!(integrator, ::NonAdaptiveController, alg, q)
integrator.dt
end

function step_reject_controller!(integrator, ::NonAdaptiveController, alg)
end

# Standard integral (I) step size controller
"""
IController()
Expand Down Expand Up @@ -308,11 +336,6 @@ end
return dt_factor
end

# checks whether the controller should accept a step based on the error estimate
@inline function accept_step_controller(integrator, controller::AbstractController)
return integrator.EEst <= 1
end
ranocha marked this conversation as resolved.
Show resolved Hide resolved

@inline function accept_step_controller(integrator, controller::PIDController)
return integrator.qold >= controller.accept_safety
end
Expand Down
20 changes: 1 addition & 19 deletions src/integrators/integrator_utils.jl
Original file line number Diff line number Diff line change
Expand Up @@ -222,7 +222,7 @@ function _loopfooter!(integrator)
end
integrator.last_stepfail = true
integrator.accept_step = false
elseif integrator.opts.adaptive
else
q = stepsize_controller!(integrator, integrator.alg)
integrator.isout = integrator.opts.isoutofdomain(integrator.u, integrator.p, ttmp)
integrator.accept_step = (!integrator.isout &&
Expand Down Expand Up @@ -255,24 +255,6 @@ function _loopfooter!(integrator)
else # Reject
integrator.stats.nreject += 1
end
elseif !integrator.opts.adaptive #Not adaptive
integrator.stats.naccept += 1
integrator.tprev = integrator.t
integrator.t = if has_tstop(integrator)
tstop = integrator.tdir * first_tstop(integrator)
if abs(ttmp - tstop) <
100eps(float(integrator.t / oneunit(integrator.t))) * oneunit(integrator.t)
tstop
else
ttmp
end
else
ttmp
end
integrator.last_stepfail = false
integrator.accept_step = true
integrator.dtpropose = integrator.dt
handle_callbacks!(integrator)
end
if integrator.opts.progress && integrator.iter % integrator.opts.progress_steps == 0
t1, t2 = integrator.sol.prob.tspan
Expand Down
6 changes: 5 additions & 1 deletion src/solve.jl
Original file line number Diff line number Diff line change
Expand Up @@ -364,7 +364,11 @@ function DiffEqBase.__init(
end

if controller === nothing
controller = default_controller(_alg, cache, qoldinit, beta1, beta2)
if adaptive
controller = default_controller(_alg, cache, qoldinit, beta1, beta2)
else
controller = NonAdaptiveController()
end
Comment on lines +366 to +370
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This will introduce an additional type instability since the type of controller is determined by the value of adaptive. I cannot really judge how severe this issue is right now. Others probably have more to say about this.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is there a type instability ?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

julia> function foo(; adaptive = true)
           if adaptive
               return 1.0
           else
               return 1
           end
       end
foo (generic function with 1 method)

julia> @code_warntype foo(adaptive = true)
MethodInstance for Core.kwcall(::@NamedTuple{adaptive::Bool}, ::typeof(foo))
  from kwcall(::NamedTuple, ::typeof(foo)) @ Main REPL[1]:1
Arguments
  _::Core.Const(Core.kwcall)
  @_2::@NamedTuple{adaptive::Bool}
  @_3::Core.Const(foo)
Locals
  adaptive::Union{}
  @_5::Bool
Body::Union{Float64, Int64}
1 ─       Core.NewvarNode(:(adaptive))
│         Core.NewvarNode(:(@_5))
│   %3  = Core.isdefined(@_2, :adaptive)::Core.Const(true)
└──       goto #3 if not %3
2 ─       (@_5 = Core.getfield(@_2, :adaptive))
└──       goto #4
3 ─       Core.Const(:(@_5 = true))
4%8  = @_5::Bool%9  = (:adaptive,)::Core.Const((:adaptive,))
│   %10 = Core.apply_type(Core.NamedTuple, %9)::Core.Const(NamedTuple{(:adaptive,)})
│   %11 = Base.structdiff(@_2, %10)::Core.Const(NamedTuple())
│   %12 = Base.pairs(%11)::Core.Const(Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}())
│   %13 = Base.isempty(%12)::Core.Const(true)
└──       goto #6 if not %13
5 ─       goto #7
6 ─       Core.Const(:(Base.kwerr(@_2, @_3)))
7%17 = Main.:(var"#foo#1")(%8, @_3)::Union{Float64, Int64}
└──       return %17

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Note the Union{Float64, Int64}

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, but is there no type instability due to the fact that the type of the controller depends also on the algorithm is adaptive is true ?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ChrisRackauckas @oscardssmith What do you think about this ?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

it is a minor problem. Are there settings of the normal PI controllers that will make them not adapt ever? It feels like there should be some combination of 0 and Inf that make it never adapt.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I do not know, is it related to my changes ?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A better strategy would be to add traits to controllers for overrides_adaptive(controller) or something, which would allow the default controller to stay the same (making it better for compilation) while allowing the alternative controllers to bypass the if adaptive checks by making them if adaptive || overrides_adaptive(controller).

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ChrisRackauckas I thought over this but for each things I tried, I do not succeed to prevent from a type instability. Can you explain a little more ? I think I have missunderstood something but I do not know what.

end

save_end_user = save_end
Expand Down
Loading