# Callback function to set negative values to zero

**URL:** <https://discourse.julialang.org/t/callback-function-to-set-negative-values-to-zero/111902>\
**Category:** Modelling & Simulations\
**Created:** [March 20, 2024, 8:49pm UTC](https://discourse.julialang.org/t/callback-function-to-set-negative-values-to-zero/111902 "2024-03-20T20:49:06Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![Nick1](https://avatars.discourse-cdn.com/v4/letter/n/47e85d/32.png) [@Nick1](https://discourse.julialang.org/u/Nick1)\
**Post date:** [March 20, 2024, 8:49pm UTC](https://discourse.julialang.org/t/callback-function-to-set-negative-values-to-zero/111902/1 "2024-03-20T20:49:06Z")

</div>

I have an ODE system of 60 equations and run multiple simulations with different sets of parameters. I want to prevent the outputs from taking any negative values since it is not physically plausible. I currently use something like that, based on this topic: [Setting really small values in ODE to zero](https://discourse.julialang.org/t/setting-really-small-values-in-ode-to-zero/45693)

```julia
condition0(u,t,integrator) = any( (u .< 0.0))
function affect0!(integrator)
   integrator.u[(integrator.u .< 0.0)] .= 0;
end
cb0 = DiscreteCallback(condition0,affect0!)

```

However, that way the computational time increases significantly.

**Without using the callback function**

> 0.831594 seconds (845.83 k allocations: 51.515 MiB, 8.17% gc time, 92.77% compilation time: 25% of which was recompilation)

**Using the callback function**

> 16.886240 seconds (40.82 M allocations: 2.280 GiB, 16.01% gc time, 11.04% compilation time: 11% of which was recompilation)

What would be a better way to tackle this issue? I have around 160 parameters so it is difficult to know which parts of the space cause these irregularities.

Thanks!

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 21, 2024, 12:21am UTC](https://discourse.julialang.org/t/callback-function-to-set-negative-values-to-zero/111902/2 "2024-03-21T00:21:48Z")

</div>

> [@Nick1](#):
>
> What would be a better way to tackle this issue?

I’m just guessing here, but could you add a [barrier function](https://en.wikipedia.org/wiki/Barrier_function) that prevents the solution from going negative?

e.g. instead of du/dt = f(u,t), use

\frac{du}{dt} = f(u,t) + \mathrm{barrier}.(u)

where `barrier.(u)` is an elementwise application of something like `barrier(x) = μ/x^2` (for some strength `μ`) that pushes the solutions away from zero?

Or alternatively multiply `f(u,t) * barrier.(u)` for some function like `barrier(x) = exp(-μ/x)` that goes smoothly to zero as u \to 0^+ but \to 1 for large u.

That is, if negative `u` are “not physically plausible”, add something to your physics that enforces this.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 21, 2024, 12:31am UTC](https://discourse.julialang.org/t/callback-function-to-set-negative-values-to-zero/111902/3 "2024-03-21T00:31:13Z")

</div>

Shampine et al. have a [2005 paper on “Non-negative solutions of ODEs”](https://doi.org/10.1016/j.amc.2004.12.011) that discusses some techniques.

One important point is that it is isn’t enough to simply set `u = 0` when it goes non-positive — you should also make sure your right-hand side is 0 at this point, so that the equations don’t try to keep pushing `u` to the negative regime. i.e. your right-hand side `f(u,t)` should be replaced with

```julia
F = f(u,t)
return @. ifelse(u <= 0, max(0, F), F)

```

Of course, this introduces a discontinuity into your right-hand-side, but with continuous callbacks this is manageable I guess. I haven’t read Shampine’s paper closely, but it seems to be about refining the handling of this discontinuity.

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [March 21, 2024, 12:36am UTC](https://discourse.julialang.org/t/callback-function-to-set-negative-values-to-zero/111902/4 "2024-03-21T00:36:52Z")

</div>

Somewhat equivalent to adding a barrier function is reparametrizing variable to make negative `u` impossible. For example u = e^w and finding the appropriate new form for the differential equation in w.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 21, 2024, 12:43am UTC](https://discourse.julialang.org/t/callback-function-to-set-negative-values-to-zero/111902/5 "2024-03-21T00:43:27Z")

</div>

> [@Dan](#):
>
> For example u = e^w and finding the appropriate new form for the differential equation in www.

That turns du/dt = f(u,t) into dw/dt = f(e^w, t) / e^w for a single variable (and the corresponding elementwise thing for multiple variables).

It looks like you run into a problem, however, if the system really “wants” to go negative: that is, if f(u,t) \< 0 as u \to 0^+. In this case, your w equations will have a divergence where it wants to go to w \to -\infty faster and faster (the right-hand side will diverge exponentially to -\infty as w gets more negative).

I suspect that this diverging solution will cause problems with your numerical ODE solver.

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [March 21, 2024, 12:45am UTC](https://discourse.julialang.org/t/callback-function-to-set-negative-values-to-zero/111902/6 "2024-03-21T00:45:46Z")

</div>

Doing both a barrier function and reparametrization might be a good solution. The barrier keeping the solution from -\infty and the reparametrization preventing discrete/edge effects from pulling u negative.

---

<div class="post-metadata">

**Author:** ![kkakosim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kkakosim/32/21546_2.png) [@kkakosim](https://discourse.julialang.org/u/kkakosim)\
**Post date:** [March 21, 2024, 5:17am UTC](https://discourse.julialang.org/t/callback-function-to-set-negative-values-to-zero/111902/7 "2024-03-21T05:17:40Z")

</div>

Have you tried implementing the [PositiveDomain()](https://docs.sciml.ai/DiffEqCallbacks/stable/step_control/#DiffEqCallbacks.PositiveDomain) in [DiffEqCallbacks.jl](https://docs.sciml.ai/DiffEqCallbacks/stable/), it was suggested to me for a [similar problem](https://discourse.julialang.org/t/zero-order-extinction-reaction/108749/2), and it appears faster than the “manual” callback-integrator.

However, I noted that my system is indeed keep pushing to negative, and while the u stays close to 0, there a small negative perturbations that cause other u to exceed other physical limits (e.g., exceed mass conservation). I have not tried yet the combinations mentioned earlier by @Dan and @stevengj.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 21, 2024, 12:24pm UTC](https://discourse.julialang.org/t/callback-function-to-set-negative-values-to-zero/111902/8 "2024-03-21T12:24:38Z")

</div>

> [@kkakosim](#):
>
> Have you tried implementing the [PositiveDomain()](https://docs.sciml.ai/DiffEqCallbacks/stable/step_control/#DiffEqCallbacks.PositiveDomain) in [DiffEqCallbacks.jl](https://docs.sciml.ai/DiffEqCallbacks/stable/), it was suggested to me for a [similar problem](https://discourse.julialang.org/t/zero-order-extinction-reaction/108749/2),

Oh, brilliant, that implements exactly the Shampine method I linked! It also mentions the \mathrm{ifelse}(u \le 0, \max(0, f(u,t)), f(u,t)) modification.

---

<div class="post-metadata">

**Author:** ![Nick1](https://avatars.discourse-cdn.com/v4/letter/n/47e85d/32.png) [@Nick1](https://discourse.julialang.org/u/Nick1)\
**Post date:** [March 26, 2024, 11:52pm UTC](https://discourse.julialang.org/t/callback-function-to-set-negative-values-to-zero/111902/9 "2024-03-26T23:52:33Z")

</div>

Thanks @kkakosim! That was helpful!