How to debug ODE

Luigi_Marongiu · November 22, 2020, 10:39am

Hello,
I have a set of ODEs to model a bacterial infection:

function kinInfect!(du, u, p, t)
    μ , φ, η, β, ω, χ, ε = p
    H = h = 0
    #=
    du[1]   = susceptible
    du[2]   = infected
    du[3]   = phages
    μ       = growth rate (replication coef), a
    φ       = adsorption rate (transmission coef), b
    η       = lysis rate, k
    β       = burst size, L
    ω       = outflow (decay rate phage), m
    H       = host response against bacterium
    h       = host response against phage
    =#
    x = χ(t, ε)
    du[1] = (μ * u[1]) - (φ * u[1] * u[3]) - (H * u[1])
    du[2] = (μ * u[2]) + (φ * u[1] * u[3]) - (η * u[2]) - (H * u[2])
    du[3] = (x * η * β * u[2]) - (φ * u[1] * u[3]) - (ω * u[3]) - (h * u[3])
end
Chi(x, y) = (x < y) ? 0.0 : 1.0
# parms
tmax    = 20.0    # time span 0-tmax
mu      = 0.5     # growth rate (replication coef), a
phi     = 1e-7    # adsorption rate (transmission coef), b
eta     = 5       # lysis rate, k
beta    = 100     # burst size, L
omega   = 5       # outflow (decay rate phage), m
tmax    = 20.0    # time span 0-tmax
t_in    = 2.5     # time of inoculum
v_in    = 1e8     # amount of inoculum
s0      = 1000.0  # initial susceptible population
i0      = 0.0     # initial infected population
v0      = 0.0     # initial phage population
tspan   = (0.0, tmax)
parms   = [mu, phi, eta, beta, omega, Chi, t_in]
u0      = [s0, i0, v0]
# delay infection
condition(u, t, integrator) = t==t_in           # time of inoculum
affect!(integrator) = integrator.u[3] += 1e8   # amount of inoculum
cb = DiscreteCallback(condition,affect!)
# implement model
prob = ODEProblem(kinInfect!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=cb, tstops=[1000])

the function Chi is there to neutralize the first term of du[3] until the inoculum time is reached. The first time I ran it I got:

But the second time I got:

I don’t remember having changed anything between the first and second time (apart from removing the vertical line).
How can I debug an ODE function? Are there some techniques?
Thank you

MatFi · November 23, 2020, 1:16pm

I guess that it’s the plot(...) command messing it up (probably an erroneous assignment t_in=... )

tomaklutfu · November 23, 2020, 5:10pm

I don’t know the reason but one thing to care about is that you have to stop at 2.5sec(or any time with this kind of discontinuity) to correctly get the changes in Chi(x, y) = (x < y) ? 0.0 : 1.0. Otherwise, the solver allowed to pass that time by large margin if accepted by its error metrics.

Luigi_Marongiu · November 23, 2020, 9:19pm

that may be truly the thing. How would I do that? Tx

thisrod · November 23, 2020, 9:59pm

The t==t_in looks dodgy. Unless you get lucky with rounding errors, that will never be true at a discrete time step, and the inoculation won’t happen. The plotted phage number never gets to 1e8, so it looks like the inoculation isn’t happening.

If I were doing this, I would solve for 2.5s, add 1e8 to the returned solution, then solve again for the remaining time, with that as the initial condition, and with the parameter chi set to 1 instead of 0.

tomaklutfu · November 24, 2020, 3:54am

