Callback bug fixes

[pengwyn]

In using VectorContinuousCallback in my code I came across several bugs, leading to exceptions, which I think I've fixed here.

1. The check integrator.vector_event_last_time == event_idx compared an Int to a Vector{Int} so it would always be false. I think it should be idx instead of event_idx.
2. The diff_t used in the nudging of the bottom t to find a different sign for bottom and top was using machine epsilon of the type, not the actual time. In my particular case, t was very large causing this to fail.
3. The check between find_callback_time and bisection was inconsistent. The nudger in find_callback_time was happy with a zero value for bottom_t, but bisection wanted a non-zero value.

I wanted to write some tests for this but ran out of time and figured that these are simple enough fixes. All the existing tests of DiffEqBase and OrdinaryDiffEq pass on my machine.

[kanav99]

This is probably one of the issues I was looking for..

[ChrisRackauckas]

yeah, this is one of those things where it's like... huh, yeah... that's umm... wow how did anything ever work before? 😆

[kanav99]

@ChrisRackauckas I am not good with eps, but should this be okay? I think this would be okay

[ChrisRackauckas]

yes, this is much better. I didn't catch this before, but you almost never want to do eps(typeof(t)) and always eps(t) since floating point is relative.

[kanav99]

@pengwyn thank you so much for your contribution! I agree with all 3 of the points you made.

All the existing tests of DiffEqBase and OrdinaryDiffEq pass on my machine.

Awesome! I assume your code works well with these changes too.

[devmotion]

Maybe computing zero(sign_top) before the loop would be slightly more efficient.

[pengwyn]

Shouldn't this be inlined by the compiler? Not that it is hard to move it outside the loop, but it is slightly less readable. Actually to contradict myself, I ran @code_warntype on it and it couldn't infer the result of zero_func.

Originally I thought that this duplicated check (between bisection and find_callback_time) would be best represented once as a function like:

different_sign(x::T,y::T) = sign(x) * sign(y) >= zero(T)

However, this wouldn't work in the scenario of Unitful quantities (there's a bit of asymmetry there between zero(x), one(x) and oneunit(x)). But this would work:

different_sign(x::T,y::T) = sign(x) * sign(y) >= zero(sign(x))

as long as zero(sign(x)) is inlined and the sign(x) calculation elided.

Alternatively, it could be argued that sign(x) is always one of [0, +1, -1] and so we could put directly:

different_sign(x::T,y::T) = sign(x) * sign(y) >= 0

[pengwyn]

I just realised that was a brain dump and not a proper reply. My question is should I a) keep the same style (and move the zero(..) out of the for loop)? or b) turn it into a separate function and put that into both find_callback_time and bisection?

[ChrisRackauckas]

Actually to contradict myself, I ran @code_warntype on it and it couldn't infer the result of zero_func.

Oh no, how come?

[pengwyn]

It turned out to be 3 bugs, one in my external code and two more in the callbacks.jl file. I'll open another pull request about them.

[devmotion]

Same here, maybe one could move zero(sign_top) out of the loop.

[ChrisRackauckas]

IIRC LLVM usually will handle this case if it's type-stable.

[pengwyn]

@pengwyn thank you so much for your contribution! I agree with all 3 of the points you made.

@kanav99 I'm glad you find it useful! I hope it helps with #564 too.

Awesome! I assume your code works well with these changes too.

Yeap, it's hitting the right benchmarks in my physics problem on everything I've tried so far!

[ChrisRackauckas]

We should check how many of https://github.com/SciML/DifferentialEquations.jl/labels/callbacks are fixed by this, and get tests on them. I assume that the fix of making the eps(t) instead of eps(typeof(t)) would fix more than a couple of them.