# Zombies in biological ODE : why is my solver not sticking to zero?

**URL:** <https://discourse.julialang.org/t/zombies-in-biological-ode-why-is-my-solver-not-sticking-to-zero/90409>\
**Category:** Numerics\
**Tags:** ode, precision, differentialequation, diffeqcallbacks\
**Created:** [November 17, 2022, 5:34pm UTC](https://discourse.julialang.org/t/zombies-in-biological-ode-why-is-my-solver-not-sticking-to-zero/90409 "2022-11-17T17:34:01Z")\
**Posts on this page:** 6\
**Page:** 2

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [November 19, 2022, 3:33pm UTC](https://discourse.julialang.org/t/zombies-in-biological-ode-why-is-my-solver-not-sticking-to-zero/90409/22 "2022-11-19T15:33:44Z")

</div>

> [@iago-lito](#):
>
> FWIU now, using a callback to force this “absorbing” behaviour of the value `0` is actually a creepy way of _changing_ what `dBdt!` actually means, and making the ODE system actually simulated _different_ from the system described by `dBdt!` alone.

It’s essentially a discrete way of implementing the “attracting” modification I proposed before. But I think it makes your system more realistic. When the number of rhinos goes to 1, then the rhino becomes extinct due to inability to reproduce…

In any case, I found in my zombies model that explicit timestepping Tsit5 worked well and a callback with Z\<1e-6 zeroing out the zombie population made sense (imagine zombies measured as fraction of initial population of 1M… Z below this threshold means less than 1 zombie.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [November 19, 2022, 8:53pm UTC](https://discourse.julialang.org/t/zombies-in-biological-ode-why-is-my-solver-not-sticking-to-zero/90409/23 "2022-11-19T20:53:42Z")

</div>

No need for the zombie voodoo. I’ve already written down what’s happening here for your colleague. See the issue:

> <https://github.com/SciML/OrdinaryDiffEq.jl/issues/1651>
>
> This was noted in the MWE https://github.com/HelgavonLichtenstein/Food-Webs by @…HelgavonLichtenstein. 
> 
> I created a debugging script https://github.com/HelgavonLichtenstein/Food-Webs/pull/2 to see what's going on. The end result is captured in:
> 
> \`\`\`julia
> (A\\b)\[101\] # -2.6218322568017613e-23
> (sparse(A)\\b)\[101\] # -0.0 
> \`\`\`
> 
> Basically, the linear solver introduces a numerical error which then revives a dead quantity. It's dead because:
> 
> \`\`\`julia
> function condition(out, bioS, t, integrator)
> # extinction threshold for species
> out\[1:(num\_plant+num\_animal)\] .= bioS\[(num\_nutrient+1):(index\_vessel-1)\] .- 1e-6
> end
> 
> function affect!(integrator, event\_idx)
> integrator.p\[1\]\[event\_idx+num\_nutrient\] = 0.0
> integrator.u\[event\_idx+num\_nutrient\] = 0.0
> @show event\_idx\[\], integrator.t
> end
> \`\`\`
> 
> enforces \`u == 0\` and \`du == 0\` for that component, yet it can still revive because:
> 
> \`\`\`julia
> if repeat\_step
> linres = dolinsolve(integrator, cache.linsolve; A = nothing, b = \_vec(linsolve\_tmp),
> du = cache.fsalfirst, u = u, p = p, t = t, weight = weight, solverdata = (;gamma = dtgamma))
> else
> linres = dolinsolve(integrator, cache.linsolve; A = W, b = \_vec(linsolve\_tmp),
> du = cache.fsalfirst, u = u, p = p, t = t, weight = weight, solverdata = (;gamma = dtgamma))
> end
> 
> @inbounds @simd ivdep for i in eachindex(u)
> k1\[i\] = -linres.u\[i\]
> end
> integrator.destats.nsolve += 1
> 
> @inbounds @simd ivdep for i in eachindex(u)
> u\[i\] = uprev\[i\] + a21\*k1\[i\]
> end
> @show u\[101\], uprev\[101\], k1\[101\]
> \`\`\`
> 
> in the numerical error of \`k1\` introduced by the factorization. Sparse matrices can impose this as a true zero, though that can lose performance for smaller sizes. The right answer is probably just doing \`u\[101\] = 0\` in the ODE, at least for performance.

Basically, the issue is that LU-factorization based solving, i.e. `A\b`, is not guaranteed to preserve zeros. This means that there can be a true zero, but `A\b` can give `2e-16` or something close to `eps(Float64)` like that. When that happens, you get a `ki = 2e-16` where it should be zero, and that gives a positive (or negative) derivative to a “dead” term.

There are many solutions to it. For one, if you use an explicit integrator then LU-factorization is never used so you’re good. Or if you change to a Krylov method like @isaacsas mentioned then you’re fine. Or you can use the trick from the issue to use mutation in the `f` function:

```julia
    for i in (num_nutrient+1):(index_vessel-1)
        if iszero(p[1][i])
            bioS[i] = 0.0 #safety against https://github.com/SciML/OrdinaryDiffEq.jl/issues/1651
        end
    end

```

In almost all cases, mutating `u` is a bad idea. In this specific case, it’s actually a good way to improve numerical accuracy.

---

<div class="post-metadata">

**Author:** ![iago-lito](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/iago-lito/32/31924_2.png) [@iago-lito](https://discourse.julialang.org/u/iago-lito)\
**Post date:** [November 20, 2022, 9:39am UTC](https://discourse.julialang.org/t/zombies-in-biological-ode-why-is-my-solver-not-sticking-to-zero/90409/24 "2022-11-20T09:39:29Z")

</div>

Hi again, and many thanks to @kcqpo and @ChrisRackauckas for nailing this down \<3

So IIUC now, here are the eventual answers to our questions:

- `dBdt!` gets called with nonzero values of `B[2]` because some internal operation used by `solve()` is **not zero-stable** (for performance reasons I guess), namely something called “LU-factorization / `A\b`” that we need to get familiar with. This being said:

- Mutating `DiffeqCallbacks` are not the worst tools for the job, but definitely not the only ones. Alternate approaches to avoiding zombies in our situation include:

- Stop the solver as soon as `B[i] < 0` is detected and spawn a new simulation with the variable `B[i]` removed. This would trigger a freeze on every extinction event. But in large systems with `len(B) > 100` and massive extinctions, performances would actually _increase_ as extinctions go because the system would shrink with time and `dBdt!` would not remain busy herding zeroes.

Not sure exactly which approach to pick now, but we do have very interesting new chunks of information to crunch. In any case KUDOS to y’all folks, you are what makes this community awesome! \<3

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [November 20, 2022, 2:31pm UTC](https://discourse.julialang.org/t/zombies-in-biological-ode-why-is-my-solver-not-sticking-to-zero/90409/25 "2022-11-20T14:31:49Z")

</div>

> [@iago-lito](#):
>
> Or is this approach “unsafe” in the _very hard sense_ that it eventually violates low-level assumptions of Julia’s compiler, which in turn invokes LLVM’s Undefined Behaviour?

I’m going to hazard a strong guess and say you don’t need to worry about this. The issue is a mathematical one not memory safety etc.

Based on what we’ve learned here I would say you should use explicit methods, and for realism also have a callback that checks to see whether any population falls below a threshold for “one organism” (however you define this in your parameterization) and sets those parameters to 0.0. this is an understandable scientifically motivated process. It’s not just that extinct species should stay extinct but also that species should go explicitly extinct if their population falls below 1 organism.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [November 21, 2022, 1:13am UTC](https://discourse.julialang.org/t/zombies-in-biological-ode-why-is-my-solver-not-sticking-to-zero/90409/26 "2022-11-21T01:13:36Z")

</div>

> [@iago-lito](#):
>
> Or is this approach “unsafe” in the _very hard sense_ that it eventually violates low-level assumptions of Julia’s compiler, which in turn invokes LLVM’s Undefined Behaviour? So mutating `B` in any situation (including ours) could corrupt low-level code optimizations and lead to _non-reproducible_ crashes of the REPL, memory corruption, unexplainable bugs or even silent wrong trajectories? If yes, then I’m afraid there is no way we should use it at all in professional software 😱

No, I never said any of this and this is simply not true. I don’t know where this idea would even come from.

> [@iago-lito](#):
>
> Is this approach “unsafe” in the _soft sense_ that it only violates high-level assumptions of `DifferentialEquations`? So mutating `B` in the general case could corrupt the solver’s internal state and lead to _reproducible_ wrong simulations or crashes of `DifferentialEquations`? If yes, then I trust you that we find ourselves in a particular situation where we can safely assume that mutating an artefactual, small positive number, resulting from only LU-factorization to zero is not actually state-breaking (at least in the current implementation) and we can go for it.

Yes, it’s soft in this sense. It violates the assumptions of any adaptive ODE solver. I show exactly this point in my JuliaCon talk about debugging simulation codes. Here’s where I give a high level overview of the adaptivity process:

[![](https://global.discourse-cdn.com/julialang/original/3X/5/d/5d8be9be2d861eecc910a9bdda15b5b546dd839f.jpeg "How to debug Julia simulation codes (ODEs, optimization, etc.!)") ](https://www.youtube.com/watch?v=g-iOOhh2U6o&t=1187)

Because you do not always go forward in time due to this rejection behavior, but mutation actually modifies the caches, this means that after rejection you will have mutated the value of `u` would could, in theory, change the behavior of the equation. In a general sense, this could lead to cases where a single `f` is not well-defined in the sense of `u' = f(u)` since the way that you have accepted/rejected steps changes what values come out of `f`. That is then precisely the issue I talk about at this point of the talk:

[![](https://global.discourse-cdn.com/julialang/original/3X/5/d/5d8be9be2d861eecc910a9bdda15b5b546dd839f.jpeg "How to debug Julia simulation codes (ODEs, optimization, etc.!)") ](https://www.youtube.com/watch?v=g-iOOhh2U6o&t=932)

where I say, such a case does not actually define an ODE, since `u' = f(u)` does not exist because it’s not `f(u)` but rather some `f(how steps are computed)` and does not have a generally unique answer (and I show a case where the ODE does not have a unique answer from ODE solvers).

However, this specific case of mutation is deterministic in the sense that for any reasonable `dt`, it will always send the value down to zero, and therefore you get the same solution whether you previously rejected a step or if the last step was an acceptance. So while it does not necessarily define a strict `u' = f(u)` ODE, it is “close enough” (for reasonable `dt`) that it then ends up defining an ODE with a unique solution in practice, which then fixes the solver behavior.

---

<div class="post-metadata">

**Author:** ![iago-lito](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/iago-lito/32/31924_2.png) [@iago-lito](https://discourse.julialang.org/u/iago-lito)\
**Post date:** [November 21, 2022, 9:00am UTC](https://discourse.julialang.org/t/zombies-in-biological-ode-why-is-my-solver-not-sticking-to-zero/90409/27 "2022-11-21T09:00:32Z")

</div>

> [@dlakelan](#):
>
> I’m going to hazard a strong guess and say you don’t need to worry about this. The issue is a mathematical one not memory safety etc.

> [@ChrisRackauckas](#):
>
> No, I never said any of this and this is simply not true. I don’t know where this idea would even come from.

Great ^ ^ Sorry I had this slight worry just because of some _“lack of safety in mutation”_ phrase used [there](https://github.com/SciML/OrdinaryDiffEq.jl/issues/1651#issuecomment-1102624067). (context: I am personally coming from a Rust background where associating the two words “mutate” and “unsafe” really means scary things ^ ^")

Happy to read that this is nothing to worry about 🙂 Thank you again! And especially for this fascinating talk \<3

[Previous page](https://discourse.julialang.org/t/zombies-in-biological-ode-why-is-my-solver-not-sticking-to-zero/90409.md?page=1)