Basin fraction continuation API + recurrences method

[Datseris]

This PR

1. Defines the API for basins fractions continuation algorithms
2. Implements the recurrences-based method for it
3. Adds some tests for the method, but more need to be added later.

[Datseris]

@awage I need your help here. Please try to run the file test\basins\basins_continuation_tests.jl. It will print some info. Some of the info is:

[ Info: Starting basins fraction continuation.
[ Info: p = [1.0, 1.0, 1.0]
[ Info: AttractorsViaRecurrences found new attractor with id: 1
[ Info: AttractorsViaRecurrences found new attractor with id: 2
[ Info: AttractorsViaRecurrences found new attractor with id: 3
p = [1.0, 1.0, 0.9]
[ Info: AttractorsViaRecurrences found new attractor with id: 4
[ Info: AttractorsViaRecurrences found new attractor with id: 5
[ Info: AttractorsViaRecurrences found new attractor with id: 6
rmap = Dict{Int16, Int16}(5 => 3, 4 => 2, 6 => 1)
fs = Dict(2 => 0.39, 3 => 0.31, 1 => 0.3)
p = [1.0, 1.0, 0.8]
[ Info: AttractorsViaRecurrences found new attractor with id: 7
[ Info: AttractorsViaRecurrences found new attractor with id: 8
[ Info: AttractorsViaRecurrences found new attractor with id: 9

At each parameter I am trying to reset the internal AttractorsViaRecurrences code. But it doesn't work. The new attractors at the new parameters are oddly starting their enumeration from the previous attractors, as if my reset failed completely. Can you please help? what am I missing? This is how I reset:

function reset!(mapper::AttractorsViaRecurrences)
    empty!(mapper.bsn_nfo.attractors)
    mapper.bsn_nfo.basins .= 0
    # TODO: Why doesn't this actually set the attractor starting labels to 1...???
    mapper.bsn_nfo.prev_label = 0
    return
end

The code at the moment also errors, but I want to figure out first how to do the resetting properly.

[awage]

To algorithm starts with these parameters:

    state = :att_search
    current_att_label = 2 
    visited_cell = 4
    consecutive_match = 0 
    consecutive_lost = 0 
    prev_label = 0 

[Datseris]

@awage it works now however I need your help a bit more. I think, and you can tell me what you think as well, that we should NOT reset the "current attractor label" to 1. So that if we have e.g., 3 attractors, but then we have two, but then three again (with the two always the same), the new third attractor should get the label "4". However, I cannot succeed in this. Have a look at the current test file test\basins\basins_continuation_tests.jl. It's the magnetic pendulum and towards the end (the keys are printed) the third attractor goes in and out of existence because it has really small basin and our random sampling doesn't find it. Hower the keys are always 1, 2, 3 and I never see a key "4". Why...?

[awage]

I will have a try. Do you know a value of the parameters where this basin appears?

[Datseris]

Just run the test file it is set up to print info!

[awage]

Tell me if I am wrong:

• It starts with a normal basins_fractions estimation with the first parameter in the list.
• The parameter changes and previous attractors are used as seeds to match them against the current ones.
• The program launch a set of random ics to sample the state space and store the new attractors if any.

I have launched the test file with different tolerances and number of random ics and I have always the same result:

k = [1, 2, 3]
k = [1, 2, 3]
.
. 
.
k = [1, 2, 3]
k = [2, 3]
k = [2, 3]
.
.
.
k = [2, 3]

It seems that the attractor 1 disapears at the end. I have also checked the position of the attractors and they are coherent values. From what I see here the algorithm does its job perfectly for this example.

Do you have any trace of the behavior you mentioned?

[Datseris]

Yes, I didn't get my point correctly across. What you will notice somewhere in the middle of the parameter space is that this happens:

k = [1, 2, 3]
k = [2, 3]
k = [1, 2, 3]
k = [2, 3]
...

The attractor "1" (or whatever number it is) appears and reapperas just before it completely vanishes. That's because its fraction is so small that we may, or may not, randomly sample an initial condition leading to it.

However, I argue that the correct behavior is that it should do:

k = [1, 2, 3]
k = [2, 3]
k = [2, 3, 4]
k = [2, 3]
...

I.e., the appearance of a new attractor should get a HIGHER index, 4, not a LOWER index 1. I do not know why the code doesn't do this, even though I've set it up to do it.

[awage]

I have the answer: the basins fraction doesn't represent some of the basins correctly. When I looked at the list of attractors there is sometimes a mismatch between the number of attractors and the number of keys in the the basins fractions.

First you have the keys of the attractors and then the keys of the basin fractions. You will see that sometimes one of the basin is not sampled at all:

collect(keys(k)) = Int16[2, 3, 1]
collect(keys(k)) = Int16[2, 3, 1]
collect(keys(k)) = Int16[2, 3, 1]
collect(keys(k)) = Int16[2, 3, 1]
collect(keys(k)) = Int16[2, 3, 1]
collect(keys(k)) = Int16[2, 3, 1]
collect(keys(k)) = Int16[2, 3, 1]
collect(keys(k)) = Int16[2, 3, 1]
collect(keys(k)) = Int16[2, 3, 1]
collect(keys(k)) = Int16[2, 3, 1]
collect(keys(k)) = Int16[2, 3, 1]
collect(keys(k)) = Int16[2, 3, 1]
collect(keys(k)) = Int16[2, 3, 1]
collect(keys(k)) = Int16[2, 3, 1]
collect(keys(k)) = Int16[2, 3, 1]
collect(keys(k)) = Int16[2, 1]
k = [1, 2, 3]
k = [1, 2, 3]
k = [1, 2, 3]
k = [1, 2, 3]
k = [1, 2]
k = [1, 2, 3]
k = [1, 2]
k = [1, 2, 3]
k = [1, 2]
k = [1, 2, 3]
k = [1, 2]
k = [1, 2]
k = [1, 2, 3]
k = [1, 2]
k = [1, 2]
k = [1, 2]

I think the algorithm works correctly. A way to test the numbering it is to shuffle the parameters so that sometimes you have two basins and sometimes three.

[awage]

I have tried with a list of random parameters and the algorithm does not work as intended. I think the problem is in the matching of the attractors not in the recurrence algorithm.

[awage]

Have a look at this sequence:

[ Info: AttractorsViaRecurrences found new attractor with id: 4
[ Info: AttractorsViaRecurrences found new attractor with id: 5
[ Info: AttractorsViaRecurrences found new attractor with id: 6
rmap = match_attractor_ids!(current_attractors, prev_attractors; metric, threshold) = Dict{Int16, Int16}(5 => 3, 4 => 2, 6 => 1)
fs = Dict(2 => 0.494, 3 => 0.017, 1 => 0.489)
p = [1.0, 1.0, 0.14652488546165196]
[ Info: AttractorsViaRecurrences found new attractor with id: 7
[ Info: AttractorsViaRecurrences found new attractor with id: 8
rmap = match_attractor_ids!(current_attractors, prev_attractors; metric, threshold) = Dict{Int16, Int16}(7 => 2, 8 => 1)
fs = Dict(2 => 0.495, 1 => 0.505)
p = [1.0, 1.0, 0.27967666953502873]
[ Info: AttractorsViaRecurrences found new attractor with id: 9
[ Info: AttractorsViaRecurrences found new attractor with id: 10
[ Info: AttractorsViaRecurrences found new attractor with id: 11
rmap = match_attractor_ids!(current_attractors, prev_attractors; metric, threshold) = Dict{Int16, Int16}(11 => 1, 10 => 1, 9 => 2)

The last operation matches the two attractors 10 and 11 to the previous attractor 1. Do you know why?

[Datseris]

I have no idea where this code is coming from, nor what is threshold :D I'll need more context before I can say anything.

[awage]

Ok, I think the problem is this threshold. By default it is Inf. If I set:

continuation = RecurrencesSeedingContinuation(mapper; threshold = 1.)

Then the result is:

collect(keys(k)) = Int16[2, 3, 1]
collect(keys(k)) = Int16[2, 3, 1]
collect(keys(k)) = Int16[2, 1]
collect(keys(k)) = Int16[2, 12, 1]
collect(keys(k)) = Int16[2, 12, 1]
collect(keys(k)) = Int16[2, 12, 1]
collect(keys(k)) = Int16[2, 12, 1]
collect(keys(k)) = Int16[2, 12, 1]
collect(keys(k)) = Int16[2, 1]
collect(keys(k)) = Int16[29, 2, 1]
k = [1, 2, 3]
k = [1, 2, 3]
k = [1, 2]
k = [1, 2, 12]
k = [1, 2, 12]
k = [1, 2, 12]
k = [1, 2, 12]
k = [1, 2]
k = [1, 2]
k = [1, 2, 29]

The numbering is a bit odd but it is the intended behavior!!

[Datseris]

but what code gives you this??? The code that currently exists in the pr never gives an ID higher than 3, and that is exactly my problem.

[awage]

I pushed a commit so you can view where the code changes.

[Datseris]

but it doesn't make any sense to use random parameters. attractor matching works by distance, that's why we are calling the method "continuation". It "continues" the attractors in the previous parameter to the next one that is very close. if the parameters are randomized instead of sorted then there is no real way to match the attractors.

Can we have a call to talk about this efficiently?

[awage]

I cannot call you today. Tomorrow morning maybe.

Forget about the random parameters, it is just a way to make attractors appearing and disapearing artificially.

[awage]

I have a limited internet connection these days, I cannot call you. I pushed a commit with a discrete system. You will see that the algorithm works but the numbering must be changed.

Also there is a problem when a discrete integrator is used with set_parameter!(mapper.integ, ....). The program throws an error.

[Datseris]

You will see that the algorithm works but the numbering must be changed.

But I don't understand what this means. Which algorithm works? There are three algorithms here:

1. Find attractors using recurrences, which is tested extensively
2. Match attractors with match_attractor_ids! function, which is also tested extensively
3. Continuate basin fractions, which we try to test here

Which of these work, which do not (if any), and what do you mean that the "numbering must be changed"? I see when I run the henon example that sometimes we get one and sometimes two attractors. Is this a problem of our already published paper on recurrences, that our algorithm (1) finds false attractors? Or is this okay? I also see that in the last 5 parameters for henon the fractions and the attractors have different keys. This is definitely something wrong, but I don't know yet which of the (1, 2, 3) algorithms is the cause.

One way this may happen is because the seeds go to attractors, but the random sampling doesn't sample one of the two. Therefore, I must ensure that the seeded trajectories are also counted in the fractions of the attractors. The code doesn't do that at the moment, but it is a large change.

Also there is a problem when a discrete integrator is used with set_parameter!(mapper.integ, ....) . The program throws an error.

I've fixed that in DynamicalSystemsBase 2.7.4. BTW if you have seen the error already please paste the error message instead of saying that it throws the error. This way I can just go and fix the error, instead of modifying the code to the state you had it before, running it, seeing the error, and then fixing it.

[awage]

I think all the algorithms work correctly. I will try to check more thoroughly the Hénon example (I cannot tell when I will do it).

The problem with the numbering can be solved:

• The numbering of the attractors in the mapper is not set to 1 after each change of parameter (this is the current_att_label variable).
• When an attractor cannot be matched against a previous attractor in the list, the algorithm assigns the current attractor number of the mapper (this number can be quite high).

Since the numbers are arbitrary, we can rename the keys sequentially at the end of the routine.

One way this may happen is because the seeds go to attractors, but the random sampling doesn't sample one of the two. Therefore, I must ensure that the seeded trajectories are also counted in the fractions of the attractors. The code doesn't do that at the moment, but it is a large change.

This is more difficult to solve. Maybe it is sufficient to add a key with the fraction equal to zero if this attractor exists but the basin is not sampled.

I've fixed that in DynamicalSystemsBase 2.7.4. BTW if you have seen the error already please paste the error message instead of saying that it throws the error. This way I can just go and fix the error, instead of modifying the code to the state you had it before, running it, seeing the error, and then fixing it.

Got it! Many thanks!

One thing that should be tested are other metrics for the matching of attractors (features, Hausdorff distance ...)

[awage]

I have checked the continuation with the Henon system and it is working. I don't know how to set up a reliable test for this function. What should be tested? The keys of the attractors? The basins fractions?

[Datseris]

Since the numbers are arbitrary, we can rename the keys sequentially at the end of the routine.

Yes, that's the way to go.

This is more difficult to solve

Yes, but I know how to solve it. Doing it is all the remains now!

I don't know how to set up a reliable test for this function.

Don't worry, I'll do this.

We need a test for an analytically resolved system, like a standard system where we know it has a bifurcation from e.g., a fixed point to a limit cycle. In this test I'll test a different "similarity function" that maps a limit cycle to itself and a fixed point to itself.

[Datseris]

I'm working on this now, please open a parallel branch if you are also working on it!

[Datseris]

Now we are "renumbering" an attractor 232 or so into attractor 1. But attractor 1 existed in the past. So this means, we have "matched" these two attractors even though they may have nothing to do with each other. Only in the magnetic penduylum case they do happen to have something to do with each other.

So this is currently wrong. The renumbering must use the next positive integer higher than existing keys,, in this case 4.

[Datseris]

I've fixed the missmatch between attractors and fractions having different keys.

[Datseris]

I think we need to change the key used from Int16 to at least Int32 or even just Int. In the magnetic pendyulum we reach attractor ID 284, and typemax(Int16) is only 32,000 or so. Maybe a system with a lot of attractors, or very fine parameter sampling, will reach this limit very easily.

[Datseris]

There is one more bug: even if show_progress = false the algorithm always prints the first attractors found. I do not understand why.

[awage]

I will not work on this until tonight. you can work on this branch.

[awage]

Is there anything I can do at this stage? I can work a bit tonight and a little bit next week.

[Datseris]

The code works as expected as far as I can tell. You can come up with some good tests! You don't have to add anything in the PR. If you come up with a good test just post it here.

[maximilian-gelbrecht]

@Datseris are you planning to merge this soon? then I'd wait with starting to work on a MCBB implementation until this is merged, or does it still take some work?

[Datseris]

@maximilian-gelbrecht you shouldn't wait for me to merge this. At the moment we will need some time to come up with good tests. However, this shouldn't really hold you back. What you need is:

1. the API defined in src\basins\continuation\continuation_basins_fractions_api.jl
2. the way we currently find attractors via clustering without varying a parameter, in src\basins\mapping\attractor_mapping_featurizing.jl
3. The clustering code in the clustering subfolder. Actually, for this, I think you don't have to care.

What I recommend is for you to start a branch off- of this PR's branch. Then, you can open a new PR into this branch. In the description of the PR you can describe the way you think of making the code, so we can give you feedback before you write it, in case we notice something is being more complicated than necessary. Does this sound good? this way you can also work at your own pace without being dependent on this PR. Once this PR is merged your PR gets rebased into main anyways.

[maximilian-gelbrecht]

Yeah, sure, when it's mainly just tests missing, then I can also base the fork / PR on this branch. I'll do that.

[Datseris]

Side remark to myself: I see issues with continuity in this approach where only the exactl previous parameter slice is used. E.g., just consider the traditional logistic map bifurcation diagram. Periodic windows each get their own attractor, but chaotic attractor in the different windows should be the same. How to ensure that...? I am thikning of a keyword argument "previous_slices::Int" that keeps this amount of slices prior and compares to that...? We'll see. this is a point to worry about in the future, after merging this PR.

[Datseris]

Holy shit isn't that cool. I've set it up so that the "distance" in "matching attractors via distance" can also be any arbitrary user-defined function. In the henon map example I use as "distance" the difference of the base-2-logarithm of teh length of each attractor. With threshold = 1, This basically matches attractors whose "period" is within a factor of two. So the period 7 attractor and the chaotic attractor get properly matched into different attractors. HAR HAR HAR

[Datseris]

@awage you should put here your kuramoto test file

[Datseris]

Alright, I am merging this in. The interface is finished as far as I can tell and is extendable. So, I'm merging this in to keep things clean, and I'll open a new issue for the remaining todos.

[awage]

Alright, I am merging this in. The interface is finished as far as I can tell and is extendable. So, I'm merging this in to keep things clean, and I'll open a new issue for the remaining todos.

I have been away for a while. I will to update the test file and understand your changes.