Fix Rootfinding in Callbacks #350

kanav99 · 2019-10-13T19:43:10Z

Ref: SciML/DifferentialEquations.jl#516 (comment)

Fixes SciML/DifferentialEquations.jl#516 , Fixes SciML/DifferentialEquations.jl#551 , Fixes bifurcationkit/BifurcationKit.jl#9

ChrisRackauckas · 2019-10-14T02:18:33Z

src/callbacks.jl

@@ -503,6 +503,10 @@ function find_callback_time(integrator,callback::VectorContinuousCallback,counte
              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))


might as well put this in the loop

set the atol and rtol to 0

ChrisRackauckas · 2019-10-14T02:18:57Z

src/callbacks.jl

@@ -503,6 +503,10 @@ function find_callback_time(integrator,callback::VectorContinuousCallback,counte
              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))
+            sign_bottom_θ = sign(zero_func(bottom_θ))
+            while sign(zero_func(Θ)) != sign_bottom_θ


going 1 by 1 can be really slow. It would be better to just put the tolerance to like 0

What do you mean by tolerance? Just discard all the callback tolerances and put atol = 0 in find_zero call?

It's no good, still the same error.

how many times is it looping?

src/callbacks.jl

kanav99 · 2019-10-14T17:25:11Z

In the first callback event, the ball goes from +ve to -ve, so essentially zero_func(θ) should be positive as we prevfloat the root. But actually it's negetive. A single prevfloat is not able to pull θ back to positive. Issue is over here. But it's still a mystery why it runs fine on debugging.

Am I really right about this? Because the tests fail when I limit the loop. That means tests used to pass without the loop. It may also mean that our tests were not good enough.

ChrisRackauckas · 2019-10-15T09:43:18Z

src/callbacks.jl

@@ -435,7 +435,7 @@ function find_callback_time(integrator,callback::ContinuousCallback,counter)
            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))
+          Θ = prevfloat(find_zero(zero_func, (bottom_θ,top_Θ), Roots.AlefeldPotraShi(), atol = 0))


ChrisRackauckas · 2019-10-15T10:04:28Z

It looks like when you need zero tolerances, Roots.jl sidesteps the algorithm choice:

https://github.com/JuliaMath/Roots.jl/blob/2921a2dba2023de3f25681383009b24cde3370b2/src/bracketing.jl#L318-L324

https://github.com/JuliaMath/Roots.jl/blob/2cb8b737de06496ee45cd2b67b6eae1f57dd4532/src/find_zero.jl#L151-L192

Let's set the 4 tolerances to zero and remove the algorithm choice and see what we get. I would be surprised if this causes a major performance difference, since most of the condition calls are likely due to checking if there is a root, not the few times that a root is found. However, this should make it a lot more robust and guarantee a 1 floating point number difference from the zero.

kanav99 · 2019-10-15T14:10:17Z

I guess it worked. I will compare the ncondition and see if there is any noticable regression. And running OrdinaryDiffEq tests now.
Also, do we need the loop now?

ChrisRackauckas · 2019-10-15T14:24:17Z

I don't know if we need the loop, since that should make it so prevfloat is always one to the left. Maybe have it only possibly prevfloat once, and otherwise error? Since we will want feedback on how well this actually works, and that shouldn't increase the cost in any noticeable way to do one more condition call.

kanav99 · 2019-10-15T14:34:24Z

I just bumped a println(sol.destats.ncondition) after every solve call here.
On master:

On this branch:

Only the fourth call, bouncing ball example, becomes double in count.

ChrisRackauckas · 2019-10-15T14:45:05Z

It's double because the bouncing ball example is such a simple ODE that it literally has an event every single step. We saw that in this paper:

https://www.sciencedirect.com/science/article/abs/pii/S0965997818310251

I think it's safe to say that's an edge case, where if you literally have a root every single step then using the more robust method uses twice as many condition calls, where the condition is trivial 🤷‍♂ . We're probably fine.

