Hi, I’m trying to find basin fractions of an s-species (s-dimensional) generalized Lotka-Volterra model using DynamicalSystems. My model is
d_t\vec{u}=\vec{\lambda}+\vec{r}\cdot\vec{u}\cdot (1-\hat{\alpha}\vec{u})
where \hat{\alpha}_{ii}=1 and \hat{\alpha}_{ij,i\neq j} are distributed c\cdot\text{Gamma}(1/\sigma^2, \sigma^2). My setup, which works fine for some parameters (more on this later), is
using DifferentialEquations
using Distributions
using DynamicalSystems
using LinearAlgebra
using Random
const SIGMA = 0.5 #sigma
const DISPERSAL = 1.0e-6 #scalar lambda
const GROWTH_RATE = 1.0 #scalar r
const MIN_POP = 1.0e-9 #threshold below which to impose extinction
# function to generate interaction matrix
function interaction_matrix(σ, s, c)
# initializes s-by-s matrix with entries ~ Gamma
α = c*rand(Gamma(1/σ^2,σ^2), (s,s))
α[diagind(α)] .= 1.0 #sets diagonals to 1
return α
end
# dynamic rule function
function GeneralizedLotkaVolterra!(du, u, p, t)
λ, r, α = p
du .= λ+r.*u.*(1 .-α*u)
return nothing
end
When I run basins_fractions() with 9 species
n_species = 9 #number of species s
ij_strength = 1.2 #interaction strength c
SEED = 1234
Random.seed!(SEED)
λ = DISPERSAL*ones(n_species) #vector lambda
r = GROWTH_RATE*ones(n_species) #vector r
α = interaction_matrix(SIGMA, n_species, ij_strength)
p = [λ, r, α]
u0 = rand(n_species) #initial conditions
xg = range(0, 1; length=101) #grid dimension
grid = ntuple(_->xg, n_species)
# checks if species fall below extinction threshold
condition(u, t, integrator) = any(u.<MIN_POP)
# sets species below extinction threshold to 0
function affect!(integrator)
integrator.u[integrator.u.<MIN_POP] .= 0
end
cb = DiscreteCallback(condition, affect!)
diffeq = (alg=AutoVern7(Rodas4()),
reltol=1.0e-8,
abstol=1.0e-8,
callback=cb)
glv = CoupledODEs(GeneralizedLotkaVolterra!, u0, p; diffeq)
mapper = AttractorsViaRecurrences(glv, grid)
sss, = statespace_sampler(;
min_bounds=minimum.(grid),
max_bounds=maximum.(grid)
)
fractions = basins_fractions(mapper, sss; N=10_000)
I quickly get a reasonable result
Mapping initial conditions to attractors: 100%|██████████| Time: 0:00:13
Dict{Int64, Float64} with 3 entries:
2 => 0.0914
3 => 0.1504
1 => 0.7582
However, changing only the first line to use 10 species,
n_species = 10 #number of species s
ij_strength = 1.2 #interaction strength c
SEED = 1234
Random.seed!(SEED)
...
fractions = basins_fractions(mapper, sss; N=10_000)
I instead get this error message I’m having trouble understanding
ArgumentError: range must be non-empty
Stacktrace:
[1] SamplerRangeNDL
@ /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/Random/src/generation.jl:333 [inlined]
[2] Sampler
@ /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/Random/src/generation.jl:189 [inlined]
[3] Sampler
@ /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/Random/src/generation.jl:418 [inlined]
[4] Sampler (repeats 2 times)
@ /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/Random/src/Random.jl:140 [inlined]
[5] rand!(rng::TaskLocalRNG, A::Vector{CartesianIndex{10}}, X::CartesianIndices{10, NTuple{10, Base.OneTo{Int64}}})
@ Random /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/Random/src/Random.jl:267
[6] rand
@ /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/Random/src/Random.jl:280 [inlined]
[7] rand
@ /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/Random/src/Random.jl:281 [inlined]
[8] rand
@ /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/Random/src/Random.jl:284 [inlined]
[9] automatic_Δt_basins(ds::CoupledODEs{true, 10, OrdinaryDiffEq.ODEIntegrator{true, CompositeAlgorithm{Tuple{Vern7{Static.False,…}, Rodas4{0,true,…}},…}, Vector{Float64}, Nothing,…}, Vector{Array{Float64}}}, grid::NTuple{10, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}; N::Int64)
@ Attractors ~/.julia/packages/Attractors/Y7iTL/src/mapping/attractor_mapping_recurrences.jl:280
[10] automatic_Δt_basins
@ ~/.julia/packages/Attractors/Y7iTL/src/mapping/attractor_mapping_recurrences.jl:271 [inlined]
[11] initialize_basin_info(ds::CoupledODEs{true, 10, OrdinaryDiffEq.ODEIntegrator{true, CompositeAlgorithm{Tuple{Vern7{Static.False,…}, Rodas4{0,true,…}},…}, Vector{Float64}, Nothing,…}, Vector{Array{Float64}}}, grid::NTuple{10, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Δtt::Nothing, sparse::Bool)
@ Attractors ~/.julia/packages/Attractors/Y7iTL/src/mapping/attractor_mapping_recurrences.jl:219
[12] AttractorsViaRecurrences(ds::CoupledODEs{true, 10, OrdinaryDiffEq.ODEIntegrator{true, CompositeAlgorithm{Tuple{Vern7{Static.False,…}, Rodas4{0,true,…}},…}, Vector{Float64}, Nothing,…}, Vector{Array{Float64}}}, grid::NTuple{10, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}; Δt::Nothing, sparse::Bool, force_non_adaptive::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ Attractors ~/.julia/packages/Attractors/Y7iTL/src/mapping/attractor_mapping_recurrences.jl:121
[13] AttractorsViaRecurrences(ds::CoupledODEs{true, 10, OrdinaryDiffEq.ODEIntegrator{true, CompositeAlgorithm{Tuple{Vern7{Static.False,…}, Rodas4{0,true,…}},…}, Vector{Float64}, Nothing,…}, Vector{Array{Float64}}}, grid::NTuple{10, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}})
@ Attractors ~/.julia/packages/Attractors/Y7iTL/src/mapping/attractor_mapping_recurrences.jl:118
[14] top-level scope
@ In[4]:29
Surprisingly, 12 species works just fine, and
n_species = 12 #number of species s
ij_strength = 1.2 #interaction strength c
SEED = 1234
Random.seed!(SEED)
...
fractions = basins_fractions(mapper, sss; N=10_000)
gives
Mapping initial conditions to attractors: 100%|██████████| Time: 0:00:21
Dict{Int64, Float64} with 2 entries:
2 => 0.0022
1 => 0.9978
Again with 10 species, changing the seed from 1234 to 1243
n_species = 10 #number of species s
ij_strength = 1.2 #interaction strength c
SEED = 1243
Random.seed!(SEED)
...
fractions = basins_fractions(mapper, sss; N=10_000)
throws the same error. Similar errors occur for 11 and 14 dimensions but not 13 or 15.
How might I fix this? Thanks!
please help, thanks!