Massive reform of all basins-related code

[Datseris]

cc @awage and @KalelR . Also see the latest status quo of code in PR #237

so I've been pondering over our basins code, and many things are suboptimal. Not in terms of performance, but in terms of code clarity, things being where they do not really belong, and code duplication. So I will outline here a massive code reform that is also inspired by the new AttractorMapper functionality, which will itself be given to the basin_fractions function.

Once I'm finished working on PR #237 things will become even more clear, but I'm writing my notes here anyways.

Definition of Proximity Version. By "proximity version" I refer to the algorithm that maps initial conditions to attractors based on the distance of the current trajectory state from the attractors. The first attractor that is < ε from the current state gets mapped as the attractor of the initial condition. This is called the "supervised method" so far in our code.

---

• [ ] The functionality of the basins_of_attraction where "any dynamical system, including stroboscopic or poincare maps" can be passed is nice and apparently also useful for the proxity version or other future code aspects. "Handling any kind of dynamical system" should be detached from basins_of_attraction.
• [ ] For this, we would have to create a StroboscopicMap struct, similarly with the existing PoincareMap. That would clarify the code quite a bit and not make us do all this special handling for stroboscopic maps in the basins function. Instead, we would be able to handle it exactly like a discrete system. One needs to make the struct and define integrator, reinit!, step!, get_state methods for it. That's it. This should be done in DynamicalSystemsBase.jl.
• [ ] Similarly, we can define a "ProjectedSystem", which is a dynamical system whose state and dimension are taken under projection. The complete_state and all of this kind of business in basins_of_attraction should really be put into the "ProjectedSystem". Then integrator(::ProjectedSystem; diffeq) returns a wrapper around a standard ODE integrator, that also has the complete state and project state functions. Then, step!, reinit! and get_state do what they are supposed to do via simple wrapping. If we do this , then in all of the basins code, there would be NO CODE that cares about what the actual dynamical system is. ProjectedSystem should also go into DynamicalSystemsBase.jl and needs the same methods as StroboscopicMap.
• [ ] Proximity version really has nothing to do with a grid. Its code needs to be detached completely from the current basins_of_attraction. It should be fully implemented by itself in the attractor mapping by proximity functionality. The basins_of_attraction can then use the clean code of the proximity version if a user has provided attractors. I am currently working on this and almost done.
• [ ] The "supervised" functionality of bSTAB seems to be worse/inferior to the proximity version. What are its benefits, if any? It still needs a featurizer, and it still does clustering. @KalelR can you comment on occasions where the supervised method of bSTAB would be preferred over what I call here "proximity version"?

---

As you can see, this is a lot of work, but it will be a huge advantage to the library to do these changes. It will make everything simpler, and easier to extend. Also, it will create transferable functionality. E.g., at the moment this nice handling of projected systems can only be used within the basins, but it may be useful in other algorithms as well. Any help is welcomed!

[KalelR]

The first attractor that is < ε from the current state gets mapped as the attractor of the initial condition.

Is it the first attractor, or the one closest to the current state?

The supervised method in basin_fractions follows directly from the bSTAB paper, where they use a first nearest neighbors algorithm ( kNN, with k=1). So it identifies the attractor ("template" in the code) which is closest to the current state ("features"). That would be the difference to "proximity version", if the latter finds the first attractor, not the closest. And its an important difference, finding the closest attractor makes the most sense Id say. If proximity version finds the closest attractor, then there should be no difference. The only practical scenario I can think for this supervised scenario is the obvious, when the user knows the features defining the attractors, and possibly wants to be sure these are ones the algorithm finds/uses.

Otherwise, I really have no clue if this 1NN is less performant than the "proximity version". What I can say from the tests I made is that it is really quick. Much quicker than the unsupervised method, if there are lots of ics. It also allows the user to pass other distance metrics and a maximum distance threshold. I guess "proximity version" could implement these without much difficulty, so I don't see any real advantage to the proximity version in this performance sense, except kNN is already implemented.

[awage]

Is it the first attractor, or the one closest to the current state?

