Add a new mode for basins_of_attraction

[awage]

This PR add a new mode to compute the basins of attraction.

The idea is to get the structure BasinsInfo and to identify individual initial conditions. The algorithm can still identify new attractors and the structure holds all the important information. The basins are stored into a SparseArray since we may not need to compute all the basins. Changes are minimals and do not break anything. However there is maybe an issue with an AbstractArray declaration (I pointed this out in the code).

The function basins_of_attraction takes care of all the initialization of the system. The usage is very simple:

ds = Systems.henon_iip(zeros(2); a = 1.4, b = 0.3)
xg = yg = range(-2.,2.,length = 70)
grid = (xg,yg)
bsn_nfo, integ = basins_of_attraction(grid, ds; tracking_mode = true)
att = ChaosTools.get_color_point!(bsn_nfo, integ, [xg[1], yg[1]])

This is not intended for end-users but rather it to identify attractors in #219. This is experimental and we should test if it brings any benefit to the bStab method. @KalelR, I will check out the branch and see how it can work with #219.

[Datseris]

@awage I don't understand this PR yet... can you please give a high level description of what you are trying to do, like what changes from e.g., what we say in the paper? And how does the code work?

Here are some further comments and/or questions:

---

The idea is to get the structure BasinsInfo and to identify individual initial conditions.

So we still have an underlying grid to identify attractors via recurrences (Section II-A in our paper), and after we identify attractors we can feed any initial condition through the main loop of our algorithm. Yes, that sounds totally fine.

---

I am assuming that the function att = ChaosTools.get_color_point!(bsn_nfo, integ, [xg[1], yg[1]]) gets the index of the attractor from the initial condition [xg, yg].

---

Why do we need SparseArrays, and how are they used? You are aware that Julia only has support for sparse matrices, not higher dimensional arrays, right?

---

The naming system of the current code really makes it difficult to read and understand if one only knows our paper. We should really be sticking to the same naming conventions we use in our paper, throughout the code. This means for example that:

• We shouldn't be using words such as "color" or "draw", as we are no longer thinking in terms of drawing a basin as the initial work of Yorke. In the paper we talk about the "label" of an attractor, and we "estimate" the basins (effortlessly as well ;) )
• Tracking is a name we don't use in the paper. Is it a completely different concept?

---

I see in your code example that you already call the full basins_of_attraction function with a "tracking mode". So I guess this "tracking mode" prepares some kind of structure

---

A totally different operation of a function that is decided by a keyword argument is not good API design I think. Whatever we do in this PR, it should be a new function name. Its goal its apparently not to calculate the basis of attraction of a system, but rather to prepare some kind of "structure" that maps initial conditions to attractors, right?

---

Final question: So I guess the benefit of all this approach is that you can use a much sparser grid to create this "tracking mode", which doesn't have as heavy dimensional limitations as the normally dense grid. But, I'm not sure if this really saves a lot. Because even in 10 dimensions, 21 points along each dimension are still way too much, yet it is barely dense enough to separate attractors in state space (or maybe it isn't I actually don't know).

[awage]

First the high level description:

The functionality that I propose to include in the library is to match a single initial condition with an attractor. The mode of operation is identical to the description of the paper, there is no change here. The difference is that the function gives back to the user the structure BasinInfo and the integrator object. With these two object you can call the low level functions of basins_of_attraction that identify and match the attractors.

I have several ideas for the use of this mode:

