# 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.

Zooming in on the large negative part near` t=3250`

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.

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

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