Skip to content

Fix Rootfinding in Callbacks #350

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.

Sign up for GitHub

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

Merged
merged 15 commits into from
Feb 15, 2020
Merged

Fix Rootfinding in Callbacks #350

Show file tree
Hide file tree
Changes from 11 commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
605dac4
recursively run prevfloat in rootfind
kanav99 Oct 13, 2019
d7e766d
limiting loop in rootfind
kanav99 Oct 14, 2019
2fdc0e0
set atol to 0
kanav99 Oct 15, 2019
de3bb6e
set rtol to 0
kanav99 Oct 15, 2019
da47cfa
set rootfinding to best accuracy
kanav99 Oct 15, 2019
8f1a03b
remove prevfloat loop
kanav99 Oct 15, 2019
08f0d44
add test for prevfloat working
kanav99 Oct 15, 2019
6033b86
Revert "remove prevfloat loop"
kanav99 Oct 15, 2019
47ca048
Merge branch 'master' into kg/cbfix
kanav99 Feb 13, 2020
990416a
try fixing rootfinding
kanav99 Feb 13, 2020
ffde953
change tolerance to 0
kanav99 Feb 13, 2020
6535e9e
Cleanup find_zero_bracket
kanav99 Feb 13, 2020
7800631
fix typo and add example
kanav99 Feb 13, 2020
8feac3f
Change algorithm for Rootfinding
kanav99 Feb 13, 2020
8da011f
hopefully fix the rootfinding
kanav99 Feb 14, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
31 changes: 29 additions & 2 deletions src/callbacks.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -522,6 +522,32 @@ function findall_events(affect!,affect_neg!,prev_sign,next_sign)
findall(x-> ((prev_sign[x] < 0 && affect! !== nothing) || (prev_sign[x] > 0 && affect_neg! !== nothing)) && prev_sign[x]*next_sign[x]<=0, keys(prev_sign))
end

function find_zero_bracket(fs, x0; kwargs...)

x = Roots.adjust_bracket(x0)
T = eltype(x[1])
F = Roots.callable_function(fs)
method = Roots.Bisection()
state = Roots.init_state(method, F, x)
options = Roots.init_options(method, state; kwargs...)

# check if tolerances are exactly 0
# iszero_tol = iszero(options.xabstol) && iszero(options.xreltol) && iszero(options.abstol) && iszero(options.reltol)

# if iszero_tol
# if T <: FloatNN
# return Roots.find_zero(F, x, Roots.BisectionExact(); kwargs...)
# else
# return find_zero(F, x, Roots.A42(); kwargs...)
# end
# end

find_zero(method, F, options, state, Roots.NullTracks())

state.xn0, state.xn1

end

function find_callback_time(integrator,callback::ContinuousCallback,counter)
event_occurred,interp_index,Θs,prev_sign,prev_sign_index,event_idx = determine_event_occurance(integrator,callback,counter)
if event_occurred
Expand Down Expand Up @@ -555,10 +581,11 @@ function find_callback_time(integrator,callback::ContinuousCallback,counter)
iter = 1
while sign(zero_func(bottom_θ)) == sign_top && iter < 12
bottom_θ *= 5
iter += 1
end
iter == 12 && error("Double callback crossing floating pointer reducer errored. Report this issue.")
Copy link
Member

Choose a reason for hiding this comment

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

I wonder if we still need to do all of this, or if we just immediately throw an error if we see this.

end
Θ = prevfloat(find_zero(zero_func, (bottom_θ,top_Θ), Roots.AlefeldPotraShi(), atol = callback.abstol/100))
Θ, _ = find_zero_bracket(zero_func, (bottom_θ,top_Θ), atol = 0, rtol = 0, xatol = 0, xrtol = 0)
integrator.last_event_error = ODE_DEFAULT_NORM(zero_func(Θ),integrator.t+integrator.dt*Θ)
end
#Θ = prevfloat(...)
Expand Down Expand Up @@ -626,7 +653,7 @@ function find_callback_time(integrator,callback::VectorContinuousCallback,counte
end
iter == 12 && error("Double callback crossing floating pointer reducer errored. Report this issue.")
end
Θ = prevfloat(find_zero(zero_func, (bottom_θ,top_Θ), Roots.AlefeldPotraShi(), atol = callback.abstol/100))
Θ, _ = find_zero_bracket(zero_func, (bottom_θ,top_Θ), atol = 0, rtol = 0, xatol = 0, xrtol = 0)
if Θ < minΘ
integrator.last_event_error = ODE_DEFAULT_NORM(zero_func(Θ),integrator.t+integrator.dt*Θ)
end
Expand Down
22 changes: 22 additions & 0 deletions test/downstream/event_detection_tests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,7 @@ function condition(u,t,integrator) # Event when event_f(u,t) == 0
u[1]
end
function affect!(integrator)
@test integrator.u[1] >= 0
integrator.u[2] = -integrator.u[2]
end
cb2 = ContinuousCallback(condition,affect!)
Expand All @@ -55,3 +56,24 @@ p = 9.8
prob = ODEProblem(f,u0,tspan,p)
sol = solve(prob,Tsit5(),callback=cb2)
@test minimum(sol') > -40

function vcondition!(out,u,t,integrator)
out[1] = u[1]
out[2] = u[2]
end

function vaffect!(integrator, event_idx)
@test integrator.u[1] >= 0.0
if event_idx == 1
integrator.u[2] = -integrator.u[2]
else
integrator.p = 0.0
end
end

u0 = [50.0,0.0]
tspan = (0.0,15.0)
p = 9.8
prob = ODEProblem(f,u0,tspan,p)
Vcb = VectorContinuousCallback(vcondition!,vaffect!, 2 , save_positions=(true,true))
sol = solve(prob,Tsit5(), callback=Vcb)