# Attractors.basins_fractions throws "ArgumentError: range must be non-empty"

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!

I’ve found another way to trigger this error. When I use the same setup as above with 17 dimensions and dimension 1e-2 grid cells

n_species = 17 #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 get a reasonable result:

Mapping initial conditions to attractors: 100%|██████████| Time: 0:00:35
Dict{Int64, Float64} with 3 entries:
2 => 0.201
3 => 0.0146
1 => 0.7844


However, when I use 1e-3 grid cells,

n_species = 17 #number of species s
ij_strength = 1.2 #interaction strength c
SEED = 1234
Random.seed!(SEED)
...
xg = range(0, 1; length=1001) #grid dimension
grid = ntuple(_->xg, n_species)
...
fractions = basins_fractions(mapper, sss; N=10_000)


I get the error:

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{17}}, X::CartesianIndices{17, NTuple{17, 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, 17, OrdinaryDiffEq.ODEIntegrator{true, CompositeAlgorithm{Tuple{Vern7{Static.False,…}, Rodas4{0,true,…}},…}, Vector{Float64}, Nothing,…}, Vector{Array{Float64}}}, grid::NTuple{17, 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, 17, OrdinaryDiffEq.ODEIntegrator{true, CompositeAlgorithm{Tuple{Vern7{Static.False,…}, Rodas4{0,true,…}},…}, Vector{Float64}, Nothing,…}, Vector{Array{Float64}}}, grid::NTuple{17, 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, 17, OrdinaryDiffEq.ODEIntegrator{true, CompositeAlgorithm{Tuple{Vern7{Static.False,…}, Rodas4{0,true,…}},…}, Vector{Float64}, Nothing,…}, Vector{Array{Float64}}}, grid::NTuple{17, 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, 17, OrdinaryDiffEq.ODEIntegrator{true, CompositeAlgorithm{Tuple{Vern7{Static.False,…}, Rodas4{0,true,…}},…}, Vector{Float64}, Nothing,…}, Vector{Array{Float64}}}, grid::NTuple{17, 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[39]:29


Nevertheless, when I use the even finer 1e-4 grid cells,

n_species = 17 #number of species s
ij_strength = 1.2 #interaction strength c
SEED = 1234
Random.seed!(SEED)
...
xg = range(0, 1; length=10001) #grid dimension
grid = ntuple(_->xg, n_species)
...
fractions = basins_fractions(mapper, sss; N=10_000)


I get a reasonable result quantitatively different from but qualitatively similar to the one I obtained with 1e-2 grid cells:

Dict{Int64, Float64} with 3 entries:
2 => 0.111
3 => 0.0134
1 => 0.8756


Using mapper = AttractorsViaRecurrences(glv, grid; Δt=0.1) (adding the kwarg Δt=0.1) seems to fix it. Running automatic_Δt_basins(glv, grid) throws the same error. I would still like help getting some intuition on this. Thanks!

Hey!

From the stacktrace you can see that the error occurs initially in the function automatic_Δt_basins. As the name suggests, the function finds a nice value for Δt for the integration. So from this we can understand why your solution works: by defining yourself the value of Δt, the function is entirely skipped and the error avoided. This is completely fine, and a nice solution.

But to understand the error, I’ve looked at the function code and found that the error is quite strange. It boils down to just:

a = CartesianIndices((101, 101, 101, 101, 101, 101, 101, 101, 101, 101))
N = 5000
rand(a, N)


which gives ArgumentError: range must be non-empty.

I have no clue why this occurs, but it is unrelated to the Attractors.jl or related packages. @Datseris do you have any clues?

Anyways, Josh, if you can just give the Δt value you should be good to go!

1 Like