`VectorContinuousCallback` not picking up threshold

[colebrookson]

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 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]

[ChrisRackauckas]

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.

This looks false. The threshold crosses are recognized. The issue is that the ODE doesn't make sense. You can see this by instrumenting 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

and running it gives:

(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]

so the ODE runs, it sets the value to zero, but your ODE gives a negative derivative to those values so they go negative. You probably want to flip a parameter at those time points to make du[4] =0 etc.