How to change u to a specific value once a value (not tstep) threshold is crossed?

Hello again,

Is there a way to change u (in my case bioS) in a differential equation to a specific value once a value threshold has been crossed?

I am running a simulation of a food web and once a species in that food web reaches a biomass of <0.000001, I want it to be set to 0.0. This is the extinction threshold in my model. The model runs but after inspecting individual species I’ve noticed that sometimes bioS becomes very small (<10^-7) or even negative. bioS is currently an array.

What is the correct way to set the extinction threshold? I specify this rule at the beginning of my function:

# calculating change in biomass of nutrients, plants, animals, and vessel
function dBdT(dB, bioS, p, t)
     D, S, K, v, r, X, num_nutrient, num_plant, num_animal, num_vessel_types, bm, foodWeb, b, q, ω, c, h, e, target_species, catch_max, b0V, μ, elasticity, vessel_attack_type = p

    # extinction threshold
    @. bioS[bioS<0.000001] = 0.0
[...]

Just adding some images to explain what I mean.

image

Zooming in on the large negative part near t=3250

image

Create a DiscreteCallback or a ContinuousCallback that sets the value to zero when the threshold is crossed. DiscreteCallback will do it at the end of a step, while ContinuousCallback has more overhead but will do it at exactly the point of the crossing. Which one is better depends on the application.

I’ve made the following callback but run into an error.

bioS is what I call u, and it is an Array. Julia is not happy with something, and I assume it is me not correctly specifying that if any value of bioS is <= 0.000001, that value of bioS becomes 0.0. What is the correct way to specify this?

function condition(bioS, t, integrator)
    bioS[bioS] .<= 0.000001
end

function affect!(integrator)
    for c in full_cache(integrator)
        bioS[bioS] .= 0.0
    end
end

extinction_cb = ContinuousCallback(condition,affect!)

p = (D, S, K, v, r, X, num_nutrient, num_plant, num_animal, num_vessel_types, bm, foodWeb, b, q, ω, c, h, e, target_species, catch_max, b0V, μ, elasticity, demand_foraging)
prob = ODEProblem(dBdT,bioS,tspan,p, callback=extinction_cb)

