Discrete callback is triggered twice each time

General Usage
,
1

I am working on a Lotka-Volterra type competition system. I am trying to specify a rule such that at pre-specified time-points the (one) species with the lowest abundance at that time-step is eliminated. I am using a discrete callback to achieve that. I am getting stuck at this: the callback works but it eliminates two species at a time. I am still struggling with Julia Types and I suspect that the problem relates to that but I have not been able to isolate it.

I would appreciate any pointers. Here is a simple example that reproduces the problem.

using DifferentialEquations
using Plots

function lotka_volterra!(du, u, p, t)
    r = [0.5, 0.6, 0.7, 0.8]
    alpha = [
        1.0 0.1 0.2 0.3;  
        0.1 1.0 0.4 0.5;  
        0.2 0.4 1.0 0.6;  
        0.3 0.5 0.6 1.0   
    ]
    for i in 1:4
        du[i] = u[i] * r[i] * (1 - sum(alpha[i, j] * u[j] for j in 1:4))
    end
end

times = [100, 200]

condition(u, t, integrator) = t ∈ times

function affect!(integrator)
    # indentify nonzero elements of u
    nonzero = findall(!iszero, integrator.u)
    # find the one with the smallest value
    minu = argmin(integrator.u[nonzero])
    idx = nonzero[minu]
    println("Eliminating species with index $idx at timestep $(integrator.t)")
    integrator.u[idx] = 0.0 
end

cb = DiscreteCallback(condition, affect!)

u0 = [0.5, 0.4, 0.3, 0.2]  # Initial populations of species 1, 2, 3, 4

tspan = (0.0, 250.0)

prob = ODEProblem(lotka_volterra!, u0, tspan, callback = cb)

sol = solve(prob, tstops = times)

plot(sol)

Screenshot from 2025-03-17 15-24-20
Screenshot from 2025-03-17 15-24-20601×396 17.4 KB

Outside of the solver the relevant code works as expected:

u = [2.0, 0.0, 1.0, 3.0]
nonzero = findall(!iszero, u)
minu = argmin(u[nonzero])
idx = nonzero[minu]
u[idx]
2

I still do not understand the behavior. Continuous callback respects the one species rule but eliminates it a couple of steps after the specified timestep (good enough for my use case).

cb = ContinuousCallback(condition, affect!)

Screenshot from 2025-03-17 22-50-56
Screenshot from 2025-03-17 22-50-56603×401 18.5 KB

3

A note unrelated to your question: it might be nicer to use a PresetTimeCallback for this.

4

For your first example: it might be that for some reason your callback is triggered twice each time, which I think would explain this behavior. You could easily check this by adding a print statement to the callback (printing t would be informative).

5

Also, are you sure this code does what you want:

u = [2.0, 0.0, 1.0, 3.0]
nonzero = findall(!iszero, u)
minu = argmin(u[nonzero])
idx = nonzero[minu]
u[idx]

Here you take an index of a subset of u to index into u.

Try

idx = findmin(x -> iszero(x) ? Inf : x, u)[1]
6

This describes the issue better than I did. However I would like to figure out why it happens. I do have a print statement in my example that demonstrates this.

It does. Yours does not; it returns the lowest nonzero value, not its index :wink:

7

Oops, you need the second index at the end then

8

I think sometimes ODE solvers do ‘length 0 timesteps’, I think especially to deal with discontinuities like this. Maybe @ChrisRackauckas knows a way around this, or you store the last t when your callback was called in your parameters and don’t perform the elimination if the time didn’t change

9

In order for that discontinuity to be perfectly vertical, you need to have two points at that time in order to capture the left continuous and right continuous value. This is done by the argument save_positions = (true,true) in the DiscreteCallback constructor, which means save before and save after the callback. If you want to remove one of those saves, simply change it to save_positions = (false,true), for example to only save the post-callback value.

10

I have tried this with true,true, false,true and vice versa. It still eliminates two species a time…

11

It’s not about when the output is saved, it’s about when the callback is triggered. Here’s what I changed:

function affect!(integrator)
    (; u, p, t) = integrator
    if t > p.t_last_elimination
        idx = findmin(x -> iszero(x) ? Inf : x, u)[2]
        println("Eliminating species with index $idx at timestep $t")
        integrator.u[idx] = 0.0
        p.t_last_elimination = t
    end
end

times = [100, 200]
cb = PresetTimeCallback(times, affect!)

mutable struct Parameters
    t_last_elimination::Float64
end

p = Parameters(0.0)

prob = ODEProblem(lotka_volterra!, u0, tspan, p, callback=cb)
sol = solve(prob) # PresetTimeCallback times are automatically added as tstops

This gives what you want, eliminating one species at a time at the right times:

image
image655×436 17.7 KB

12

Yep, that does it. Thanks a lot!

1 Like
13

I am returning to this because I am still struggling with it. Instead of a simplified example, I am posting the code that is pretty close to my actual case. I tried to remove stuff that are not essential but this could have been simplified more.