• You can evaluate initial conditions at arbitrary points (the grid for ic's and the grid for identification are now separate things).
• Basins fraction estimation based on a sampling of the basins (maybe random walks are worth trying).
• Basin entropy estimation with random sampling
• Identification of stable and unstable manifolds of saddles.

I had to change some details of the algorithm such as the numbering of unidentified cells from 1 to 0. It does not change anything else to the code and it has no impact on the results. The reason for this change is the possible use of a sparse array that are initialized with zeros instead of ones. BUT, I was not aware that SparseArrays only support matrices. It is a big bummer for me, I had hopes to be able to sample high dimensional state space. In this case there is absolutely no need for a sparse matrix and this part should be rolled back.

The code works in the same way as basins_of_attraction does. But the user has to iterate the initial condition "manually".

I will answer the other comments later in the day.

[awage]

So we still have an underlying grid to identify attractors via recurrences (Section II-A in our paper), and after we identify attractors we can feed any initial condition through the main loop of our algorithm. Yes, that sounds totally fine.

This is the main idea of this functionality.

I am assuming that the function att = ChaosTools.get_color_point!(bsn_nfo, integ, [xg[1], yg[1]]) gets the index of the attractor from the initial condition [xg, yg].

This is the index of the basin or attractor if the initial condition is in a attractor cell.

Why do we need SparseArrays, and how are they used? You are aware that Julia only has support for sparse matrices, not higher dimensional arrays, right?

I answered earlier. It seems that SparseArrays are not efficient in high dimension and use of Dictionary of keys is an alternative. There are some package that aim at this functionallity (like SimpleSparseArrays) but it is not in the standard library and rather experimental. I will leave this part out and leave it as a side project. bSTAB seems a better tool for high dims systems.

The naming system of the current code really makes it difficult to read and understand if one only knows our paper. We should really be sticking to the same naming conventions we use in our paper, throughout the code. This means for example that:

• We shouldn't be using words such as "color" or "draw", as we are no longer thinking in terms of drawing a basin as the initial work of Yorke. In the paper we talk about the "label" of an attractor, and we "estimate" the basins (effortlessly as well ;) )

Totally agree I will rename functions and variables accordingly.

* Tracking is a name we don't use in the paper. Is it a completely different concept?

This is a bad name. I lacked of imagination the other day. ic_labeling or ic_identification fits better.

I see in your code example that you already call the full basins_of_attraction function with a "tracking mode". So I guess this "tracking mode" prepares some kind of structure

The returned structure is the same than basins_of_attraction.

A totally different operation of a function that is decided by a keyword argument is not good API design I think. Whatever we do in this PR, it should be a new function name. Its goal its apparently not to calculate the basis of attraction of a system, but rather to prepare some kind of "structure" that maps initial conditions to attractors, right?

Totally agree. I will prepare a wrapper around basins of attraction with the same input options. The interface of basins_of_attraction is just perfect.

Final question: So I guess the benefit of all this approach is that you can use a much sparser grid to create this "tracking mode", which doesn't have as heavy dimensional limitations as the normally dense grid. But, I'm not sure if this really saves a lot. Because even in 10 dimensions, 21 points along each dimension are still way too much, yet it is barely dense enough to separate attractors in state space (or maybe it isn't I actually don't know).

I think we can use this functionality to sample basins up to dimension 6 or 7 (I tried it with the Kuramoto oscillator). But as I said in the previous message it will be useful for other purposes.

[Datseris]

the grid for ic's and the grid for identification are now separate things

I've been thinking about it a lot when we were writing the paper. But I have a crucial challenge to this whole idea: What do we gain by having the grid for basins and grid for attractors different? Does it actually help anywhere? The limitation of basins_of_attraction is memory. If either grid is larger, well, we don't really care about the smaller one. The larger grid is what limits everything.

Yeah, I'm not sure I follow your motivation for these grids being different. What I totally understand is having this separate functionality, which maps an initial condition to an attractor. But in the end, it uses 1 grid, no matter what this grid is, as in the same grid both the attractors and their basins are stored. If we have any basin points, we can use them even in random sampling. A trajectory that hits 60 basin cells in a row can still be mapped to the attractor as with anything else. So all in all, you would just have a coarse grid, but I still see no reason for this grid to be separate for basins and for attractors.

---

Basin fraction estimation with random sampling

Yes, this can me used with the bSTAB method, but the bSTAB method itself remains unafected. The design of bSTAB is: "the user provides a function that given an initial condition it returns a vector of features". Well, the users can use this special tracking mode as this function and the vector of features will infact have only one feature, the attractor label.

So this means that we shouldn't need to actually do anything in #219 .

---

the numbering of unidentified cells from 1 to 0.

Besides the fact that SparseArrays cannot be used in the end here, I think it is a good idea to keep this change. 0 is rather special, yet it is unused in our algorithm. 1 on the other hand is "the next free odd number", i.e., the next new attractor label.

[awage]

What do we gain by having the grid for basins and grid for attractors different? Does it actually help anywhere?

Suppose you have attractors in the square [-5, 5]x[-5, 5] and you are interested only in a smaller part of this region you can use this method with a second grid or a random sampling of the smaller region. One grid stores only the attractors while the basins are stored in your custom grid. There is already a mode for this operation in basins_of_attraction but it needs two pass. In fact both functionalities are very similar, the difference is that you have the automatic detection of attractors in one case.

Also, in the attractor grid, the basins are not updated since the initial conditions are different. You don't have the speed improvement of hitting the same basins in a row.

The other benefit that I had in mind was using the sparse matrix for attractors only to avoid memory limitations. But lets forget about this part.

I know that these are not strong reasons to have to separate grids. Anyway this functionality can be useful.

[Datseris]

I'll need a bit of time before I can return here (more than a week probably), so just letting you know in advance. It is good to think about this thoroughly and let it simmer as much as possible before merging anything.

[awage]

All right! I will leave it for now and have a though about parallel computing.

[awage]

I leave here some ideas and results on the parallelization of the code.

• It is possible to get multithreading working with the second mode of operation. I used the tools developed in this PR:

using DynamicalSystems

function multi_threaded_basins(bsn_nfo, integ, grid)
    nthreads = Threads.nthreads()
    bsn = zeros(Int16, map(length, grid))
    bsn_nfo_v = [deepcopy(bsn_nfo) for i in 1:nthreads]
    integ_v = [deepcopy(integ) for i in 1:nthreads]
    I = CartesianIndices(bsn)
    Threads.@threads for ind in I
        u0 = ChaosTools.generate_ic_on_grid(grid, ind)
        bsn[ind] = get_label_ic!(bsn_nfo_v[Threads.threadid()], integ_v[Threads.threadid()], u0)
    end
    return (bsn .- 1) .÷ 2
end

F = 0.27; ω = 0.1;  # smooth boundary
res = 500
ds = Systems.duffing(zeros(2), ω = ω, f = F, d = 0.15, β = -1)
xg = yg = range(-2.2, 2.2, length = res)
grid = (xg, yg)
@time basins, att = basins_of_attraction(grid, ds; T=2*pi/ω, diffeq = (;reltol = 1e-9, alg = Vern9()), show_progress = false)
bsn_nfo, integ = ic_labelling(ds; attractors = att, T=2*pi/ω, diffeq = (;reltol = 1e-9, alg = Vern9()))
@time bsn2 = multi_threaded_basins(bsn_nfo, integ, grid)

Timing for two threads: 30.383229 seconds (2.00 M allocations: 137.854 MiB, 0.16% gc time) 14.088567 seconds (4.52 M allocations: 370.467 MiB, 0.35% gc time)

Distributed and pmap can be also used in the same way.

• I think it is possible to get a parallel version of the "normal" mode of operation. It is almost the same code as before but each thread operates on a part of the complete grid. After the computation the attractors must be matched between threads (the label may be changed from one worker to the other). The results are encouraging but need more work.

Everything is at a early stage still.

[Datseris]

Yes, the first method is certainly parallelizable as there aren't any race conditions. There isn't any shared array that different initia lconditions would write onto.

I don't think it is possible to parallelize the normal version. "each thread operates on a part of the complete grid", but you don't get to decide which part of the grid a trajectory visits, that's something we don't know in advance and hence can't plan. Different initial conditions would write different entries into the basins array if each started first. For me this means the method isn't parallelizable. But I'm interested to see what you have coded so far.

[awage]

Each thread have an independent copy of the grid. It is as dumb as computing n smaller basins. It is a waste of memory, but it is maybe the price to pay to get better in speed.

This code works but I have only tested it on 3 examples:

function multi_threaded_basins(bsn_nfo, integ, grid)
    nthreads = Threads.nthreads()
    bsn = zeros(Int16, map(length, grid))
    bsn_nfo_v = [deepcopy(bsn_nfo) for i in 1:nthreads]
    integ_v = [deepcopy(integ) for i in 1:nthreads]
    I = CartesianIndices(bsn)
    Threads.@threads for ind in I
        if bsn_nfo_v[Threads.threadid()].basins[ind] == 0
            u0 = ChaosTools.generate_ic_on_grid(grid, ind)
            bsn_nfo_v[Threads.threadid()].basins[ind] = get_label_ic!(bsn_nfo_v[Threads.threadid()], integ_v[Threads.threadid()], u0)
        end
    end

    # Remove attractor labels from basins
    for k in 1:nthreads
        ind = iseven.(bsn_nfo_v[k].basins)
        bsn_nfo_v[k].basins[ind] .+= 1
        bsn_nfo_v[k].basins .= (bsn_nfo_v[k].basins .- 1) .÷ 2
    end

    # Match labels between basins.
    for k in 2:nthreads
        match_attractors!(bsn_nfo_v[k-1].basins, bsn_nfo_v[k-1].attractors, bsn_nfo_v[k].basins, bsn_nfo_v[k].attractors)
    end

    # Reconstruct the whole basins
    for k in 1:nthreads
        for ind in I
            if bsn_nfo_v[k].basins[ind] != 0
                bsn[ind] = bsn_nfo_v[k].basins[ind]
            end
        end
    end
    return bsn
end

# Continuous system
ω = 0.1617; F = 0.395 # fractal boundary
res = 200
ds = Systems.duffing(zeros(2), ω = ω, f = F, d = 0.15, β = -1)
xg = yg = range(-3., 3., length = res)
grid = (xg, yg)
bsn_nfo, integ = ic_labelling(ds; grid = grid, T = 2*pi/ω, diffeq = (;reltol = 1e-10, alg = Vern9()))
@time bsn2 = multi_threaded_basins(bsn_nfo, integ, grid)
# 44.148590 seconds (77.60 M allocations: 1.577 GiB, 0.63% gc time, 1.08% compilation time)

#@time basins, att = basins_of_attraction(grid, ds; T = 2*pi/ω, diffeq = (;reltol = 1e-10, alg = Vern9()), show_progress = true)
#  84.830231 seconds (78.47 M allocations: 1.689 GiB, 0.39% gc time, 2.17% compilation time)

The gain is roughly 2 on two threads.

[Datseris]

@awage as a heads up please do not do parallelization additions in the original basins implementation, otherwise this PR will be too complicated to review. Let's please keep this PR to (1) adding the new tracking mode/feature and (2) better naming for the algorithm variables in teh source code. Any parallelization for the original method should be done in a separate PR. Of course, parallelization of the new method feel free to add in this PR.

[awage]

We are on the same page! Parallelization is always a delicate thing and it deserve another PR. I left the comments in this PR just to share some advances in this direction.

[awage]

I finished the modification in this PR. There are many little changes on naming and also the new interface for initial condition labelling.

[Datseris]

I've finished working on #219 finally. That was a big project. When I have free time later I will work here. My conclusion from working on #219, is that whatever method we have here, it should be a different function entirely than what #219 does. In #219 the clustering and approximate features are important. In our method, they aren't. Attractors are identified exactly and hence one does not need to be saving trajectories.

One important question: how is the basin fraction calculation happening in this PR? Do we (1) first find the attractors by letting the machine run for long times, or are the attractors found on the fly while calculating fractions? I.e. a list of attractors is being kept that increases as the simulation runs and some new initial conditions that go to a different attractor are identified. My guess is that the latter would be better for performance?

[awage]

I've finished working on #219 finally. That was a big project. When I have free time later I will work here. My conclusion from working on #219, is that whatever method we have here, it should be a different function entirely than what #219 does. In #219 the clustering and approximate features are important. In our method, they aren't. Attractors are identified exactly and hence one does not need to be saving trajectories.

Exactly, this is why we evolve the trajectory step by step until we decet the attractor or the basin.

One important question: how is the basin fraction calculation happening in this PR? Do we (1) first find the attractors by letting the machine run for long times, or are the attractors found on the fly while calculating fractions? I.e. a list of attractors is being kept that increases as the simulation runs and some new initial conditions that go to a different attractor are identified. My guess is that the latter would be better for performance?

The fractions are not computed as such, but it can be implemented straightforwardly. To answer you question the attractors are detected on the fly using the finite state machine with the normal mode. But the second mode can also be used if the attractors have been supplied.

The function get_label_ic! match the initial condition to an attractor stored in a structure BasinsInfo. All the operations are the same as computing the entire basin using basins_of_attraction but the information on the basins in BasinsInfo is not updated. The reason is that the grid and the initial conditions do not match necessarily.

Here is a summary of the changes in this PR:

• renaming of some variables for clarity of the algorithm.
• The unlabelled initial conditions are now 0 by default instead of 1.
• The second mode of operation that compare distances from state point to attractors have been modified such that the finite state machine is skipped. A side effect is that the performances improved.
• Function ic_labelling to initialize the structure BasinsInfo and the integrator. This is basins_of_attraction under the hood and the two mode of IC labelling are possible (on the fly attractor detection and distance matching).
• Function get_label_ic! that match the IC to an attractor. This function can be used to compute the basins fractions in parallel given a bunch of ICs.
• Some minor changes in estimate_basins!

[Datseris]

Okay, I've had some time here now. I propose a different API which is a bit more intuitive and will make our work easier down the line for integrating with the new features like random sampling.

Define a struct

struct AttractorMapper{B<:BasinInfo, I, A}
    basin_info::B
    integ::I
    attractors::Dict{Int8, A}
end

An instance of AttractorMapper is created using the code of basins_of_attraction with this new keyword ic_lab_mode = true. (notice that basin_info has the grid inside it)

then the labelling interface should be:

mapper = AttractorMapper(.....)
label = mapper(u0) # for a given initial condition it gives you the label. Literally is `get_label_ic!`

Importantly, this new functionality should not be part of basins_of_attraction, especially not via a keyword argument.. Let's keep things simple in the front end API. Of course, the source code of AttractorMapper is like 20 lines or so, as it re-uses everything from basins_of_attraction, but the point is: for a user, this new mode / mapping functionality should just stay separate. And for us, as developers, the cleaner we can make this separation, the better. If we need to copy 5 lines of code from basins_of_attraction to actually initialize the AttractorMapper, then that's okay. We really just need to create the BasinsInfo structure and the integrator. I really want to avoid having multiple functionalities in one function name. It will make our work harder down the line.

In fact, I'd even argue that generating the BasinInfo and the integrator should be a function in itself, that both basins_of_attraction and this new AttractorMapper will call.

@awage if you are unfamiliar with making Function Like Objects (which is what this AttractorMapper is), have a look here: https://docs.julialang.org/en/v1/manual/methods/#Function-like-objects , or just ask me questions.

---

Everything else in the PR is okay!

[Datseris]

Actually @awage let me do the changes I suggested, you've done already a lot of work on this :)

[Datseris]

@awage I've done the changes I outlined. Please do have a look and let me know if you like this code. I'm still not sure about the name AttractorMapper, perhaps we can find a better one :)

