# PositiveDomain callback and isoutofdomain implementation

**URL:** <https://discourse.julialang.org/t/positivedomain-callback-and-isoutofdomain-implementation/11097>\
**Category:** New to Julia\
**Tags:** diffeq\
**Created:** [May 23, 2018, 2:28pm UTC](https://discourse.julialang.org/t/positivedomain-callback-and-isoutofdomain-implementation/11097 "2018-05-23T14:28:55Z")\
**Posts on this page:** 18\
**Page:** 1

<div class="post-metadata">

**Author:** ![swainarpit](https://avatars.discourse-cdn.com/v4/letter/s/7feea3/32.png) [@swainarpit](https://discourse.julialang.org/u/swainarpit)\
**Post date:** [May 23, 2018, 2:28pm UTC](https://discourse.julialang.org/t/positivedomain-callback-and-isoutofdomain-implementation/11097/1 "2018-05-23T14:28:55Z")

</div>

Hi,

I have been trying to solve a system of differential equations in Julia. I have run into a few problems and am not clear on what I am doing wrong.

Problem definition:  
The system I am working on is a model of the immune system. This is represented by 9 differential equations. One of the parameters is time-dependent and is a step function representing the lifetime of an infection in the immune system. This makes the system discontinuous at the time of entry and removal of the infection. To address this, I used the options tstops and callbacks in the solver. This system seems to be stiff and so I used the Rosenbrock23 algorithm to solve this.

Problems:  
Because I do not know the parameters, they are assigned random values. For independent simulations (for different parameter sets), I sometimes run into negative solutions. I tried using the callback PositiveDomain but that does not help. Even for simulations where Rosenbrock23 algorithm handles the integration quite well, adding the callback PositiveDomain gives negative solution and aborts with the warning:

WARNING: dt \<= dtmin. Aborting. If you would like to force continuation with dt=dtmin, set force\_dtmin=true

I understand that the step size can’t be reduced anymore and that’s why it’s aborting but I don’t understand why adding the callback PositiveDomain would do this. Or is my implementation faulty?  
In cases, where I do get negative values for some variables, isoutofbound solver option also doesn’t seem to make any difference.

What am I doing wrong here?

The code: [Pastiebin.com 5b0560a1b8e93](https://www.pastiebin.com/5b0560a1b8e93)

I am using DifferentialEquations version 4.5.0

---

<div class="post-metadata">

**Author:** ![mauro3](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mauro3/32/292_2.png) [@mauro3](https://discourse.julialang.org/u/mauro3)\
**Post date:** [May 23, 2018, 2:38pm UTC](https://discourse.julialang.org/t/positivedomain-callback-and-isoutofdomain-implementation/11097/2 "2018-05-23T14:38:26Z")

</div>

It’s not quite clear to me but I suspect that the differential equations allow the solution becoming negative (for some parameters). If that is indeed the case, then `PostiveDomain` callback cannot help. As you seem to be fitting parameters, I’d suggest that you have a callback which triggers as a variable reaches zero, stop the integration there, and discard the parameters as unsuitable.

---

<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:** [May 23, 2018, 2:40pm UTC](https://discourse.julialang.org/t/positivedomain-callback-and-isoutofdomain-implementation/11097/3 "2018-05-23T14:40:29Z")

</div>

> [@mauro3](#):
>
> It’s not quite clear to me but I suspect that the differential equations allow the solution becoming negative (for some parameters).

That’s what I suspect too. `alg` is specified incorrectly, but even with that fixed there’s an issue. Interestingly

```julia
sol=solve(prob,Rosenbrock23(),abstol=1e-12,reltol=1e-12,tstops=dis,callback=cb)

```

no matter how low I set the tolerance, the ODE was going negative for me at exactly `t=320`, which is exactly when the parameters are changed:

```julia
function affect!(integrator)
  @show "parameter change",integrator.t
    k = find(x->x==integrator.t,h1);
    for i = 1:length(k)
        integrator.p[10][k[i]] = rand();
    end
    k = find(x->x==integrator.t,h2);
    for i = 1:length(k)
        integrator.p[10][k[i]] = 0.0;
    end
end

```

```julia
("parameter change", integrator.t) = ("parameter change", 320.0)

```

---

<div class="post-metadata">

**Author:** ![swainarpit](https://avatars.discourse-cdn.com/v4/letter/s/7feea3/32.png) [@swainarpit](https://discourse.julialang.org/u/swainarpit)\
**Post date:** [May 23, 2018, 2:44pm UTC](https://discourse.julialang.org/t/positivedomain-callback-and-isoutofdomain-implementation/11097/4 "2018-05-23T14:44:54Z")

</div>

Hi, thanks for the reply.

I am not fitting the parameters here. I have several clones of a single type of cell and they all vary from each other in their properties. In the function f, the loop i creates n different clones. That is why I have defined random parameter values within a range before hand and am not changing it during the simulation.

---

<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:** [May 23, 2018, 2:48pm UTC](https://discourse.julialang.org/t/positivedomain-callback-and-isoutofdomain-implementation/11097/5 "2018-05-23T14:48:34Z")

</div>

Nope, it looks like there’s something up with the `PositiveDomain` callback.

```julia
sol=solve(prob,Rosenbrock23(),tstops=dis,callback=antigen)
sol=solve(prob,Rosenbrock23(),tstops=dis,isoutofdomain=(y,p,t)->any(x->x<0,y),callback=antigen)

```

These are fine.

---

<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:** [May 23, 2018, 2:54pm UTC](https://discourse.julialang.org/t/positivedomain-callback-and-isoutofdomain-implementation/11097/6 "2018-05-23T14:54:08Z")

</div>

I opened an issue to track this for PositiveDomain. In the meantime, `isoutofdomain` seems to be a good option.

[https://github.com/JuliaDiffEq/DiffEqCallbacks.jl/issues/44](https://github.com/JuliaDiffEq/DiffEqCallbacks.jl/issues/44)

---

<div class="post-metadata">

**Author:** ![swainarpit](https://avatars.discourse-cdn.com/v4/letter/s/7feea3/32.png) [@swainarpit](https://discourse.julialang.org/u/swainarpit)\
**Post date:** [May 23, 2018, 2:55pm UTC](https://discourse.julialang.org/t/positivedomain-callback-and-isoutofdomain-implementation/11097/7 "2018-05-23T14:55:06Z")

</div>

> [@ChrisRackauckas](#):
>
> `alg` is specified incorrectly

Ah! But, it still worked?

> [@ChrisRackauckas](#):
>
> These are fine

Is isoutofdomain working? I simulated a parameter set with just the Rosenbrock23 and found negative solution. So, I kept the parameter set fixed and added the isoutofdomain, yet it did not seem to work. It was still going to the negative values.

---

<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:** [May 23, 2018, 2:56pm UTC](https://discourse.julialang.org/t/positivedomain-callback-and-isoutofdomain-implementation/11097/8 "2018-05-23T14:56:48Z")

</div>

> [@swainarpit](#):
>
> I simulated a parameter set with just the Rosenbrock23 and found negative solution. So, I kept the parameter set fixed and added the isoutofdomain, yet it did not seem to work. It was still going to the negative values.

Can you post an MWE where that is happening? It didn’t go negative for the parameters I checked. A minimal working example should take out the randomness. A good check is when `sum(sol .< 0)`.

> [@swainarpit](#):
>
> But, it still worked?

Yes, this is something I want to address sooner rather than later.

> <https://github.com/SciML/DifferentialEquations.jl/issues/275>
>
> @dpsanders wants a better system for warning users about type choices (specifica…lly, when integers are used with adaptive integrators). We also want a way to automatically convert functions to input vectors https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/251 . @tkf suggested that we could do \`init\` in DiffEqBase and then have that lower to a \`\_init\`. Also, \`solve\` can be declared at a higher level as well and only overridden as necessary since it's just generic and uses \`init\` internally:
> 
> https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/solve.jl#L1-L9
> 
> This can be implemented by doing an algorithm-generic \`init(prob::...,alg::DEAlgorithm)\`. Implementations just opt in by dispatching on \`\_init\` instead of \`init\` itself. This can be implemented without being breaking at all, since this just won't be used until the solver packages themselves make the change to \`\_init\` and/or \`\_solve\`.

---

<div class="post-metadata">

**Author:** ![swainarpit](https://avatars.discourse-cdn.com/v4/letter/s/7feea3/32.png) [@swainarpit](https://discourse.julialang.org/u/swainarpit)\
**Post date:** [May 23, 2018, 3:14pm UTC](https://discourse.julialang.org/t/positivedomain-callback-and-isoutofdomain-implementation/11097/9 "2018-05-23T15:14:42Z")

</div>

> [@ChrisRackauckas](#):
>
> Can you post an MWE where that is happening?

My fault. It was because of the incorrect use of alg definition. Replacing it with Rosenbrock23() corrected it. That means I was simulating with Tsit5() till now?

---

<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:** [May 23, 2018, 3:22pm UTC](https://discourse.julialang.org/t/positivedomain-callback-and-isoutofdomain-implementation/11097/10 "2018-05-23T15:22:45Z")

</div>

> [@swainarpit](#):
>
> Replacing it with Rosenbrock23() corrected it.

That shouldn’t make a difference, and I cannot recreate that. Can you get an example where the auto-alg is getting a negative? That should just be impossible with `isoutofdomain`.

```julia
sol=solve(prob,tstops=dis,isoutofdomain=(y,p,t)->any(x->x<0,y),callback=antigen)
sum(sol .< 0)

```

> [@swainarpit](#):
>
> That means I was simulating with Tsit5() till now?

It’s hard to know since it’s dependent on a few things. My guess would be is it would’ve been using `AutoTsit5(Rosenbrock23())`.

---

<div class="post-metadata">

**Author:** ![swainarpit](https://avatars.discourse-cdn.com/v4/letter/s/7feea3/32.png) [@swainarpit](https://discourse.julialang.org/u/swainarpit)\
**Post date:** [May 23, 2018, 3:50pm UTC](https://discourse.julialang.org/t/positivedomain-callback-and-isoutofdomain-implementation/11097/11 "2018-05-23T15:50:42Z")

</div>

> [@ChrisRackauckas](#):
>
> That shouldn’t make a difference, and I cannot recreate that.

The attached code has a parameter set where the alg specification (alg=Rosenbrock23 vs Rosenbrock23()) is the difference between positive and negative solutions.

> [@ChrisRackauckas](#):
>
> Can you get an example where the auto-alg is getting a negative? That should just be impossible with `isoutofdomain` .

I will check but I don’t think I will. The command seems quite clear, I was only doubtful of the way I had implemented it. If that is correct, I doubt I can find a violation case. But if I do, I will attach a link here.

Code: [https://www.pastiebin.com/5b058c394dfc8](https://www.pastiebin.com/5b058c394dfc8)

---

<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:** [May 23, 2018, 4:27pm UTC](https://discourse.julialang.org/t/positivedomain-callback-and-isoutofdomain-implementation/11097/12 "2018-05-23T16:27:56Z")

</div>

Oh I see what it is. Checking

```julia
plot(sol,denseplot=false)

```

it seems that the interpolation goes negative (and is somewhat unstable) with the RK method, but the values are all positive. This can happen with `isoutofdomain` but cannot happen with `PositiveDomain` since `isoutofdomain` only checks the steps. In fact, this small dip may be what’s tripping up `PositiveDomain()` since it tries to maximize steps and then interpolate back pre-zero to speed up (vs `isoutofdomain` which just rejects any step outside of the domain)

But you can check that lower tolerances do make the “full interpolation” positive:

```julia
sol=solve(prob,Rosenbrock23(),abstol=1e-12,reltol=1e-12,tstops=dis,callback=antigen)
A = sol(tspan[1]:0.01:tspan[2])
sum(A.<0)

```

---

<div class="post-metadata">

**Author:** ![swainarpit](https://avatars.discourse-cdn.com/v4/letter/s/7feea3/32.png) [@swainarpit](https://discourse.julialang.org/u/swainarpit)\
**Post date:** [May 24, 2018, 2:07pm UTC](https://discourse.julialang.org/t/positivedomain-callback-and-isoutofdomain-implementation/11097/13 "2018-05-24T14:07:05Z")

</div>

> [@ChrisRackauckas](#):
>
> It’s hard to know since it’s dependent on a few things. My guess would be is it would’ve been using `AutoTsit5(Rosenbrock23())` .

I checked again. Because alg=Rosenbrock23 is the wrong syntax, the solver takes the default option as when no algorithm is specified. While using the default tolerances, the composite algorithm used for this case was Tsit5() and Rosenbrock23(), whereas by lowering the abstol and reltol to 1e-12, the solver chose to use Vern9() and Rodas5() as the composite algorithm.

> [@ChrisRackauckas](#):
>
> it seems that the interpolation goes negative (and is somewhat unstable) with the RK method, but the values are all positive.

Ah! The interpolation. I tried out Rosenbrock23(), AutoTsit5(Rosenbrock23()) and the default (i.e. the composite of Vern9() and Rodas5()) with and without lower tolerances (1e-12) and found that the interpolation in cases where Tsit5() is used with the default tolerance goes negative. So, yes, I think you are right when you say the interpolation with the RK method is at fault here.

Thank you for the help!

What can be done for the PositiveDomain()?

---

<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:** [May 24, 2018, 2:25pm UTC](https://discourse.julialang.org/t/positivedomain-callback-and-isoutofdomain-implementation/11097/14 "2018-05-24T14:25:35Z")

</div>

> [@swainarpit](#):
>
> What can be done for the PositiveDomain()?

It can get fixed for this example or we can diagnose why it’s properly failing. I am not sure which of the two it will end up being.

---

<div class="post-metadata">

**Author:** ![swainarpit](https://avatars.discourse-cdn.com/v4/letter/s/7feea3/32.png) [@swainarpit](https://discourse.julialang.org/u/swainarpit)\
**Post date:** [May 24, 2018, 2:35pm UTC](https://discourse.julialang.org/t/positivedomain-callback-and-isoutofdomain-implementation/11097/15 "2018-05-24T14:35:39Z")

</div>

I think a more universal solution would be the way to go. But fixing this problem might be a start. My negativity question has already been answered, I was using the wrong syntax! So, I will accept your last reply as the answer.

Also, I was interested in the performance of this code as I was looking to simulate a larger system (~2000 ODEs) and as it is now, the simulation fails because of memory. Do I need to start a new topic for this?

---

<div class="post-metadata">

**Author:** ![mauro3](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mauro3/32/292_2.png) [@mauro3](https://discourse.julialang.org/u/mauro3)\
**Post date:** [May 24, 2018, 2:38pm UTC](https://discourse.julialang.org/t/positivedomain-callback-and-isoutofdomain-implementation/11097/16 "2018-05-24T14:38:28Z")

</div>

> [@swainarpit](#):
>
> Do I need to start a new topic for this?

I’d say yes.

---

<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:** [May 24, 2018, 2:53pm UTC](https://discourse.julialang.org/t/positivedomain-callback-and-isoutofdomain-implementation/11097/17 "2018-05-24T14:53:10Z")

</div>

> [@swainarpit](#):
>
> Also, I was interested in the performance of this code as I was looking to simulate a larger system (~2000 ODEs) and as it is now, the simulation fails because of memory. Do I need to start a new topic for this?

Yup, new topic but feel free to ping me. My first guess would be that you’re saving the whole continuous solution which may not be a good idea in for 2000 x 2000.

---

<div class="post-metadata">

**Author:** ![swainarpit](https://avatars.discourse-cdn.com/v4/letter/s/7feea3/32.png) [@swainarpit](https://discourse.julialang.org/u/swainarpit)\
**Post date:** [May 29, 2018, 10:40am UTC](https://discourse.julialang.org/t/positivedomain-callback-and-isoutofdomain-implementation/11097/18 "2018-05-29T10:40:17Z")

</div>

For whoever is following this thread and is interested in memory management:

> [@Efficient memory allocation for large system of ODEs](https://discourse.julialang.org/t/efficient-memory-allocation-for-large-system-of-odes/11226):
>
> Hi, I have a system of differential equations that mimic the clonal population in an immune system. For whoever is curious, the problem definition can be found in the post: At present, I am simulating with a small number of clones and thus, am able to deal with the memory required. But, for higher number of clones, the performance of my present code is pretty detrimental. Say, for 2000 clones (variable n), the code runs for about 15 hours and stops prematurely because of the shortage of memo…