base_simulation=solve(prob, alg = CVODE_BDF())
BoundsError: attempt to access 133-element Array{Float64,1} at index [[7.240610761864229, 8.157086597412007, 7.197820763387717, 2.7759222193200817, 1.5273135178916486, 7.06140775166741, 8.552645764565877, 7.7183298686176, 9.201408232291614, 7.245864412676511  …  6.703892924721959, 4.065422576584616, 7.276238849282553, 2.4681265775503336, 8.443596981833313, 7.481300380521403, 3.66466910401394, 5.883428618363283, 8.073258561352732, 0.00029764921209184553]]
throw_boundserror(::Array{Float64,1}, ::Tuple{Array{Float64,1}}) at abstractarray.jl:538
checkbounds at abstractarray.jl:503 [inlined]
_getindex at multidimensional.jl:669 [inlined]
getindex at abstractarray.jl:981 [inlined]
condition(::Array{Float64,1}, ::Float64, ::Sundials.CVODEIntegrator{Array{Float64,1},Tuple{Float64,Array{Float64,1},Array{Float64,2},Array{Float64,1},Array{Float64,1},Array{Float64,1},Int64,Int64,Int64,Int64,Array{Float64,1},Array{Float64,2},Array{Float64,2},Float64,Array{Float64,1},Float64,Array{Float64,2},Array{Float64,1},Int64,Float64,Float64,Float64,Float64,typeof(demand_foraging)},Sundials.Handle{Sundials.CVODEMem},ODESolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Nothing,ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Tuple{Float64,Array{Float64,1},Array{Float64,2},Array{Float64,1},Array{Float64,1},Array{Float64,1},Int64,Int64,Int64,Int64,Array{Float64,1},Array{Float64,2},Array{Float64,2},Float64,Array{Float64,1},Float64,Array{Float64,2},Array{Float64,1},Int64,Float64,Float64,Float64,Float64,typeof(demand_foraging)},ODEFunction{true,typeof(dBdT),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Symbol,ContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing},Tuple{Symbol},NamedTuple{(:callback,),Tuple{ContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}}}},DiffEqBase.StandardODEProblem},CVODE_BDF{:Newton,:Dense,Nothing,Nothing},DiffEqBase.HermiteInterpolation{Array{Float64,1},Array{Array{Float64,1},1},Array{Array{Float64,1},1}},DiffEqBase.DEStats},CVODE_BDF{:Newton,:Dense,Nothing,Nothing},ODEFunction{true,typeof(dBdT),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Sundials.FunJac{ODEFunction{true,typeof(dBdT),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,Nothing,Tuple{Float64,Array{Float64,1},Array{Float64,2},Array{Float64,1},Array{Float64,1},Array{Float64,1},Int64,Int64,Int64,Int64,Array{Float64,1},Array{Float64,2},Array{Float64,2},Float64,Array{Float64,1},Float64,Array{Float64,2},Array{Float64,1},Int64,Float64,Float64,Float64,Float64,typeof(demand_foraging)},Nothing,Nothing,Array{Float64,1},Nothing,Nothing,Nothing},Nothing,Sundials.DEOptions{DataStructures.BinaryHeap{Float64,DataStructures.LessThan},DataStructures.BinaryHeap{Float64,DataStructures.LessThan},CallbackSet{Tuple{ContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}},Tuple{}},Float64,Float64,typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE)},Array{Float64,1},Tuple{Int64},Array{Float64,1},Sundials.LinSolHandle{Sundials.Dense},Sundials.MatrixHandle{Sundials.DenseMatrix},Nothing}) at untitled-2b415c32029fe4bb83d76568e874f43e:6
determine_event_occurance at callbacks.jl:456 [inlined]
find_callback_time(::Sundials.CVODEIntegrator{Array{Float64,1},Tuple{Float64,Array{Float64,1},Array{Float64,2},Array{Float64,1},Array{Float64,1},Array{Float64,1},Int64,Int64,Int64,Int64,Array{Float64,1},Array{Float64,2},Array{Float64,2},Float64,Array{Float64,1},Float64,Array{Float64,2},Array{Float64,1},Int64,Float64,Float64,Float64,Float64,typeof(demand_foraging)},Sundials.Handle{Sundials.CVODEMem},ODESolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Nothing,ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Tuple{Float64,Array{Float64,1},Array{Float64,2},Array{Float64,1},Array{Float64,1},Array{Float64,1},Int64,Int64,Int64,Int64,Array{Float64,1},Array{Float64,2},Array{Float64,2},Float64,Array{Float64,1},Float64,Array{Float64,2},Array{Float64,1},Int64,Float64,Float64,Float64,Float64,typeof(demand_foraging)},ODEFunction{true,typeof(dBdT),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Symbol,ContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing},Tuple{Symbol},NamedTuple{(:callback,),Tuple{ContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}}}},DiffEqBase.StandardODEProblem},CVODE_BDF{:Newton,:Dense,Nothing,No...

I have never used callbacks in Julia, so I am not an expert…

To me, your construct bioS[bioS] looks strange (i.e., in both function condition and function affect!). A question:

  • Is bioS is a vector?

According to the documentation, it seems like the callback is triggered if the return value from the condition function becomes zero in the time interval.

Would the following condition function work?

function condition(bioS, t, integrator)
    minimum(bioS[:] .>= 0.000001)
end

I don’t know if this works… In the condition function above, if a single element in vector bioS drops below 0.000001 (or better: 1e-6; your number 0.000001 is 10^{-6}, not 10^{-7}), then the condition of that element becomes false which is equal to 0. Thus, minimum(bioS[:].>1e-6)becomes false or 0 if a single value of bioS drops below 1e-6.

So – I don’t know if false is accepted to trigger a callback.

https://docs.juliadiffeq.org/latest/features/callback_functions/

Indeed it’s strange.

function condition(bioS, t, integrator)
    any(bioS .<= 0.000001)
end

function affect!(integrator)
  bioS[bioS .<= 0.000001] .= 0.0
end

extinction_cb = DiscreteCallback(condition,affect!)

is the easiest one. Then for exactness:

function condition(bioS, t, integrator)
    bioS .- 0.000001
end

function affect!(integrator,event_index)
  bioS[event_index] = 0.0
end

extinction_cb = VectorContinuousCallback(condition,affect!,length(u0)) # whatever the length is

And I noticed that the docstrings in the docs here were especially bad, so I just fixed that: typedef was an awfully useless description of how to build the type, … by ChrisRackauckas · Pull Request #440 · SciML/DiffEqBase.jl · GitHub

1 Like

Yep. I found the example… I just didn’t scroll down in the documentation…

OK… so the condition function should evaluate to true if discrete callback is used, and to zero if continuous (or vector continuous) callback is used.

Yes, one of them is a true/false at the end of the time step: if true at the end of the time step, then do…

The other is a rootfinding condition: at exactly the time point at which condition == 0, then do … so this will reverse time a bit and make the condition happen precisely at the right time point.

The first is to add control flow, the latter is for things like the bouncing ball, where you need to keep the ball above the ground at all time!

1 Like

“… bouncing ball”, or the water level in a bath tub which at some stage starts to overflow, I guess. Although, continuous detection is more critical for the bouncing ball :-).