Oh, by the way does get_label_ic! return the "correct" labels, 1, 2, 3, 4 or does it return the "only odd number" labels?

[awage]

Nice! Now the two functions are clearly separated. Functions-like objects are a good way to attach a method to type. I will keep this trick in mind.

I changed some small details and fixed the test.

Maybe the last thing I would add is a convenience function to get the attractors from the mapper.

[Datseris]

I've done large changes and everything is broken now, I just don't have the strength to code more for tonight. Will summarize stuff soon. Please don't push any commit to this code until I finish the changes :)

[awage]

I've done large changes and everything is broken now, I just don't have the strength to code more for tonight. Will summarize stuff soon. Please don't push any commit to this code until I finish the changes :)

You are in charge. If you need any help please tell.

[Datseris]

Okay @awage and @KalelR , I've finished the design here. Haven't corrected tests, but doesn't matter for now. Please read the new design and provide feedback. The new design is in the following files: attractor_mapping.jl first. Then, the attractor_mapping_proximity.jl and _recurrences.jl. Finally, the basin_fractions_clustering.jl.

The benefits of this approach:

• Unified api for calculating basins of attraction irrespectively of the method
• Significant elimination of code duplication.
• Simplification of the clustering algorithm docstring
• It all makes more sense from a pure dynamical systems perspective: initial conditions are mapped to attractors with some methodology. The methodology is just an instance of the abstract type AttractorMapper.

