Regime estimation for parameter space of small system

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:

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

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

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

set_parameter!()

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

Many Thanks for your time and effort people!

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

that’d be new to me.

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, and yes try different tolerances.

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

sol.retcode == :Terminated

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.

Not sure what you mean.

Hi ChrisRackauckas,
what a great and fast response.

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

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

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) is usually a good tool for this.

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

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

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:

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

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

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.


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.

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