Basins of Attraction of irregularly spaced grid

[mdepitta]

How can I do a search/characterization of the basins of attraction on an irregular grid? It is, in fact, not an uncommon situation to have basins of attraction that are robust against some state variable perturbations rather than others. So it is convenient to have a search on, say, N points on variables that are known to have little impact on trajectories' ending points within a reasonable interval, but on 2xN or more points on those that are critical, and show high variability.

Ideally say, we would like to have:

    using DynamicalSystems, OrdinaryDiffEq
  
    ds = ContinuousDynamicalSystem(my_model, u0, p)
    
    ## Say "my_model" has a state space of three variables 
    
    x1 = range(0,10,length=N_1)
    x2 = range(0,70,length=N_2)
    x3 = range(0.01,5,length=N_3)
    grid = (x1,x2,x3)     


    # # Calculate basins of attraction (old)
    reltol = 1e-6
    abstol = 1e-6
    diffeq = (alg = Vern9(), reltol = reltol, abstol = abstol)
    attractors = ... # Some assignment from known attractors of my_model 
    mapper = AttractorsViaProximity(mevmo, attractors; diffeq=diffeq)
    basins, = basins_of_attraction(mapper,grid)

When I tried to run something similar to this I get in IJulia:

MethodError: no method matching generate_ic_on_grid(::Tuple{Matrix{Float64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, ::CartesianIndex{4})

Closest candidates are:
generate_ic_on_grid(::Tuple{Vararg{T, B}}, ::Any) where {B, T} at ~/.julia/packages/ChaosTools/PHPDF/src/basins/attractor_mapping.jl:145

[awage]

Hi! Once you have the mapper you can feed the initial conditions to the mapper in any order you like.

For example

    y1r = range(-1, 1.5, length = res)
    y2r = range(-0.3, 0.3, length = res)
    ics = [ u0(y1,y2) for y1 in y1r, y2 in y2r]
    bsn = [ mapper(u) for u in ics]

u0 is a function that takes two inputs and creates a Nd vector on the state space. This way I have a custom projection.

But you can define any grid you fancy and iterate over the initial conditions. If the grid is very large you can also use on the fly statistics to compute the basins volumes or whatever measure you need.

[mdepitta]

Thanks, @awage. I have a quick follow-up question on this. I noticed that it must be grid::Tuple of a specific kind to generate the ICs grid. Eventually, looking at basins_of_attraction.jl, the function responsible for generating the ICs is generate_ic_on_grid under attractor_mapping.jl. I guess my little knowledge of Julia so far prevents me from understanding how I can pass the grid tuple by mixing true vectors and StepRangeLen. Meaning: if I give this grid the basins_of_attraction(mapper, grid) will run:

    gr = range(0.01, 70; length=npts*2)
    fr = range(1.0, 4.5; length=npts*2)
    ar = range(1, 40; length=npts*2)
    cr = range(1, 250; length=npts*2)
    grid = (gr, fr, ar, cr)

But if I provide this one instead, grid_ic_on_grid will fail:

    gr = vcat(exp10.(range(log10(0.01), log10(0.74), length=npts)),exp10.(range(log10(0.75), log10(72),length=npts)))
    fr = range(1.0, 4.5; length=npts*2)
    ar = range(1, 40; length=npts*2)
    cr = range(1, 250; length=npts*2)
    grid = (gr, fr, ar, cr)

[Datseris]

Notice that Attractors.jl is a separate module now that has several updates and modifications from the version of basins/tipping/etc. currently present in the feature frozen DynamicalSystems.jl documentation. I've transferred the issue to the appropriate repo. I'd recommend using Attractors.jl directly for now.

[Datseris]

@mdepitta the problem is that the grid needs to be a tuple of AbstractRanges. Your vcat(exp10.(range(log10(0.01), log10(0.74), length=npts)),exp10.(range(log10(0.75), log10(72),length=npts))) isn't a range, it is a Vector.

From a computing perspective, there is no reason for things to be ranges. The code would work exactly the same given vectors or ranges or practically any iterable. Scientifically, I am not sure however if anything outside a range would make sense. You want to have a tesselation of the stapte space so that in the limit of box size of tesselation going to 0 you get the exact answer of the basins of attractions.

On the other hand, maybe we should allow any arbitrary vector to be given as a grid for a dimension, provided the vector is sorted? In the source code it is a rather trivial change of AbstractRange to AbstractVector. But I am not sure how much sense this would make for the recurrences algorithm....? @awage thoughts?

[Datseris]

For the future; Please, when you say "doing X fails", always also paste the actual error message, so that the devs don't have to speculate. At the moment I am speculating about your actual error.

[Datseris]

Also, hi, and welcome to Julia and the DynamicalSystems.jl framework, and thanks a lot for opening issues, this is very helpful :)

[mdepitta]

@Datseris Thanks for the welcoming words!
I see your point. However, you could still have a valid result if the non-uniform tessellation satisfies the limit (although some additional constraints need to be assumed).

For my case of non-uniform / mixed data type grid I came up with the following temporary solution: that is an edit of the estimate_by_proximity(), which was originally in basins_of_attraction (DynamicalSystems.jl version 2..3.2):

## These are my variable ranges
npts = 50
gr = vcat(exp10.(range(log10(0.01), log10(32), length=npts)),range(log10(35), log10(72),length=npts))
fr = range(1.0, 4.5; length=npts*2)
ar = exp10.(range(log10(1.0), log10(40.0), length=npts*2))
cr = exp10.(range(log10(1.0), log10(250.0), length=npts*2))

## Refined 
grid = (gr, fr, ar, cr)

# # Calculate basins of attraction (old)
reltol = 1e-6
abstol = 1e-6
diffeq = (alg = Vern9(), reltol = reltol, abstol = abstol)

## Find basins by proximity
## In my case attractors are known and I read them from data 
mapper = AttractorsViaProximity(mevmo, attractors; diffeq=diffeq)

## This is a modified version of `estimate_by_proximity`
## Allocate basins (this is a modified for-cycle from basins_of_attraction.estimate_basins_proximity to account for special ICs)
basins = zeros(Int16, map(length, grid))
progress = ProgressMeter.Progress(
    length(basins); desc = "Basins of attraction: ", dt = 1.0
)
for (k,ind) in enumerate(CartesianIndices(basins))
    ProgressMeter.update!(progress, k)
    ## ICs
    y0 = SVector{4}([gr[ind[1]],fr[ind[2]],ar[ind[3]],cr[ind[4]]])  #<------------------- This is the modified line, instead of using `generate_ic_on_grid`
    ## Associated basins to that ics
    basins[ind] = mapper(y0)
end

[Datseris]

I don't understand what basins_of_attraction.estimate_by_proximity() means. It isn't Julia code nor syntax, it reads like Python to me. Can you please clarify where you are referring to? Are you referring to the source code of Attractors.jl? If so, which function ,in which file? Are you referring to your own code? If so, which code.

(Also, please use julia after using triple ticks: ```julia, as this makes the code have highlighting)

[awage]

On the other hand, maybe we should allow any arbitrary vector to be given as a grid for a dimension, provided the vector is sorted? In the source code it is a rather trivial change of AbstractRange to AbstractVector. But I am not sure how much sense this would make for the recurrences algorithm....? @awage thoughts?

This change would not break anything and I don't see any possible regression on performances. The function _generate_ic_on_grid does exactly what its name says. If the grid is regular or not doesn't matter.