Using one seed (123), I get the desired behaviour, but with other seeds things fall appart:

using DifferentialEquations
using Plots
using Random
using Distributions
using Statistics
using DataFrames 

# number of time points (365 days, 500 years)
num_time_points = 365 * 50
# generate the time points (in days)
time_points = 0:1:num_time_points-1  

# number of species
N = 16
# store random samples
samples_df = DataFrame(lncR=Float64[], EaR=Float64[], EhR=Float64[], ThR=Float64[],
                       lncP=Float64[], EaP=Float64[], EhP=Float64[], ThP=Float64[]
)

################
Random.seed!(123)
################


# random samples for the parameters
for i in 1:N
    lncR = rand(Normal(-0.09, 0.23))  
    EaR = rand(Normal(1.07, 0.04336814)) 
    EhR = rand(Normal(2.62, 0.1581662))
    ThR = rand(Normal(307.23, 0.7296052))

    lncP = rand(Normal(1.25, 0.22))  
    EaP = rand(Normal(0.74, 0.02806174)) 
    EhP = rand(Normal(6.10, 0.5714391))
    ThP = rand(Normal(306.89, 0.8597097))
    
    
    push!(samples_df, (lncR, EaR, EhR, ThR, lncP, EaP, EhP, ThP))
end


function schoolfield_high(lnc, Ea, Eh, Th, temp, Tc)
    Tc = 273.15 + Tc
    k = 8.62e-05

    boltzmann_term = lnc + log(exp(Ea/k * (1/Tc - 1/temp)))
    inactivation_term = log(1/(1 + exp(Eh/k * (1/Th - 1/temp))))

    return exp(boltzmann_term + inactivation_term)
end


# initial conditions (populations of the species)
u0 = fill(1, N)


# time span for the simulation
tspan = (0, num_time_points-1)
# carrying capacity
K = fill(10.0, N)

# extinction threshold
ext = 1e-2  # any value below this will be set to 0
function condition(out, u, t, integrator) 
    out .= u .- ext # out is the condition or value we want u to change at
end
function affect!(integrator, event_index)
    integrator.u[event_index] = 0.0
end
extinction = VectorContinuousCallback(condition, affect!, length(u0))


# the times where the second callback should trigger
times = 365 .* 2.5 * collect(1:N)

condition2(u, t, integrator) = t ∈ times

# the affect function for the discrete callback
function affect2!(integrator)
    (; u, p, t) = integrator
    if t > p.t_last_elimination

    nonzero = findall(!iszero, integrator.u)
    

    # if there are no nonzero elements, return to stop the callback
        if isempty(nonzero)
            println("No nonzero elements left, stopping callback.")
            return
        end

        minu = argmin(integrator.u[nonzero])
        idx = nonzero[minu]
        # or, randomly select one of the nonzero elements
        #idx = nonzero[rand(1:length(nonzero))]
        println("Eliminating species at index $idx at timestep $(integrator.t) with mass $(integrator.u[idx])")
        integrator.u[idx] = 0.0 # Set that species' population to 0
        p.t_last_elimination = t
    end
end

mutable struct Parameters
    t_last_elimination::Float64
end
p = Parameters(0.0)

# the DiscreteCallback for eliminating the species with the lowest population
extinction2 = DiscreteCallback(condition2, affect2!)


function competition!(du, u, p, t)
    
    T = 12.0

    alpha = zeros(N, N)

    for i in 1:N
            sample = samples_df[i, :]
            lncR_i = sample.lncR
            EaR_i = sample.EaR
            EhR_i = sample.EhR
            ThR_i = sample.ThR
    
            lncP_i = sample.lncP
            EaP_i = sample.EaP
            EhP_i = sample.EhP
            ThP_i = sample.ThP
    
            # temperature-dependent growth rate for species i
            r_T_i = schoolfield_high(lncP_i, EaP_i, EhP_i, ThP_i, T + 273.15, 20)
            # temperature-dependent respiration rate for species i
            x_T_i = schoolfield_high(lncR_i, EaR_i, EhR_i, ThR_i, T + 273.15, 20)

            for j in 1:N
                sample_j = samples_df[j, :]
                
                lncP_j = sample_j.lncP
                EaP_j = sample_j.EaP
                EhP_j = sample_j.EhP
                ThP_j = sample_j.ThP
    
                # temperature-dependent growth rate for species j
                r_T_j = schoolfield_high(lncP_j, EaP_j, EhP_j, ThP_j, T + 273.15, 20)
    
                if j == i
                    alpha[i, j] = 1.0
                else
                # calculate the competition coefficient a_ij
                    alpha[i, j] = .4 * (r_T_i / r_T_j)^.2
                end
            end
    

        du[i] = u[i] * (r_T_i * (1 - sum(alpha[i, j] * u[j] for j in 1:N) / K[i]) - x_T_i)

    end
end

prob = ODEProblem(competition!, u0, tspan, p, callback = CallbackSet(extinction
                                                                 ,extinction2
                                                                  ), saveat = 1)
