Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Basins of Attraction of irregularly spaced grid #35

Open
mdepitta opened this issue Jan 24, 2023 · 9 comments
Open

Basins of Attraction of irregularly spaced grid #35

mdepitta opened this issue Jan 24, 2023 · 9 comments

Comments

@mdepitta
Copy link

mdepitta commented Jan 24, 2023

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
Copy link
Contributor

awage commented Jan 25, 2023

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
Copy link
Author

mdepitta commented Jan 26, 2023

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 Datseris transferred this issue from JuliaDynamics/DynamicalSystems.jl Jan 27, 2023
@Datseris
Copy link
Member

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
Copy link
Member

Datseris commented Jan 27, 2023

@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
Copy link
Member

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
Copy link
Member

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

@mdepitta
Copy link
Author

mdepitta commented Jan 29, 2023

@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
Copy link
Member

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
Copy link
Contributor

awage commented Jan 30, 2023

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.

Datseris added a commit that referenced this issue Aug 15, 2023
* One issue in function store_attractor! to be discussed, rest works quite well, both for regular and irregular

* tests

* changelog

* typo fix

* mfs test fail fix

* Update src/mapping/attractor_mapping_recurrences.jl

Co-authored-by: George Datseris <[email protected]>

* sketch of example

* sketch of example

* Update docs/src/examples.md

Co-authored-by: George Datseris <[email protected]>

* fixes

* typo

* create an example that showcases everything correctly

---------

Co-authored-by: George Datseris <[email protected]>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants