# Regime estimation for parameter space of small system

**URL:** <https://discourse.julialang.org/t/regime-estimation-for-parameter-space-of-small-system/41588>\
**Category:** Performance\
**Tags:** diffeq\
**Created:** [June 17, 2020, 12:46pm UTC](https://discourse.julialang.org/t/regime-estimation-for-parameter-space-of-small-system/41588 "2020-06-17T12:46:01Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![SidxA](https://avatars.discourse-cdn.com/v4/letter/s/7993a0/32.png) [@SidxA](https://discourse.julialang.org/u/SidxA)\
**Post date:** [June 17, 2020, 12:46pm UTC](https://discourse.julialang.org/t/regime-estimation-for-parameter-space-of-small-system/41588/1 "2020-06-17T12:46:01Z")

</div>

Hello julia community,  
I’m relatively new to julia and have a few questions regarding a dynamical system  
representing a simple ecological model of five coupled equations.  
The goal is to let the system integrate for some number of steps  
and then look at qualities such as individual variances, wether it has been terminated  
or the dominant lyapunov exponent of the last elements, repetitive for different parameters.  
Here are my first steps towards the main aims:

```julia
using DifferentialEquations
using DynamicalSystems
using Statistics
function build(p1,p2)
    function evolution!(du,u,p,t)
     du[1] = u[1]*( 2.5*(1-u[1]) - (2.5*u[2])/(1+5*(u[1]+u[2]+u[3]) ) )
     du[2] = u[2]*( 2.5*(1-u[2]) - (5*u[2])/(1+5*(u[1]+u[2]+u[3]) ) )
     du[3] = u[3]*( 2.5*(1-u[3]) - (7.5*u[2])/(1+5*(u[1]+u[2]+u[3]) ) )
     du[4] = u[4]*( ((2.5*u[1]+5*u[2]+7.5*u[3])/(1+5*(u[1]+u[2]+u[3]) )) - ((1*u[5])/(1+2*u[4])) - p[1] )
     du[5] = u[5]*( (1*u[4])/(1+2*u[4]) - p[2] )
        return nothing
    end
    function condition(u,t,integrator)
      minimum(u)-10^-4
    end
    affect!(integrator) = terminate!(integrator)
    cb = ContinuousCallback(condition,affect!)
    u0 = SVector{5}(0.5,0.5,0.5,0.5,0.5)
    ds = ContinuousDynamicalSystem(evolution!, u0, [0.1,0.1],t0 = 0.0)
    for i in p1, j in p2
        set_parameter!(ds,[i,j])
        prob = ODEProblem(ds, (0,100.0),[i,j])
        sol = solve(prob,callback=cb,dense=false)
        if sol.t[end] != 100.0
           println(sol.t[end])
        elseif var(sol[1,end-20:end]) < 10^-4 && var(sol[2,end-20:end]) < 10^-4
           println("equilibrium")
        end
        println(lyapunov(ds,100))
    end
end

```

Calling

```julia
build([0.8],[0.4])

```

only takes a couple of seconds, but allready 8\*8 different values seems never to end.  
Now there a quite a few issues I got, any help to any of them is very much appreciated, even some link if I did’nt research enough. The integration period is aimed at a lot of steps of which the most intensive beside the integration I believe is going to be the computation of the lyapunov exponent of about the last fifth. That the whole process is going to take a while to compute is okay, I just want to optimize as much as possible.

i) I had the feeling some inplace static equations had roughly the same consumption for one solve, but is this now the right choice for so many iterations?

ii) I had the feeling that the trajectory() out of DynamicalSystems.jl had a slightly different behaviour, which I’d like to test but did’nt manage to implement the terminal callback with something like

```julia
trajectory(ds,100.0, callback = cb)

```

iii) Im relatively new to numerical integration and am rather insecure about the choice of solving algorithm, is the standard Tsit5() a good choice here? Is it reasonable to try out for higher tolerances explicitely? The numbers usually stay pretty low and I believe there are no huge slopes to expect.

iv)Do I have to specify some kind of multithreading or does julia automate the computation on the different cores? Would it be good to try out some form such that the solver is done one after another using the same space and storing the evaluation somewhere else because intermediate results are not required later? Regarding that which variables are better static like the initial conditions here for instance?

v) Is there property to check whether solver got terminated maybe by the callback, easier then comparing the final time-value?

vi) Is the usage of the minimum() a good choice in the callback to define some lower threshold or would it increase the performance by putting the callback aside and looking later through the values? The point is that in some cases after the threshold is reached another equation goes through the roof ~e80.

vii) If I wanted to implement more detailed conditional computation, say check for termination and only if not for the variance of all equations, are there faster ways then simple elseif?

viii) Is the

```julia
set_parameter!()

```

a bit expendable here as I already call the parameters in the ODEproblem?

Many Thanks for your time and effort people!

---

<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:** [June 17, 2020, 1:21pm UTC](https://discourse.julialang.org/t/regime-estimation-for-parameter-space-of-small-system/41588/2 "2020-06-17T13:21:04Z")

</div>

> [@SidxA](#):
>
> i) I had the feeling some inplace static equations had roughly the same consumption for one solve, but is this now the right choice for so many iterations?

depends. It can be better just because the cache is smaller though.

> [@SidxA](#):
>
> ii) I had the feeling that the trajectory() out of DynamicalSystems.jl had a slightly different behaviour, which I’d like to test but did’nt manage to implement the terminal callback with something like
> 
> ```julia
> 
> ```

that’d be new to me.

> [@SidxA](#):
>
> iii) Im relatively new to numerical integration and am rather insecure about the choice of solving algorithm, is the standard Tsit5() a good choice here? Is it reasonable to try out for higher tolerances explicitely? The numbers usually stay pretty low and I believe there are no huge slopes to expect.

Best thing to do is just try things. There’s a bunch of suggestions in [https://diffeq.sciml.ai/latest/solvers/ode\_solve/#Recommended-Methods-1](https://diffeq.sciml.ai/latest/solvers/ode_solve/#Recommended-Methods-1), and yes try different tolerances.

> [@SidxA](#):
>
> iv)Do I have to specify some kind of multithreading or does julia automate the computation on the different cores? Would it be good to try out some form such that the solver is done one after another using the same space and storing the evaluation somewhere else because intermediate results are not required later? Regarding that which variables are better static like the initial conditions here for instance?

If using pure DiffEq here, you can use [https://diffeq.sciml.ai/latest/features/ensemble/](https://diffeq.sciml.ai/latest/features/ensemble/) which will parallelize it.

> [@SidxA](#):
>
> v) Is there property to check whether solver got terminated maybe by the callback, easier then comparing the final time-value?

`sol.retcode == :Terminated`

> [@SidxA](#):
>
> vi) Is the usage of the minimum() a good choice in the callback to define some lower threshold or would it increase the performance by putting the callback aside and looking later through the values? The point is that in some cases after the threshold is reached another equation goes through the roof ~e80.

`minimum` is fine. Yes, it sounds like the equation could have bad behavior near zero, and extinctions in ecological models are definitely an example of this.

> [@SidxA](#):
>
> ```julia
> set_parameter!()
> 
> ```
> 
> a bit expendable here as I already call the parameters in the ODEproblem?

Not sure what you mean.

---

<div class="post-metadata">

**Author:** ![SidxA](https://avatars.discourse-cdn.com/v4/letter/s/7993a0/32.png) [@SidxA](https://discourse.julialang.org/u/SidxA)\
**Post date:** [June 18, 2020, 12:37pm UTC](https://discourse.julialang.org/t/regime-estimation-for-parameter-space-of-small-system/41588/3 "2020-06-18T12:37:36Z")

</div>

Hi ChrisRackauckas,  
what a great and fast response.

ii) calling trajectory(ds,100,callback = cb) yieds

```julia
5-dimensional Dataset{Float64} with 10001 points
...
 0.938048 0.876095 0.814142 0.243634 1.26175e-12
 0.938048 0.876095 0.814142 0.243652 1.25877e-12

```

and viii)  
nevermind I need the set\_parameter!() for the lyapunov I was just curious wether I could have just called ODEProblem if that wasnt required.

So what do you believe is the most efficient way to iterate over my parameter space?

Kind regards

---

<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:** [June 18, 2020, 12:45pm UTC](https://discourse.julialang.org/t/regime-estimation-for-parameter-space-of-small-system/41588/4 "2020-06-18T12:45:15Z")

</div>

> [@SidxA](#):
>
> So what do you believe is the most efficient way to iterate over my parameter space?

[GitHub - SciML/QuasiMonteCarlo.jl: Lightweight and easy generation of quasi-Monte Carlo sequences with a ton of different methods on one API for easy parameter exploration in scientific machine learning (SciML)](https://github.com/SciML/QuasiMonteCarlo.jl) is usually a good tool for this.

---

<div class="post-metadata">

**Author:** ![SidxA](https://avatars.discourse-cdn.com/v4/letter/s/7993a0/32.png) [@SidxA](https://discourse.julialang.org/u/SidxA)\
**Post date:** [July 22, 2020, 8:36am UTC](https://discourse.julialang.org/t/regime-estimation-for-parameter-space-of-small-system/41588/5 "2020-07-22T08:36:51Z")

</div>

Good tip.  
Do you happen to know the difference between passing a constructor directly to an ODEProblem  
or giving it to a DynamicalSystem?  
Because if I increase the dimensionality of a given system

```julia
function evolute!(du,u,p,t)
 du[1] = 0.98*u[1]-0.2*u[1]*u[2]/(1+0.05*u[1])
 du[2] = -1*u[2]+0.2*u[1]*u[2]/(1+0.05*u[1])-1*u[2]*u[3]
 du[3] = -10*(u[3]-0.006)+1*u[2]*u[3]
    return nothing
end

```

to

```julia
function algeb_evol!(dr,r,p,t)
    u = @view r[:,1]
    v = @view r[:,2]
    w = @view r[:,3]
    du = @view dr[:,1]
    dv = @view dr[:,2]
    dw = @view dr[:,3]
    @. du = 0.98 * u - 0.2 * u * v / (1 + 0.05 * u)
    @. dv = - 1 * v + 0.2 * u * v / (1 + 0.05 * u) - 1 * v * w 
    @. dw = - 10 * (w - 0.006) + 1 * v * w
    dv += 0.9 .* sum((A .* v) .- transpose(A.*v),dims=2)
    dw += 0.9 .* sum((A .* w) .- transpose(A.*w),dims=2)
end

```

I can solve it directly by passing it to the solver, but I dont understand the argumenterror trying to give it to a dynamicalsystem:

```julia
ArgumentError: The state of a dynamical system *must* be <: AbstractVector/Number!

```

---

<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:** [July 22, 2020, 12:17pm UTC](https://discourse.julialang.org/t/regime-estimation-for-parameter-space-of-small-system/41588/6 "2020-07-22T12:17:26Z")

</div>

DynamicalSystems.jl only supported AbstractVector, which is something that could be fixed by them doing a `vec` on their side @Datseris

---

<div class="post-metadata">

**Author:** ![Datseris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/datseris/32/13406_2.png) [@Datseris](https://discourse.julialang.org/u/Datseris)\
**Post date:** [July 22, 2020, 12:36pm UTC](https://discourse.julialang.org/t/regime-estimation-for-parameter-space-of-small-system/41588/7 "2020-07-22T12:36:09Z")

</div>

Hi thanks for the tag, was not aware of this discussion. (@SidxA you’re more likely to get a response if you ask questions directly on GitHub, or the chat channels, at least for JuliaDynamics-related topics).

I am not sure whether the solution is as trivial as using `vec`, because the reason I explicitly ask for vectors was to play well with ForwardDiff. Maybe things have moved on since then and now putting a `vec` around everything would actually work…? Fingers crossed I guess!

Regardless, @SidxA please open a GitHub issue explicitly on DynamicalSystemsBase.jl and put there your code `algeb_evol!` (and also the full code of _exactly_ how you gave it to the `ContinuousDynamicalSystem` constructor). If this stays here I can tell you for sure I won’t remeber to do something about it…

* * *

`trajectory` is a wrapper around `solve` which explicitly passes the `saveat` argument and otherwise propagates all other arguments. (thus, the only reason to use it is as a shorthand (and this is the only reason we use it at DynamicalSystems.jl as well)).

* * *

DynamicalSystems.jl doesn’t have any multithreading, Chris’ suggestion is the way to go.

* * *

`set_parameter!` is safe (to my knowledge, if not, you have to open an issue), as it just literally sets the container values to the given ones.

* * *

> [@SidxA](#):
>
> only takes a couple of seconds, but allready 8\*8 different values seems never to end.

That is too general of a statement to know what is happening… You could use the callbacks to trigger a `println` statement, and see whether things terminate or not.

---

<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:** [July 22, 2020, 4:14pm UTC](https://discourse.julialang.org/t/regime-estimation-for-parameter-space-of-small-system/41588/8 "2020-07-22T16:14:02Z")

</div>

> [@Datseris](#):
>
> Maybe things have moved on since then and now putting a `vec` around everything would actually work…? Fingers crossed I guess!

Yes, you can take Jacobians by interpreting the high order tensor as a vector through `vec` (which is lazy)