In integrating a system with a very large number of callbacks (method similar to https://pmc.ncbi.nlm.nih.gov/articles/PMC10470263/ for a large system of biological reactions), a very small proportion of the time solve() will fail with a "Double callback crossing floating pointer reducer errored".
Either the event is skipped, or the event fires, without terminating the program.
I suspect the fix is to delay updating integrator.last_event_error until after all event indices have been processed.
This is difficult because of the size of the system, and because it is highly floating-point error dependent. However, on my system the following, very artificial example reproduces the bug.
ERROR: LoadError: Double callback crossing floating pointer reducer errored. Report this issue.
Stacktrace:
[1] error(s::String)
@ Base ./error.jl:44
[2] find_callback_time(integrator::OrdinaryDiffEqCore.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, false, Vector{Float64}, Nothing, Float64, SciMLBase.NullParameters, Float64, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, Nothing, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, ODEFunction{false, SciMLBase.FullSpecialize, typeof(ode), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, OrdinaryDiffEqCore.InterpolationData{ODEFunction{false, SciMLBase.FullSpecialize, typeof(ode), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, Nothing, OrdinaryDiffEqTsit5.Tsit5ConstantCache, Nothing}, SciMLBase.DEStats, Nothing, Nothing, Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, typeof(ode), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, OrdinaryDiffEqTsit5.Tsit5ConstantCache, OrdinaryDiffEqCore.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{VectorContinuousCallback{typeof(condition), typeof(jump!), typeof(jump!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT), Float64, Int64, Rational{Int64}, Nothing, Nothing, Int64, Tuple{}}}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, Vector{Float64}, Float64, DiffEqBase.CallbackCache{Vector{Float64}, Vector{Float64}}, DiffEqBase.DefaultInit, Nothing}, callback::VectorContinuousCallback{typeof(condition), typeof(jump!), typeof(jump!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT), Float64, Int64, Rational{Int64}, Nothing, Nothing, Int64, Tuple{}}, counter::Int64)
@ DiffEqBase ~/.julia/packages/DiffEqBase/6sydR/src/callbacks.jl:505
[3] macro expansion
@ ~/.julia/packages/DiffEqBase/6sydR/src/callbacks.jl:130 [inlined]
[4] find_first_continuous_callback
@ ~/.julia/packages/DiffEqBase/6sydR/src/callbacks.jl:125 [inlined]
[5] find_first_continuous_callback
@ ~/.julia/packages/DiffEqBase/6sydR/src/callbacks.jl:123 [inlined]
[6] handle_callbacks!
@ ~/.julia/packages/OrdinaryDiffEqCore/GMkz9/src/integrators/integrator_utils.jl:379 [inlined]
[7] _loopfooter!(integrator::OrdinaryDiffEqCore.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, false, Vector{Float64}, Nothing, Float64, SciMLBase.NullParameters, Float64, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, Nothing, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, ODEFunction{false, SciMLBase.FullSpecialize, typeof(ode), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, OrdinaryDiffEqCore.InterpolationData{ODEFunction{false, SciMLBase.FullSpecialize, typeof(ode), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, Nothing, OrdinaryDiffEqTsit5.Tsit5ConstantCache, Nothing}, SciMLBase.DEStats, Nothing, Nothing, Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, typeof(ode), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, OrdinaryDiffEqTsit5.Tsit5ConstantCache, OrdinaryDiffEqCore.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{VectorContinuousCallback{typeof(condition), typeof(jump!), typeof(jump!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT), Float64, Int64, Rational{Int64}, Nothing, Nothing, Int64, Tuple{}}}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, Vector{Float64}, Float64, DiffEqBase.CallbackCache{Vector{Float64}, Vector{Float64}}, DiffEqBase.DefaultInit, Nothing})
@ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/GMkz9/src/integrators/integrator_utils.jl:284
[8] loopfooter!
@ ~/.julia/packages/OrdinaryDiffEqCore/GMkz9/src/integrators/integrator_utils.jl:248 [inlined]
[9] solve!(integrator::OrdinaryDiffEqCore.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, false, Vector{Float64}, Nothing, Float64, SciMLBase.NullParameters, Float64, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, Nothing, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, ODEFunction{false, SciMLBase.FullSpecialize, typeof(ode), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, OrdinaryDiffEqCore.InterpolationData{ODEFunction{false, SciMLBase.FullSpecialize, typeof(ode), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, Nothing, OrdinaryDiffEqTsit5.Tsit5ConstantCache, Nothing}, SciMLBase.DEStats, Nothing, Nothing, Nothing, Nothing}, ODEFunction{false, SciMLBase.FullSpecialize, typeof(ode), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, OrdinaryDiffEqTsit5.Tsit5ConstantCache, OrdinaryDiffEqCore.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{VectorContinuousCallback{typeof(condition), typeof(jump!), typeof(jump!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT), Float64, Int64, Rational{Int64}, Nothing, Nothing, Int64, Tuple{}}}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, Vector{Float64}, Float64, DiffEqBase.CallbackCache{Vector{Float64}, Vector{Float64}}, DiffEqBase.DefaultInit, Nothing})
@ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/GMkz9/src/solve.jl:612
[10] #__solve#49
@ ~/.julia/packages/OrdinaryDiffEqCore/GMkz9/src/solve.jl:7 [inlined]
[11] __solve
@ ~/.julia/packages/OrdinaryDiffEqCore/GMkz9/src/solve.jl:1 [inlined]
[12] #solve_call#24
@ ~/.julia/packages/DiffEqBase/6sydR/src/solve.jl:142 [inlined]
[13] solve_call
@ ~/.julia/packages/DiffEqBase/6sydR/src/solve.jl:109 [inlined]
[14] solve_up(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, ODEFunction{false, SciMLBase.FullSpecialize, typeof(ode), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Vector{Float64}, p::SciMLBase.NullParameters, args::Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}; originator::SciMLBase.ChainRulesOriginator, kwargs::@Kwargs{callback::VectorContinuousCallback{typeof(condition), typeof(jump!), typeof(jump!), typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(SciMLBase.FINALIZE_DEFAULT), Float64, Int64, Rational{Int64}, Nothing, Nothing, Int64, Tuple{}}})
@ DiffEqBase ~/.julia/packages/DiffEqBase/6sydR/src/solve.jl:578
[15] solve_up
@ ~/.julia/packages/DiffEqBase/6sydR/src/solve.jl:555 [inlined]
[16] #solve#30
@ ~/.julia/packages/DiffEqBase/6sydR/src/solve.jl:545 [inlined]
[17] top-level scope
@ ~/julia_bug/callback_test.jl:32
[18] include(mod::Module, _path::String)
@ Base ./Base.jl:306
[19] exec_options(opts::Base.JLOptions)
@ Base ./client.jl:317
[20] _start()
@ Base ./client.jl:550
in expression starting at /home/foz/julia_bug/callback_test.jl:32
Describe the bug 馃悶
In integrating a system with a very large number of callbacks (method similar to https://pmc.ncbi.nlm.nih.gov/articles/PMC10470263/ for a large system of biological reactions), a very small proportion of the time solve() will fail with a "Double callback crossing floating pointer reducer errored".
The issue seems to be that in determine_event_occurrence(integrator, callback::VectorContinuousCallback, counter), the value of integrator.last_event_error from the previous event is used to decide whether or not to slightly nudge the point at which the condition is evaluated https://github.com/SciML/DiffEqBase.jl/blob/2ff29cc394503da446f50127c31cc2b889b12b8f/src/callbacks.jl#L215
However, in function find_callback_time(integrator, callback::VectorContinuousCallback, counter), if multiple events are detected in the same timestep, then the value of integrator.last_event_error is changed during the loop over event_idx at https://github.com/SciML/DiffEqBase.jl/blob/2ff29cc394503da446f50127c31cc2b889b12b8f/src/callbacks.jl#L513. If this increases integrator.last_event_error, events with higher indices could potentially search for the event time using a nudged start time https://github.com/SciML/DiffEqBase.jl/blob/2ff29cc394503da446f50127c31cc2b889b12b8f/src/callbacks.jl#L494-L507 , when they were detected using a non-nudged start time.
Expected behavior
Either the event is skipped, or the event fires, without terminating the program.
I suspect the fix is to delay updating integrator.last_event_error until after all event indices have been processed.
Minimal Reproducible Example 馃憞
This is difficult because of the size of the system, and because it is highly floating-point error dependent. However, on my system the following, very artificial example reproduces the bug.
Error & Stacktrace鈿狅笍
Environment (please complete the following information):
using Pkg; Pkg.status()using Pkg; Pkg.status(; mode = PKGMODE_MANIFEST)versioninfo()Additional context
The comments (e.g. at https://github.com/SciML/DiffEqBase.jl/blob/2ff29cc394503da446f50127c31cc2b889b12b8f/src/callbacks.jl#L201) seem to refer to an older version of the code. This bug is likely related to #3220