If you approve this new design I'll move on to:

1. Change the source code of basins_of_attraction to initialize and use an AttractorsViaProximity mapper and re-use its source code.
2. Fix tests.

Please comment on the names: So far I have all these names AttractorsViaWhatever. I don't like them very much but couldn't think of anything better anyways.

[awage]

Nice! I will be able to start working on Monday on this code. Also I will start the interface for StroboscopicMap.

[Datseris]

Please focus on the StroboscopicMap for now as I'm not done yet here!

[Datseris]

@KalelR just for reference, featurizing the Henon map seems to create problems:

        f(A, t) = [mean(A[:, 1]), mean(A[:, 2])]
        mapper = AttractorsViaFeaturizing(ds, f; Ttr = 100)
        henon_fractions_test(mapper)

errors with

ERROR: ArgumentError: radius NaN must be ≥ 0
Stacktrace:

_dbscan(kdtree::NearestNeighbors.KDTree{SVector{2, Float64}, Euclidean, Float64}, points::Matrix{Float64}, radius::Float64; min_neighbors::Int64, min_cluster_size::Int64) at [C:\Users\m300808.julia\packages\Clustering\tt9vc\src\dbscan.jl]()

I am assuming this is because trajectories diverge to infinity and therefore some feature values become Inf, which makes your optimal radious become NaN, which makes DBSCAN fail. I am not sure how to fix this, i.e., what kind of featurizing functions I should define that do not fail...

