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)
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]
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!)
A note unrelated to your question: it might be nicer to use a PresetTimeCallback for this.
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).
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]
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 
Oops, you need the second index at the end then
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
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.
I have tried this with true,true, false,true and vice versa. It still eliminates two species a time…
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:
Yep, that does it. Thanks a lot!
1 Like