sol = solve(prob, tstops = times)

Eliminating species at index 8 at timestep 912.5 with mass 0.25347481869374444
Eliminating species at index 6 at timestep 1825.0 with mass 0.6310769829516735
Eliminating species at index 10 at timestep 2737.5 with mass 0.9805875134502837
Eliminating species at index 7 at timestep 3650.0 with mass 1.1496895611450328
Eliminating species at index 3 at timestep 4562.5 with mass 1.2897040212490765
Eliminating species at index 4 at timestep 5475.0 with mass 1.5346786416125333
Eliminating species at index 9 at timestep 6387.5 with mass 1.768617170665732
Eliminating species at index 14 at timestep 7300.0 with mass 2.2663432331842195
Eliminating species at index 1 at timestep 8212.5 with mass 2.5733422546046643
Eliminating species at index 15 at timestep 9125.0 with mass 2.9832460675253616
Eliminating species at index 13 at timestep 10037.5 with mass 3.6375138501894084
Eliminating species at index 5 at timestep 10950.0 with mass 4.803604072567314
Eliminating species at index 2 at timestep 11862.5 with mass 6.198789865552722
Eliminating species at index 11 at timestep 12775.0 with mass 8.72034083767297
No nonzero elements left, stopping callback.
No nonzero elements left, stopping callback.
No nonzero elements left, stopping callback.
No nonzero elements left, stopping callback.

Screenshot from 2025-04-02 20-04-26
Screenshot from 2025-04-02 20-04-261552×349 18.7 KB

Using a different seed (321)

Eliminating species at index 12 at timestep 912.5 with mass -2.40132276424625e-18
Eliminating species at index 1 at timestep 1825.0 with mass -1.2879681645311077e-21
Eliminating species at index 1 at timestep 2737.5 with mass -5.095014003689226e-25
Eliminating species at index 1 at timestep 3650.0 with mass -6.368769688430557e-26
Eliminating species at index 1 at timestep 4562.5 with mass 7.438509326714505e-43
Eliminating species at index 1 at timestep 5475.0 with mass -3.731700989210302e-28
Eliminating species at index 1 at timestep 6387.5 with mass -3.4984684801399734e-29
Eliminating species at index 1 at timestep 7300.0 with mass 1.7492348386921838e-29
Eliminating species at index 1 at timestep 8212.5 with mass 1.9435939326455973e-29
Eliminating species at index 1 at timestep 9125.0 with mass -1.399387870953741e-28
Eliminating species at index 1 at timestep 10037.5 with mass -1.4576956989101483e-29
Eliminating species at index 1 at timestep 10950.0 with mass -5.830780800233315e-30
Eliminating species at index 12 at timestep 11862.5 with mass 2.146825699968211e-30
Eliminating species at index 1 at timestep 12775.0 with mass -5.636422562641964e-29
Eliminating species at index 1 at timestep 13687.5 with mass -1.5937467520637673e-28
Eliminating species at index 12 at timestep 14600.0 with mass 2.047720931595638e-30

231:

Eliminating species at index 7 at timestep 912.5 with mass 0.08135520084356326
Eliminating species at index 7 at timestep 1825.0 with mass -4.499509540064276e-18
Eliminating species at index 7 at timestep 2737.5 with mass -1.0887596576718764e-23
Eliminating species at index 7 at timestep 3650.0 with mass 2.295574815477395e-25
Eliminating species at index 7 at timestep 4562.5 with mass -1.2297747206782116e-25
Eliminating species at index 7 at timestep 5475.0 with mass 4.919089239682802e-26
Eliminating species at index 7 at timestep 6387.5 with mass -5.8926701226329334e-27
Eliminating species at index 7 at timestep 7300.0 with mass -3.202531461426503e-28
Eliminating species at index 7 at timestep 8212.5 with mass -2.401901580792794e-28
Eliminating species at index 7 at timestep 9125.0 with mass -3.7779947850981226e-28
Eliminating species at index 7 at timestep 10037.5 with mass -4.0031674484660875e-29
Eliminating species at index 7 at timestep 10950.0 with mass -1.2109593080413522e-28
Eliminating species at index 7 at timestep 11862.5 with mass -1.726365719052915e-29
Eliminating species at index 7 at timestep 12775.0 with mass 1.1308952380635951e-28
Eliminating species at index 7 at timestep 13687.5 with mass -5.20412415612934e-29
Eliminating species at index 7 at timestep 14600.0 with mass 2.401899212962485e-28

There is nothing unreasonable with the parameters that are generated based on the random values. The callback just fails…

Screenshot from 2025-04-02 20-11-57
Screenshot from 2025-04-02 20-11-571552×349 12.6 KB

14

Can you open an issue to DiffEqBase? Otherwise it won’t get in the email list.

1 Like
15

Done.

16

My bad. All this because I did not explicitly specify a solver. Tsit5 gets the job done. It looks like in some cases some other solver was chosen and that was the source of the problem.