kanav99 · 2019-10-15T14:54:59Z

Okay. For the looping, I guess single prevfloat should be enough, otherwise error. Glad this error was caught, but didn't expected this error to show up in such simple problem.

ChrisRackauckas · 2019-10-15T15:03:00Z

this is the beauty of having users who write test cases which break when floating point is exact. Being off by 1 floating point number is probably the reason for half of our event handling or adaptive time stepping tests 🤣.

Let's make sure we get this test case added as a downstream test.

ChrisRackauckas · 2019-10-15T15:47:31Z

Downstream tests are now failing.

kanav99 · 2019-10-15T16:28:27Z

That means the rootfinding still isn't accurate enough to control a bouncing ball 🤣 I guess we should have some loop then. 10 times is good enough. I think that we should have our own rootfinding that works on finding a number less than the actual root.

ChrisRackauckas · 2019-10-15T16:31:37Z

wait, how far off is it? The docs there literally say it should be "exact" 🤷‍♂

This reverts commit 8f1a03b.

kanav99 · 2019-10-15T16:52:16Z

After adding the loop, I added a @show prevfloat_idx and it is mostly 0 and just 1 once. Note that we do prevfloat once outside the loop too. Its actually pretty accurate but still it fails 😛.

kanav99 · 2019-10-15T17:40:01Z

The PSOS event test mostly gives 0 prevfloat_idx but many times it is a 4 digit integer.

ChrisRackauckas · 2019-10-16T10:10:46Z

Wait, with PSOS the bisection halts very early?

ChrisRackauckas · 2019-10-20T15:28:18Z

It should be "impossible" for it to be that far away with Bisection. That would be worth to build into an MWE.

ChrisRackauckas · 2020-02-03T21:29:30Z

@kanav99 where did we leave this?

kanav99 · 2020-02-04T07:58:48Z

The issue was that the result of root was not consistent in terms of sign. We need something which either gives us sign of the lower limit of exact Root or the upper limit. Roots.jl doesn't give me this consistency. I saw your mention and yeah, if we get an interval instead of exact value, that should solve our issue.

ChrisRackauckas · 2020-02-13T18:50:49Z

src/callbacks.jl

@@ -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.")


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

kanav99 · 2020-02-14T21:30:09Z

I guess scope of the PR is complete. Some AD issues need to be addressed. I will open a seperate issue for that probably.
Please squash before merging, history of this PR is messed up 😛

ChrisRackauckas · 2020-02-15T04:59:21Z

Let's get some more AD tests before committing to it. That's a pretty core feature that we can't break.

recursively run prevfloat in rootfind

605dac4

ChrisRackauckas reviewed Oct 14, 2019

View reviewed changes

limiting loop in rootfind

d7e766d

kanav99 commented Oct 14, 2019

View reviewed changes

src/callbacks.jl Show resolved Hide resolved

set atol to 0

2fdc0e0

ChrisRackauckas reviewed Oct 15, 2019

View reviewed changes

set rtol to 0

de3bb6e

ChrisRackauckas mentioned this pull request Oct 15, 2019

Finishing with bisection? JuliaMath/Roots.jl#165

Open

set rootfinding to best accuracy

da47cfa

kanav99 added 2 commits October 15, 2019 20:34

remove prevfloat loop

8f1a03b

add test for prevfloat working

08f0d44

ChrisRackauckas approved these changes Oct 15, 2019

View reviewed changes

Revert "remove prevfloat loop"

6033b86

This reverts commit 8f1a03b.

ChrisRackauckas mentioned this pull request Oct 20, 2019

More robust discontinuity tracking SciML/DelayDiffEq.jl#159

Merged

ChrisRackauckas mentioned this pull request Feb 2, 2020

Strange result of multiplication JuliaPhysics/Measurements.jl#64

Closed

rveltz mentioned this pull request Feb 13, 2020