---

EDIT: I mean, how do you detect diverging trajectories in the unsupervised version? In the supervised its easy, sure, but what do you do in the unsupervised? I also now noticed that our original tests did not have a diverging to infinity test case. I think the clustering method is just incapable of detecting infinity. It just considers it as another attractor, and hence we should mention this in the docstring.

[Datseris]

Note to self: todo: pretty printing for attractor mappers.

[Datseris]

Note to self: todo: autimatic estimation of ε in Proximity version as the minimum distance between all attractors.

[KalelR]

I am assuming this is because trajectories diverge to infinity and therefore some feature values become Inf, which makes your optimal radious become NaN, which makes DBSCAN fail.

Yes, you're right. I just tested this to make sure. DBSCAN itself just ignores the Inf values, the error was that the optimal radius ε_optimal was NaN.

I am not sure how to fix this, i.e., what kind of featurizing functions I should define that do not fail...

Assigning a new value to Inf features is one way, but it's not ideal I think. If the reassigned value is too large, the optimal radius ε will also become large, because of the way we initially guess viable radii: feat_ranges = maximum(features, dims=2)[:,1] .- minimum(features, dims=2)[:,1] ϵ_grid = range(minimum(feat_ranges)/200, minimum(feat_ranges), length=200)

which might lead to nonoptimal results.

