# `VectorContinuousCallback` not picking up threshold

I’m having trouble getting `VectorContinuousCallback` to work as desired and I’m not sure what I’m doing wrong. I have a large system of equations, and essentially, any time any of the values cross some threshold value (in my system it’s 10e-30 but in this reprex 0.05), I want the value to go to zero.

That is, if at any point the values of `u` go below 0.05, I want the callback to take the value to zero, but right now, the solver seems to just almost ignore the callback? Not any of the crosses of the threshold are recognized.

A reprex:

``````using DifferentialEquations, Plots

function biomass_sim!(du, u, p, t)
# change = growth + gain from eating - loss from eating - loss
du[1] = 0.2*u[1] + (0.1*u[2] + 0.15*u[3]) - (0.2*u[4]) - 0.9*u[1]
du[2] = 0.2*u[2] + (0.1*u[1] + 0.05*u[3]) - (0.1*u[1] + 0.4*u[4]) - 0.5*u[2]
du[3] = 1.2*u[3] + 0 - (0.15*u[1] + 0.005*u[2]) - 1.3*u[3]
du[4] = 0.2*u[4] + (0.2*u[1] + 0.4*u[2]) - 1.9*u[1]
end

# set up extinction callback
function extinction_threshold(out,u,t,integrator)
# loop through all species to make the condition check all of them
for i in 1:4
out[i] = 0.05 - u[i]
end
end

function extinction_affect!(integrator, event_idx)
# loop again through all species
for i in 1:4
if event_idx == i
integrator.u[i] = 0
end
end
end
extinction_callback =
VectorContinuousCallback(extinction_threshold,
extinction_affect!,
4,
save_positions = (true, true),
interp_points = 1000
)

tspan = (0.0, 10.0)
u0 = [10, 10, 10, 10]

prob = ODEProblem(biomass_sim!,
u0,
tspan)
sol= solve(prob,
Tsit5(),
abstol = 1e-15,
reltol = 1e-10,
callback = extinction_callback,
progress = true,
progress_steps = 1)
plot(sol)
``````

What I want to see here is that the two values that DO cross the threshold to be < 0.05 at least visually (u3 and u4 clearly go below zero), I want those values to become zero.

The application of this is an extinction to a species, so if they go below some threshold, I want to consider them extinct and therefore not able to be consumed by other species.

I’ve tried changing the tolerances, using different solvers (I’m not married to `Tsit5()`), but have yet to find a way to do this.

Any help much appreciated!!

The full output:

``````retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 165-element Vector{Float64}:
0.0
0.004242012928189618
0.01597661154828718
0.03189297583643294
0.050808376324350105
⋮
9.758563212772982
9.850431863368996
9.94240515017787
10.0
u: 165-element Vector{Vector{Float64}}:
[10.0, 10.0, 10.0, 10.0]
[9.972478301795496, 9.97248284719235, 9.98919421326857, 9.953394015118882]
[9.89687871881005, 9.896943262844019, 9.95942008072302, 9.825050883392004]
[9.795579619066798, 9.795837189358107, 9.919309541156514, 9.652323759097303]
[9.677029343447844, 9.67768414040866, 9.872046236050455, 9.449045554530718]
⋮
[43.86029800419986, 110.54225328286441, -12.173991695732434, -186.40702483057268]
[45.33660997599057, 114.35164541304869, -12.725800474246844, -192.72257104623995]
[46.86398454218351, 118.2922142830212, -13.295579652115606, -199.25572621838901]
[47.84633546050675, 120.82634905853745, -13.661479860003494, -203.45720035095707]
``````

Answered in `VectorContinuousCallback` not picking up threshold · Issue #843 · SciML/DifferentialEquations.jl · GitHub. This was a “user error”. When you check the callback:

``````function extinction_affect!(integrator, event_idx)
# loop again through all species
@show integrator.u,event_idx
for i in 1:4
if event_idx == i
integrator.u[i] = 0
end
end
biomass_sim!(get_tmp_cache(integrator)[1], integrator.u, integrator.p, integrator.t)
@show get_tmp_cache(integrator)[1]
end

``````

It is definitely called and does exactly as intended.

``````(integrator.u, event_idx) = ([5.347462662161639, 5.731062469090074, 7.64667777801325, 0.05000000000000008], 4)
(get_tmp_cache(integrator))[1] = [-2.0231159499021523, -1.3369848518263594, -1.5954424894710222, -6.798261538038756]
(integrator.u, event_idx) = ([12.968499097445866, 30.506371944743357, 0.050000000000001314, -53.521085634736835], 3)
(get_tmp_cache(integrator))[1] = [4.676904953209599, 12.256522670471728, -2.0978067243405967, -20.548116814707996]
``````

but what this also shows is that, even if `u[4] = 0`, `du[4] < 0` and so it’s clear why it goes negative: that’s due to how the ODE is defined. You should flip a parameter or something to make the derivative = 0 if you want to keep it at zero past the callback point.