Simply adding more to tstops=[2.5,1000 #=And more of the discontinuity time=#] or by continuous callbacks.

thisrod · November 24, 2020, 6:33am

Here’s the simplest way to do what I think you want to do.

tmax    = 20.0    # time span 0-tmax
t_in    = 2.5     # time of inoculum
v_in    = 1e8     # amount of inoculum

μ      = 0.5     # growth rate (replication coef), a
φ     = 1e-7    # adsorption rate (transmission coef), b
η     = 5       # lysis rate, k
β    = 100     # burst size, L
ω   = 5       # outflow (decay rate phage), m
H = 0      # host response against bacterium
h = 0    # host response against phage

s0      = 1000.0  # initial susceptible population
i0      = 0.0     # initial infected population
v0      = 0.0     # initial phage population

function kinInfect!(du, u, _, t)
    du[1] = (μ * u[1]) - (φ * u[1] * u[3]) - (H * u[1])
    du[2] = (μ * u[2]) + (φ * u[1] * u[3]) - (η * u[2]) - (H * u[2])
    du[3] = ((t > t_in) * η * β * u[2]) - (φ * u[1] * u[3]) - (ω * u[3]) - (h * u[3])
    @assert du[3] ≈ (η * β * u[2]) - (φ * u[1] * u[3]) - (ω * u[3]) - (h * u[3])
end

prob = ODEProblem(kinInfect!, [s0, i0, v0], (0, t_in))
soln1 = solve(prob, AutoVern7(Rodas5()))
s1, i1, v1 = soln1[end]
v1 += v_in
prob = ODEProblem(kinInfect!, [s1, i1, v1], (t_in, tmax))
soln2 = solve(prob, AutoVern7(Rodas5()))

tt = [soln1.t; soln2.t]
uu = zeros(length(tt), 3)
for (j,u) in pairs([soln1.u; soln2.u])
    uu[j,:] = u
end
plot(tt, uu)

Note:

Parameters defined at the top level are still available in the body of kinInfect!
The boolean expression t>t_in can be promoted to a number, 0 if false, 1 if true. This makes a lot of simulation code simpler. Edit: But it can be left out of this code. If the initial state is in the subspace with no phages and no infections, u[2] == u[3] == 0, then it stays there.
You can use the solution of one problem as the initial value for another.
The plotting code should be plot(soln1); plot!(soln2), but plot! has a bug and replaces the existing plot instead of adding to it.

Luigi_Marongiu · December 8, 2020, 2:43pm

I tried but did not work. I re-arranged as:

function kinInfect!(du, u, p, t)
    μ, φ, η, β, ω, ε = p
    H = h = 0
    #=
    du[1]   = susceptible
    du[2]   = infected
    du[3]   = phages
    μ       = growth rate (replication coef), a
    φ       = adsorption rate (transmission coef), b
    η       = lysis rate, k
    β       = burst size, L
    ω       = outflow (decay rate phage), m
    H       = host response against bacterium
    h       = host response against phage
    =#
    X(x, y) = (x < y) ? 0.0 : 1.0
    du[1] = (μ * u[1]) - (φ * u[1] * u[3]) - (H * u[1])
    du[2] = (μ * u[2]) + (φ * u[1] * u[3]) - (η * u[2]) - (H * u[2])
    du[3] = (X(t,ε) * η * β * u[2]) - (φ * u[1] * u[3]) - (ω * u[3]) - (h * u[3])
    println("t= ", t, "> ", X(t,ε), ": ", du[3])
end 
# parameters
mu  = 0.5       # growth rate, a
phi = 1e-7      # adsorption rate, b
eta = 5         # lysis rate, k
beta = 100      # burst size, L
omega = 5       #  outflow, m
s0 = 1000       # initial susceptible population, x0
i0 = 0.0        # initial infected population, y0
v0 = 0.0        # initial phage population
t_in = 2.5      # time of inoculum, tφ
v_in = 100      # amount of inoculum, vφ
tmax = 20.0     # duration
# instantiate
tspan   = (0.0, tmax)
parms   = [mu, phi, eta, beta, omega, t_in]
u0      = [s0, i0, v0]
condition(u, t, integrator) = t==t_in           # time of inoculum
affect!(integrator) = integrator.u[3] += v_in   # amount of inoculum
cb = DiscreteCallback(condition,affect!)
prob = ODEProblem(kinInfect!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=cb, tstops=[1000])

and I got:

julia>
t= 0.02> 0.0: 0.0
t= 0.0> 0.0: 0.0
t= 0.0011528391772627539> 0.0: 0.0
t= 0.02510627541594441> 0.0: 0.0
t= 0.03765941312391662> 0.0: 0.0
t= 0.10502364904863687> 0.0: 0.0
t= 0.1405332743433084> 0.0: 0.0
t= 0.20382196654005486> 0.0: 0.0
t= 0.21327524779360946> 0.0: 0.0
t= 0.23056783545255075> 0.0: 0.0
t= 0.23056783545255075> 0.0: 0.0
t= 0.23056783545255075> 0.0: 0.0
t= 0.23500388065939562> 0.0: 0.0
t= 0.32717504217939464> 0.0: 0.0
t= 0.37547864554281657> 0.0: 0.0
t= 0.6346915537961186> 0.0: 0.0
t= 0.7713301294032545> 0.0: 0.0
t= 1.0148606280227241> 0.0: 0.0
t= 1.051236198718852> 0.0: 0.0
t= 1.117776876821525> 0.0: 0.0
t= 1.117776876821525> 0.0: 0.0
t= 1.117776876821525> 0.0: 0.0
t= 1.124538601794153> 0.0: 0.0
t= 1.2650322206698672> 0.0: 0.0
t= 1.3386598925940385> 0.0: 0.0
t= 1.7337700218279322> 0.0: 0.0
t= 1.942043929289803> 0.0: 0.0
t= 2.3132498519821483> 0.0: 0.0
t= 2.368695996757698> 0.0: 0.0
t= 2.470121871347117> 0.0: 0.0
t= 2.470121871347117> 0.0: 0.0
t= 2.470121871347117> 0.0: 0.0
t= 2.4779124745446595> 0.0: 0.0
t= 2.639783896538048> 1.0: 0.0
t= 2.7246149091335132> 1.0: 0.0
t= 3.1798458226432613> 1.0: 0.0
t= 3.419811123806132> 1.0: 0.0
t= 3.847500516672675> 1.0: 0.0
t= 3.9113834628925255> 1.0: 0.0
t= 4.028242510855667> 1.0: 0.0
t= 4.028242510855667> 1.0: 0.0
t= 4.028242510855667> 1.0: 0.0
t= 4.037152472484574> 1.0: 0.0
t= 4.2222816752185315> 1.0: 0.0
t= 4.319301257399964> 1.0: 0.0
t= 4.8399400152491> 1.0: 0.0
t= 5.114383671461072> 1.0: 0.0
t= 5.603523726846436> 1.0: 0.0
t= 5.676585412203474> 1.0: 0.0
t= 5.81023483663708> 1.0: 0.0
t= 5.81023483663708> 1.0: 0.0
t= 5.81023483663708> 1.0: 0.0
t= 5.820092888234482> 1.0: 0.0
t= 6.024921293647171> 1.0: 0.0
t= 6.132264522152217> 1.0: 0.0
t= 6.708303337160416> 1.0: 0.0
t= 7.011949956102049> 1.0: 0.0
t= 7.55313835905778> 1.0: 0.0
t= 7.633974382156477> 1.0: 0.0
t= 7.78184515611751> 1.0: 0.0
t= 7.78184515611751> 1.0: 0.0
t= 7.78184515611751> 1.0: 0.0
t= 7.792595310607519> 1.0: 0.0
t= 8.015959631677706> 1.0: 0.0
t= 8.133016869457805> 1.0: 0.0
t= 8.761184230157331> 1.0: 0.0
t= 9.092309304086916> 1.0: 0.0
t= 9.682472469951104> 1.0: 0.0
t= 9.77062373676918> 1.0: 0.0
t= 9.931876054119314> 1.0: 0.0
t= 9.931876054119314> 1.0: 0.0
t= 9.931876054119314> 1.0: 0.0
t= 9.94339866563282> 1.0: 0.0
t= 10.182812927080112> 1.0: 0.0
t= 10.30828136356051> 1.0: 0.0
t= 10.981585962999713> 1.0: 0.0
t= 11.336504173041954> 1.0: 0.0
t= 11.969073769707178> 1.0: 0.0
t= 12.063559184117928> 1.0: 0.0
t= 12.236398356820517> 1.0: 0.0
t= 12.236398356820517> 1.0: 0.0
t= 12.236398356820517> 1.0: 0.0
t= 12.248606138382897> 1.0: 0.0
t= 12.502256710845652> 1.0: 0.0
t= 12.63518588785822> 1.0: 0.0
t= 13.348527257153213> 1.0: 0.0
t= 13.724549999534794> 1.0: 0.0
t= 14.394734137049063> 1.0: 0.0
t= 14.494837945860569> 1.0: 0.0
t= 14.677954669296248> 1.0: 0.0
t= 14.677954669296248> 1.0: 0.0
t= 14.677954669296248> 1.0: 0.0
t= 14.690752829766621> 1.0: 0.0
t= 14.95667016398437> 1.0: 0.0
t= 15.096027911328429> 1.0: 0.0
t= 15.84386708814722> 1.0: 0.0
t= 16.23807461659283> 1.0: 0.0
t= 16.940669440458176> 1.0: 0.0
t= 17.045614356315234> 1.0: 0.0
t= 17.237586763370828> 1.0: 0.0
t= 17.237586763370828> 1.0: 0.0
t= 17.237586763370828> 1.0: 0.0
t= 17.250894453033002> 1.0: 0.0
t= 17.52739867156927> 1.0: 0.0
t= 17.672304625668488> 1.0: 0.0
t= 18.449917291594815> 1.0: 0.0
t= 18.859819282055923> 1.0: 0.0
t= 19.590386295643068> 1.0: 0.0
t= 19.69950935087289> 1.0: 0.0
t= 19.89912469580549> 1.0: 0.0
t= 19.89912469580549> 1.0: 0.0
t= 19.89912469580549> 1.0: 0.0
t= 19.899629072326462> 1.0: 0.0
t= 19.91010889559556> 1.0: 0.0
t= 19.915600995490593> 1.0: 0.0
t= 19.94507339686609> 1.0: 0.0
t= 19.960609146882575> 1.0: 0.0
t= 19.988298464713438> 1.0: 0.0
t= 19.992434352185413> 1.0: 0.0
t= 20.0> 1.0: 0.0
t= 20.0> 1.0: 0.0
julia>

Looks like the test worked but the condition did not…

thisrod · December 9, 2020, 1:05am

However you rearrange the code, the way you’re using condition and affect! simply won’t work. In my last post, I showed you a different way that does work.

Luigi_Marongiu · December 9, 2020, 12:24pm

Sorry, but I do not understand the code. For starter, why affect and condition don’t work, since they work in other instances. Then, if I understand the code, you are splitting the analysis in two: one before t=2.5 and one after. Why not a single one with a split in the middle (the role of having condition)? And what is the role of @assert? Thank you

Luigi_Marongiu · December 9, 2020, 12:46pm

adding more to tstops=[2.5,1000 #=And more of the discontinuity time=#] or by continuous callbacks.

How would I do that?

Luigi_Marongiu · December 9, 2020, 1:05pm

I wrote this:

using DifferentialEquations
using PyPlot

# function
function baseInfect!(du, u, p, t)
    μ, κ, φ, ω, η, β = p
    #=
    du[1] = susceptible
    du[2] = infected
    du[3] = phages
    =#
    du[1] = ((μ * u[1]) * (1 - ((u[1]+u[2])/κ))) - (φ * u[1] * u[3]) - (ω * u[1])
    du[2] = (φ * u[1] * u[3]) - (η * u[2]) - (ω * u[2])
    du[3] = (β * η * u[2]) - (φ * u[1] * u[3]) - (ω * u[3])
end

# parameters
mu    = 0.47        # maximum growth rate susceptible (B. longum)
kappa = 2.2*10^7    # maximum population density
phi   = 10.0^-9     # adsorption rate
omega = 0.05        # outflow
eta   = 1.0         # lyse rate
beta  = 50.0        # burst size
tmax  = 4000.0      # time span 0-tmax
t_in  = 1000.0      # injection time
s0    = 50000.0     # initial susceptible population
i0    = 0.0         # initial infected population
v0    = 0.0         # initial phage population
v_in  = 80.0        # amount of injection

#=
SCENARIO 1
    one host one predator
=#
# instantiate solver
u0 = [s0, i0, v0]
parms = [mu, kappa, phi, omega, eta, beta]
tspan = (0.0, tmax)
# delay infection
condition(u, t, integrator) = t==t_in           # time of inoculum
affect!(integrator) = integrator.u[3] += v_in    # amount of inoculum
cb = DiscreteCallback(condition,affect!)
# instantiate model
prob = ODEProblem(baseInfect!, u0, tspan, parms)
soln = solve(prob, AutoVern7(Rodas5()), callback=cb, tstops=[1000])

In this case, condition and affect do work. Is the problem of time scale? In this example, the injection is at t=1000, in the original problem it was 2.5. In that case, would [tomaklutfu]'s suggestion be the easiest?

Luigi_Marongiu · December 9, 2020, 1:48pm

Interestingly, I tried with t=900 or 1100, and I get this:

So I got lucky the first time choosing 1000 as t_in…

thisrod · December 10, 2020, 12:12am

That way could work too, when you set tstops=[t_in].

Luigi_Marongiu · December 10, 2020, 6:59am

Looks like it is fine this way. Thanks