This the one closest to the current state.

If proximity version finds the closest attractor, then there should be no difference. The only practical scenario I can think for this supervised scenario is the obvious, when the user knows the features defining the attractors, and possibly wants to be sure these are ones the algorithm finds/uses.

The difference is that the "proximity version" checks the distances at each time-step and stops when the distance falls bellow ε. When using the feature, you first integrate the full trajectory, compute the feature vector and then match it to the template. I think this approach can be more robust in some cases (ill defined attractors or two attractors very nearby points in the state space).

[awage]

As you can see, this is a lot of work.

You love titanic tasks!

Any help is welcomed!

Sure! Happy to contribute. Tell me what you need.

[Datseris]

That would be the difference to "proximity version", if the latter finds the first attractor, not the closest.

No. I mean, yes, technically. It maps a trajectory to the first attractor whose distance to the trajectory is < ε. But this happens DURING the time evolution of the orbit. So it is almost certainly guaratneed that the first attractor with distance < ε is also the closest. Alex described it well. Notice that by default ε is really small, 1e-3.

As far as I can tell, the proximity version is superior to the "supervised featurizing version" for the following reasons:

• It does not need to store in memory a full trajectory to produce features.
• It knows when to terminate the time evolution, which, as far as I can tell, is the biggest performance problem.
• It actually computes distances in state space, not feature space, and thus is more rigorous
• The nearest neighbor in real state space is simply more likely to be the "correct attractor", than the nearest neighbor in feature space. Why? Every attractor is guaranteed to be separated in state space. In feature space, this is only guaranteed if the user chooses appropriate features.

---

except kNN is already implemented.

Yes. this is a good point. Actually, I can imagine a scenario where one would prefer the features: for extremely high-dimensional dynamical systems, e.g. 100s or 1000s of dimensions. There our approach, which uses KDTrees, would become visibly slow, because KDTrees do not perform so well in such high dimensions. However, the featurizing would typically project the system to 1-10 dimensional space. So yes, I agree, we will keep this version.

---

You love titanic tasks!

Absolutely! haha, but that's how you make great things! Similarly to re-writing a paper 3 times until you get it right. Fun fact: I have re-written DynamicalSystems.jl twice in its early stages. Looking back, it was a great choice!!!

---

Sure! Happy to contribute. Tell me what you need.

The first two steps is to make (1) the ProjectedSystem structure and (2) the StroboscopicMap structure. The first should be a wrapper that contains the dynamical system and the complete_state business and then defines an integrator method that makes a new wrapper around new_integ = integrator(projected_system.ds) so that get_state(new_integ) = projection(get_state(new_integ.integ)). StroboscopicMap is pretty much the same as the PoincareMap.

[Datseris]

@awage some quick comments:

1. Please add the stroboscopic map code in DynamicalSystemsBase.jl
2. Make a dedicated file for it
3. Do not worry at the moment for making Stroboscopic map work with the basins code. The stroboscopic map is a feature on its own that the basins must utilize. hence, the basins needs to be adapted after teh stroboscopic code is merged in DynamicalSystemsBase.jl.

[KalelR]

I see, thanks for the explanations, guys! I have a question: for systems of oscillators with multistability there are several attractors in a bounded state space. Would "proximity version" work well in finding the correct attractor?

As far as I can tell, the proximity version is superior to the "supervised featurizing version" for the following reasons:

Yes, those are very good points.

Actually, I can imagine a scenario where one would prefer the features: for extremely high-dimensional dynamical systems, e.g. 100s or 1000s of dimensions.

Ah, excellent point. In systems of oscillators for instance some good features are very clear, like the degree of synchronization. But their values at the attractors are generally unknown, so the supervised version would not be better anyway. As far as I understand, Euclidian distance in high dimensional spaces can be problematic, which could be an issue besides just the perfomance of KDTree. So supervised version can come in handy, though I guess not as often as proximity version. Since it is already implemented, I see no reason to remove it.

[Datseris]

This hjas been implemented in the ongoing #237