`VectorContinuousCallback` not picking up threshold

General Usage
, ,
1

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)

Screenshot from 2022-01-14 02-16-12
Screenshot from 2022-01-14 02-16-12706×468 40.4 KB

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

Answered in https://github.com/SciML/DifferentialEquations.jl/issues/843. 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.

1 Like