Error with DifferentialEquations and multiple callbacks bifurcationkit/BifurcationKit.jl#9

Closed

Merge branch 'master' into kg/cbfix

47ca048

kanav99 force-pushed the kg/cbfix branch from dd78203 to 62f2253 Compare February 13, 2020 17:53

try fixing rootfinding

990416a

kanav99 force-pushed the kg/cbfix branch from 62f2253 to 990416a Compare February 13, 2020 17:55

change tolerance to 0

ffde953

ChrisRackauckas reviewed Feb 13, 2020

View reviewed changes

kanav99 added 4 commits February 14, 2020 01:06

Cleanup find_zero_bracket

6535e9e

fix typo and add example

7800631

Change algorithm for Rootfinding

8feac3f

hopefully fix the rootfinding

8da011f

ChrisRackauckas mentioned this pull request Feb 15, 2020

beef up the AD tests for ContinuousCallback SciML/OrdinaryDiffEq.jl#1034

Merged

ChrisRackauckas merged commit 29e8b76 into master Feb 15, 2020

ChrisRackauckas deleted the kg/cbfix branch February 15, 2020 15:02

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Fix Rootfinding in Callbacks #350

Fix Rootfinding in Callbacks #350

kanav99 commented Oct 13, 2019 •

edited by ChrisRackauckas

Loading

ChrisRackauckas Oct 14, 2019

ChrisRackauckas Oct 14, 2019

ChrisRackauckas Oct 14, 2019

kanav99 Oct 14, 2019

kanav99 Oct 14, 2019

ChrisRackauckas Oct 14, 2019

kanav99 Oct 14, 2019

kanav99 commented Oct 14, 2019

ChrisRackauckas Oct 15, 2019

ChrisRackauckas commented Oct 15, 2019

kanav99 commented Oct 15, 2019

ChrisRackauckas commented Oct 15, 2019

kanav99 commented Oct 15, 2019

ChrisRackauckas commented Oct 15, 2019

kanav99 commented Oct 15, 2019

ChrisRackauckas commented Oct 15, 2019

ChrisRackauckas commented Oct 15, 2019

kanav99 commented Oct 15, 2019

ChrisRackauckas commented Oct 15, 2019

kanav99 commented Oct 15, 2019

kanav99 commented Oct 15, 2019

ChrisRackauckas commented Oct 16, 2019

ChrisRackauckas commented Oct 20, 2019

ChrisRackauckas commented Feb 3, 2020

kanav99 commented Feb 4, 2020

ChrisRackauckas Feb 13, 2020

kanav99 commented Feb 14, 2020

ChrisRackauckas commented Feb 15, 2020

Fix Rootfinding in Callbacks #350

Fix Rootfinding in Callbacks #350

Conversation

kanav99 commented Oct 13, 2019 • edited by ChrisRackauckas Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kanav99 commented Oct 14, 2019

Choose a reason for hiding this comment

ChrisRackauckas commented Oct 15, 2019

kanav99 commented Oct 15, 2019

ChrisRackauckas commented Oct 15, 2019

kanav99 commented Oct 15, 2019

ChrisRackauckas commented Oct 15, 2019

kanav99 commented Oct 15, 2019

ChrisRackauckas commented Oct 15, 2019

ChrisRackauckas commented Oct 15, 2019

kanav99 commented Oct 15, 2019

ChrisRackauckas commented Oct 15, 2019

kanav99 commented Oct 15, 2019

kanav99 commented Oct 15, 2019

ChrisRackauckas commented Oct 16, 2019

ChrisRackauckas commented Oct 20, 2019

ChrisRackauckas commented Feb 3, 2020

kanav99 commented Feb 4, 2020

Choose a reason for hiding this comment

kanav99 commented Feb 14, 2020

ChrisRackauckas commented Feb 15, 2020

kanav99 commented Oct 13, 2019 •

edited by ChrisRackauckas

Loading