Yup, in those cases you don’t want the water to rise a meter above the bathtub because you had a big dt, and then suddenly overflow: you’d want a discontinuous change to your system the moment it starts to overflow.

That doesn’t mean the DiscreteCallback is useless though, it just has some very different properties. It’s great for doing things like, if x > 100, plot the output. Or things like SavingCallback are made by using it. Or you can do things like project to a manifold after each step, which is a nice way to get energy conservation and other properties. So you really just have to choose the right one for the job.

I think in this case, if negatives for a small amount of time are okay in the model, then a DiscreteCallback to zero anything that’s sufficiently small at the end of a timestep is a good solution. Just note that you’ll likely want to use saveat here instead of the interpolation, since you’ve now manually changed values. But if any slight negative could really mess up the results, then zeroing immediately when it’s below a threshold via ContinuousCallback is the right solution. It’s a modeling choice the user has to make.

2 Likes

Thank you for your reply,

I’ve changed my code based on your example. However I now run into an error I can’t find much information on. Any idea what is happening?

tmax = 10000.0
tspan = (0.0,tmax)
const target_species = 0

function condition(bioS, t, integrator)
    bioS .- 0.000001
end

function affect!(integrator,event_index)
    bioS[event_index] = 0.0
end

extinction_cb = VectorContinuousCallback(condition,affect!, length(bioS))

p = (D, S, K, v, r, X, num_nutrient, num_plant, num_animal, num_vessel_types, bm, foodWeb, b, q, ω, c, h, e, target_species, catch_max, b0V, μ, elasticity, demand_foraging)
prob = ODEProblem(dBdT,bioS,tspan,p, callback=extinction_cb)

base_simulation=solve(prob, alg = CVODE_BDF())
UndefVarError: uBottomEltype not defined
#__init#23(::Bool, ::VectorContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}, ::Float64, ::Float64, ::Array{Float64,1}, ::Array{Float64,1}, ::Int64, ::Nothing, ::Float64, ::Float64, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::String, ::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), ::Nothing, ::Bool, ::Bool, ::Nothing, ::Bool, ::Base.Iterators.Pairs{Symbol,CVODE_BDF{:Newton,:Dense,Nothing,Nothing},Tuple{Symbol},NamedTuple{(:alg,),Tuple{CVODE_BDF{:Newton,:Dense,Nothing,Nothing}}}}, ::typeof(DiffEqBase.__init), ::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Tuple{Float64,Array{Float64,1},Array{Float64,2},Array{Float64,1},Array{Float64,1},Array{Float64,1},Int64,Int64,Int64,Int64,Array{Float64,1},Array{Float64,2},Array{Float64,2},Float64,Array{Float64,1},Float64,Array{Float64,2},Array{Float64,1},Int64,Float64,Float64,Float64,Float64,typeof(demand_foraging)},ODEFunction{true,typeof(dBdT),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Symbol,VectorContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing},Tuple{Symbol},NamedTuple{(:callback,),Tuple{VectorContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}}}},DiffEqBase.StandardODEProblem}, ::CVODE_BDF{:Newton,:Dense,Nothing,Nothing}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}) at solve.jl:69
(::DiffEqBase.var"#kw##__init")(::NamedTuple{(:callback, :alg),Tuple{VectorContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing},CVODE_BDF{:Newton,:Dense,Nothing,Nothing}}}, ::typeof(DiffEqBase.__init), ::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Tuple{Float64,Array{Float64,1},Array{Float64,2},Array{Float64,1},Array{Float64,1},Array{Float64,1},Int64,Int64,Int64,Int64,Array{Float64,1},Array{Float64,2},Array{Float64,2},Float64,Array{Float64,1},Float64,Array{Float64,2},Array{Float64,1},Int64,Float64,Float64,Float64,Float64,typeof(demand_foraging)},ODEFunction{true,typeof(dBdT),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Symbol,VectorContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing},Tuple{Symbol},NamedTuple{(:callback,),Tuple{VectorContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}}}},DiffEqBase.StandardODEProblem}, ::CVODE_BDF{:Newton,:Dense,Nothing,Nothing}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}) at none:0
#init_call#440(::Base.Iterators.Pairs{Symbol,Any,Tuple{Symbol,Symbol},NamedTuple{(:callback, :alg),Tuple{VectorContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing},CVODE_BDF{:Newton,:Dense,Nothing,Nothing}}}}, ::typeof(DiffEqBase.init_call), ::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Tuple{Float64,Array{Float64,1},Array{Float64,2},Array{Float64,1},Array{Float64,1},Array{Float64,1},Int64,Int64,Int64,Int64,Array{Float64,1},Array{Float64,2},Array{Float64,2},Float64,Array{Float64,1},Float64,Array{Float64,2},Array{Float64,1},Int64,Float64,Float64,Float64,Float64,typeof(demand_foraging)},ODEFunction{true,typeof(dBdT),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Symbol,VectorContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing},Tuple{Symbol},NamedTuple{(:callback,),Tuple{VectorContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing}}}},DiffEqBase.StandardODEProblem}, ::CVODE_BDF{:Newton,:Dense,Nothing,Nothing}, ::Vararg{Any,N} where N) at solve.jl:11
#init_call at none:0 [inlined]
#init#441 at solve.jl:30 [inlined]
#init at none:0 [inlined]
#__solve#22(::Base.Iterators.Pairs{Symbol,Any,Tuple{Symbol,Symbol},NamedTuple{(:callback, :alg),Tuple{VectorContinuousCallback{typeof(condition),typeof(affect!),typeof(affect!),typeof(DiffEqBase.INITIALIZE_DEFAULT),Float64,Int64,Nothing},CVODE_BDF{:Newton,:Dense,Nothing,Nothing}}}}, ::typeof(DiffEqBase.__solve), ::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Tuple{Float64,Array{Float64,1},Array{Float64,2},Array{Float64,1},Array{Float64,1},Array{Float64,1},Int64,Int64,Int64,Int64,Array{Float64,1},Array{Float64,2},Array{Float64,2},Float64,Array{Float64,1},Float64,Array{Float64,2},Array{Float64,1},Int64,Float64,Float64,Float64,Float64,typeof(demand_foragin...