Maybe instead we should instead disconsider all NaN features, and assign them to a same cluster?

I think the clustering method is just incapable of detecting infinity

What do you think is the ideal way to deal with infinities? I would think assigning them to the same attractor makes sense, doesn't it? We can do that.

It just considers it as another attractor, and hence we should mention this in the docstring.

For sure.

[Datseris]

@awage I've pushed the code that computes ε of proximity version automatically as half the minimum distance of points across attractors (if an ε is not given). You can have a look now for why proximity and recurrences give wrong results.

[awage]

I have been thinking about the organization of the code. With the new mapper interface we can call these functions directly from basins_of_attraction. It will also allow us to strip out the proximity search code from the finite state machine part. It will be much cleaner. You just have to loop over the initial condition you need on the grid. The rest is exactly the same. Maybe it was your idea from the beginning.

I made a small diagram (I don't know for you but it is easier for me to visualize the code like that)

[Datseris]

yeap that was my idea from the start, I was just waaiting on the new system types. I will finish the code soon

[Datseris]

So, the bSTAB algorithm is much worse than we thought. It is impossible for it to get any decent mapping of the fractions for the complex Lorenz84 case, where half the time it only finds two clusters, and many times it also attributes -1 to some initial conditions. We should point this out in a next paper. Devise a method/way to quantify the "sensitivity" of each algorithm, like, if we change parameters, how badly does it perform, etc. etc,..

[Datseris]

Damn, this PR is an insane amount of work. I'll need several days to finish it. I have planned out everything to do thankfully, now only coding it remains.

[awage]

It is impossible for it to get any decent mapping of the fractions for the complex Lorenz84 case,

An extension of this method would be to automatically find the feature function that maximizes the distances between attractors in the feature space. Maybe some fancy adaptive filtering method. But this is easier said than done.

If you need some help for tests or coding I am available.

[Datseris]

I am currently re-writting the test suite entirely, because it was simply not maintainable to write tests for each invidual case. Once I am done with that, I Will ask you to add a call to basins_of_attraction in the test suite for all tests I've written and make sure that basins_of_attraction works as before, with all possibliities as before. It is crucial that this PR does not introduce breaking changes, only deprecations. Code that worked before should still work now... I'll let you know once I am done. Almost there......

[Datseris]

@awage I have written the tests in a generalizable way. I am very happy with the way I've written the tests, but there is still much work to do and I would appreciate any help because I need to turn to and release Agents.jl 5.0. The tests are in the file attractor_mapping_tests.jl. At the end of the file I outline the TODOs remaining:

TODO: Add a call to `basins_of_attraction` within the `test_basin_fractions` function.
TODO: Tests for  projected system (magnetic) and poincare map (thomas cyclical)
TODO: Tests of Lorenz84 fail, significant differences between recurrence and proximity
TODO: Tests of Duffing fail, significant differences between recurrence and proximity

anything you can hack is great help, I'll continue working on this when I get some more time!

[awage]

OK! I'll start to work on this.

[awage]

The unsupervised method of clustering features gives me a lot of trouble. Now I have a strange error with the Lorenz84 system. I leave this part for now and focus on the Poincaré map tests.

[awage]

The PoincareMap integrator code is becoming ugly. There is an ambiguous choice about the state of the integrator: is it the full state of the dynamical system or just the state on the plane. There is surely a simpler approach but right now I am too tired to think straight.

[awage]

I think I found the solution. Maybe it is too complicated: I use a projected integrator to wrap a Poincaré Map.

This way there are minimal changes to the Poincaré Map integrator interface. It seems to work but I don't know if there is a regression because of this integrator within an integrator.

[awage]

The interface is not intuitive but it works:

    pmap = poincaremap(ds, (3, 0.0), 1e6;
        rootkw = (xrtol = 1e-8, atol = 1e-8), diffeq=(reltol=1e-9,)
    )
    psys = projected_integrator(pmap, 1:2, [0.])

I think it is possible to create a better integration of both integrators together. Maybe with a keyword in the Poincaré Map interface such as projected = true so that the integrator is defined on the plane?

[Datseris]

Unfortunately I do not agree with combining or connecting the projected integrator with the Poincare map. They are fundamentally different things. The Poincare map is not projecting the state into a lower dimensional plane, but rather it intersects it with it. This is why the Poincare map has the fundamental property of 1-to-1 correspondance with the continuous dynamics. A projected system does not have this property. Hence, I think it is more intuitive if these two things are separated.

As for the code, one possible solution is a simple low level function that for all integrators returns get_state while for the poincare map it returns the "intersected" state. This changes nothing in no integrators and is a 5 lines of code change.

[Datseris]

@awage can you please revert the changes you did to the Poincare map, and starting from there we see what we is missing and what we need to do in code? It is really difficult for me to do this via the commit history, because there was also a lot of testing around. Generally speaking, if you have an idea to do anything breaking, it is always better to first discuss it, and then commit code after we have settled on a design. This will likely save time of re-thinking or re-implementing something.

[awage]

Unfortunately I do not agree with combining or connecting the projected integrator with the Poincare map. They are fundamentally different things. The Poincare map is not projecting the state into a lower dimensional plane, but rather it intersects it with it. This is why the Poincare map has the fundamental property of 1-to-1 correspondance with the continuous dynamics. A projected system does not have this property. Hence, I think it is more intuitive if these two things are separated.

As for the code, one possible solution is a simple low level function that for all integrators returns get_state while for the poincare map it returns the "intersected" state. This changes nothing in no integrators and is a 5 lines of code change.

Ok, I understand your argument. Although I still find the idea of reducing the dimension appealing. I though I was on the right track this is why I pushed the commits. I will revert to where you left it.

[Datseris]

Although I still find the idea of reducing the dimension appealing.

But I agree with that. I believe that we can have the best of both worlds: both allow a dimension reduction (within basins), but also not merge the projected system and the Poincare map. Notice that it is not possible to have a dimension reduction in the case of a generic Poincare map where the plane has arbitrary orientation. Only in the case where the plane is orthogonal to one of the dimensions of the system.

(on the other hand you might just say that the arbitrary plane orientation anyways doesn't work with the recurrences code, which is true.,.,.)

[Datseris]

I've pushed some commits I was working on a different branch (independent improvements).

Alright, lets reorganize ourselves. So, at the moment, the "problem" we have with the Poincare map is that the grid is D-1 dimensions, but the state of the map is D dimensional, right? My proposal is to just define a new function _get_state_on_grid, which for all integrators is just get_state, but for the poincare map it returns the state on the plane, which has one less dimensio0n. This internal function _get_state_on_grid is just used in the AttractorsViaRecurrences code.

This has the advantange of not breaking anything. What do you think of this approach @awage ?

[Datseris]

Just letting you know that I've tested this locally @awage . I changed absolutely nothing, and only added one more reinit! method to the PoincareMap tht works for a given state with D-1 dimensions. _Everything works, basinfractions, etc., even though I had to write no extra code.

[awage]

Just letting you know that I've tested this locally @awage . I changed absolutely nothing, and only added one more reinit! method to the PoincareMap tht works for a given state with D-1 dimensions. _Everything works, basinfractions, etc., even though I had to write no extra code.

That's wizardry and julia black magic. Thanks! I will try to have a look this weekend.

[Datseris]

That's wizardry and julia black magic. Thanks! I will try to have a look this weekend.

Wait until I push the commits first. Doing some small debugging at the moment. Will probably not have time to finish this today because i'm in london now, but will do it whenever I have some spare time. Soon soon.

[awage]

AttractorsViaRecurrences

We will need a function of this kind. Another possibility is to get the index of the variables that correspond to the grid (D-1 indices in the case of the Poincaré Map). The only place where these indices are needed is in the function basin_cell_index where you get the correspondence between the coordinates on the grid and the index of the basins matrix.

These indices can be stored in the BasinsInfo struct. What do you think?