I am choosing to do a continuous callback as I get a DomainError... Exponentiation yielding a complex result requires a complex argument. Replace x^y with (x+0im)^y, Complex(x)^y, or similar, otherwise. I know this error is caused by trying to calculate a negative number with an exponent (no negatives should be allowed).

This is being fixed by https://github.com/JuliaDiffEq/Sundials.jl/pull/233

Hi, the UndefVarError: uBottomEltype not defined has been fixed, but I still get DomainError with -4.730795533123454e-6: Exponentiation yielding a complex result requires a complex argument. Replace x^y with (x+0im)^y, Complex(x)^y, or similar.. This means some values are still turning negative.

The error is stemming from this part of my code

#Feeding Function
function FF(num_nutrient, bm, foodWeb, pred, prey, bioS, b, h, ω, q, c)
    sumF = 0.0

    for k in num_nutrient+1:length(bm)
        if foodWeb[k, pred] != 0.0
            sumF = sumF+b[k, pred]*bioS[k]^(1+q) # error highlights this section
        end
    end
    return (((ω[pred]*b[prey, pred]*bioS[prey]^(1+q))/(1+c*bioS[pred]+ω[pred]*h[pred]*sumF))/bm[pred])
end

What am I doing wrong with the callback?

tmax = 10000.0
tspan = (0.0,tmax)
const target_species = 0

function condition(bioS, t, integrator)
    bioS .- 0.000001
end

function affect!(integrator,event_index)
    bioS[event_index] = 0.0
end

extinction_cb = VectorContinuousCallback(condition,affect!, length(bioS))

p = (D, S, K, v, r, X, num_nutrient, num_plant, num_animal, num_vessel_types, bm, foodWeb, b, q, ω, c, h, e, target_species, catch_max, b0V, μ, elasticity, demand_foraging)
prob = ODEProblem(dBdT,bioS,tspan,p, callback=extinction_cb)

base_simulation=solve(prob, alg = CVODE_BDF())

Make that sumF = sumF+b[k, pred]*abs(bioS[k])^(1+q)

Thanks to @ChrisRackauckas for finding the error, here is how it works.

tspan = (0.0,1000.0) # how long your simulation runs

function condition(out,u, t, integrator) 
    out .= u .- 1e-6 # out is the condition or value you want u to change at
end

function affect!(integrator,event_index)
    bioS[event_index] = 0.0 # event_index is whatever you want u to change to
end

cb = VectorContinuousCallback(condition,affect!, length(u))

p = (#,#,#) # whatever your parameters are
prob = ODEProblem(your_function,u,tspan,p, callback=cb)

solution =solve(prob)
1 Like