Match attractors by BasinEnclosure

[KalelR]

Idea

So far in the (global) continuation methods, the attractors are matched by some arbitrarily-defined distance function between them. This is quite nice, but isn't necessarily the best method in all scenarios. I've been working on some highly multistable systems in which the attractors can be quite close in state space and also feature space. Tracking the attractors properly with some distance function has proved very difficult.
Instead, also following literature on this, I decided it makes more sense to track the attractors differently. For an attractor A at parameter p1 and attractor B at p2>p1, I would like the algorithm to match A and B if a trajectory starting on the endpoint of A converges to B when p = p2. In this sense, to consider match attractors when the basin of B includes A, or at least part of A.

Implementation idea

I have implemented one way of doing this, but I'm pretty sure it is not the most efficient. The code is also not the very well written, but due to time constraints I cannot work on improving this now. In any case, because of PR #132 I will share my code here.

The basic idea is to use the following matching function:

function distance_matching_endpoint(A,B) #After and Before
    dist =  evaluate(Euclidean(), A[1], B[end])
    return dist
end

which will match attractors if the endpoint of the previous attractor is the same as the starting point of the subsequent attractor. Then, we need to implement a function that will achieve this, by:

1. Running the continuation as normal
2. For each attractor found in the continuation, go to the first parameter p1 in which it was found; then, use its endpoint as the initial condition for a trajectory integrated at the subsequent parameter p2 - this generates what I called att_future_integrated
3. Find, in p2, the attractor that corresponds to this trajectory. To do this, match att_future_integrated with the already found attractors in p2 (in the code, I called them atts_future), using some distance function. If the attractor at p1 continues to exist at p2, this att_future_integrated should be very similar (basically the same) as one of the already-found attractors at p2 - they should be the same attractor. So this distance should basically be zero. If the attractor at p1 ceases to exist, the trajectory will converge to another attractor. Then there are two possibilities. First, the attractor to which it converged exists at p1. So here two attractors at p1 converge to the same attractor at p2, I call them coflowing attractors. Second, the attractor to which it converged emerges at p2. Here, the attractor at p1 ceases to exist, and another attractor emerges at p2.
4. If there are coflowing attractors, I consider that matching should be done between the attractor at p1 that is closest to the the attractor at p2. This is implemented in _find_coflowing and _key_legitimate.
5. After deciding the corresponding attractor at p2, replace it by the trajectory starting at the endpoint of the attractor at p1.

I'm sorry if this is too confusing! I myself was confused writing it, so I guess the explanation and code reflect that. I'll reiterate the basic idea, though: keep the continuation and matching function as they are. Only add a new function that essentially replaces the attractors with a trajectory that starts at the endpoint of the attractors at the previous parameter. Then, simply apply the matching function using distance_matching_endpoint.

Code

function replace_by_integrated_atts(ds, featurizer,atts_all, prange, pidx; T, Ttr, Δt, coflowing_threshold=0.1)
    all_keys = unique_keys(atts_all)
    atts_new = deepcopy(atts_all)
    for idx in 1:length(prange)-1
        p_future = prange[idx+1]
        ds_copy = deepcopy(ds)
        set_parameter!(ds_copy, pidx, p_future)
        atts_current = atts_new[idx]
        atts_future =  atts_new[idx+1]
        atts_future_integrated  = Dict(k=>trajectory(ds_copy, T, att_current[end]; Ttr, Δt)[1] for (k, att_current) in atts_current) 
        ts = Ttr:Δt:T
        keys_ilegitimate_all = Int64[]
        
        feats_current = Dict(k=>featurizer(att, ts) for (k, att) in atts_current)
        feats_future = Dict(k=>featurizer(att, ts) for (k, att) in atts_future)
        feats_future_int = Dict(k=>featurizer(att, ts) for (k, att) in atts_future_integrated)

        for (k_future_int, att_future_int) in atts_future_integrated
            if k_future_int ∈ keys_ilegitimate_all 
                continue
            end
            key_matching_att = _find_matching_att(featurizer, ts, att_future_int, atts_future)
            
            #before replacing, must check if att_future_int is legitimate (and doesn't come from an att that disappeared!)
            keys_coflowing = _find_coflowing(featurizer, ts, att_future_int, k_future_int, atts_future_integrated, coflowing_threshold)
            if length(keys_coflowing) == 1 #no coflowing, att is legitimate 
                att_replace = att_future_int 
            else #coflowing 
                key_legitimate = _key_legitimate(featurizer, ts, att_future_int, atts_current)
                if key_legitimate == k_future_int #check if this is legitimate, i.e. if it is close to a previousy existing att 
                    att_replace = att_future_int 
                    keys_ilegitimate = filter(x->x==key_legitimate, keys_coflowing)
                    push!(keys_ilegitimate_all, keys_ilegitimate...)
                else  #"$k_future_int is coflowing but not legitimate"
                    continue
                end
            end
            
            atts_new[idx+1][key_matching_att] = att_replace 
        end
    end
    return atts_new 
end

function keys_minimum(d)
    vals = values(d)
    ks = keys(d)
    min = minimum(vals)
    idxs_min = findall(x->x==min, collect(vals))
    keys_min = collect(ks)[idxs_min] 
end

function _find_matching_att(featurizer, ts, att_future_int, atts_future) 
    f1 = featurizer(att_future_int, ts) 
    dists_att_to_future = Dict(k=>evaluate(Euclidean(), f1, featurizer(att_future, ts)) for (k, att_future) in atts_future) 
    label_nearest_atts = keys_minimum(dists_att_to_future) 
    if isempty(label_nearest_atts)
        @error "No matching attractor. Shouldn't happen!"
    elseif length(label_nearest_atts) == 1 
        key_matching_att = label_nearest_atts[1] 
    else #coflowing attractors 
        @warn "Coflowing attractors. This will be dealt with later."
    end
    return key_matching_att
end

function _find_coflowing(featurizer, ts, att_current, k_current, atts, coflowing_threshold) #coflowing if dist(att, atts) <= coflowing_th
    f1 = featurizer(att_current, ts)
    dist_to_atts = Dict(k=>evaluate(Euclidean(), f1, featurizer(att, ts)) for (k, att) in atts) #includes it to itself
    idxs_coflowing = findall(x->x<=coflowing_threshold, collect(values(dist_to_atts)))
    keys_coflowing = collect(keys(dist_to_atts))[idxs_coflowing]
    return keys_coflowing
end

function _key_legitimate(featurizer, ts, att_current, atts_prev)
    f1 = featurizer(att_current, ts)
    dist_to_prev = Dict(k=>evaluate(Euclidean(), f1, featurizer(att_prev, ts)) for (k, att_prev) in atts_prev)
    ks_legitimate = keys_minimum(dist_to_prev)
    # @info "finding legitimate, $(dist_to_prev), $ks_legitimate"
    if length(ks_legitimate) == 1 
        return ks_legitimate[1]
    else 
        @warn "Two legitimates. Deal with it."
    end 
end

[Datseris]

Continuing the discussion here to stay focused, from #132 (comment)

it has the benefit of just adding a new function, and not changing the continuation or matching functions.

Hm, but is it really a benefit to add a new function? I am not sure. The more functions you add, the more complicated the API of a software is. At the moment we already have in place the concept of a matcher, and the function that matches (replacement_map!). If you add a new function you increase complexity by adding both a new algorithm that matches, as well as adding a new function that does the matching.

In contrast, with the proposal of PR #132 a new basin enclosure method (or any other matching method) does not have to change the continuation or matching functions. If you have a look at the declared API of SSSetMatcher, all you have to do is make a new subtype and implement replacement_map for it. And everything will "just work" in typical Julia generic code fashion.

Crucially, you don't have to implement match_sequentially! yourself. It is already implemented. This includes the "ghost"/"vanished" functionality as well! Something your current code in #133 doesn't have.

Here is how it should work:

first, the new matcher MatchByBasinEnclosure should referrence the dynamical system, as it needs to re-init it constantly, and also change its parameter. The matcher should also reference prange and pidx, both of which are currently arguments to your function. Actually, the MatchByBasinEnclosure type stores everything that is currently an input to your function except atts_all. It then implements replacement_map, using the keyword argument i to access prange[i], to reinit the ds at that parameter. The extra keyword next_id can be used for attractors that cannot be matched to any of the previous attractors.

The whole

        for (k_future_int, att_future_int) in atts_future_integrated

loop is unecessary, as this is taken care of by match_sequentially!. You only need to write the code that matches between one parameter slice and the next one.

Does this sound good to you?

[Datseris]

Also, there is no reason to have featurizer here. We have extract_attractors. The matcher takes as input the found attractors. Any mapper that finds the actual attractors should work with this matching process. Why limit it to featurizing?

[KalelR]

In contrast, with the proposal of PR #132 a new basin enclosure method (or any other matching method) does not have to change the continuation or matching functions

How would you manage to keep track of the endpoint-integrated trajectories without changing the continuation code, though? I'd be very glad if you find a better way! I just don't see it so far. I did implement a modification of the continuation that did this, but that was very involved and not worth it. The nice part about the separate function is that the roles are more well separated.

[Datseris]

Ah, now I see the problem. THere is a fundamental difference in doing this for recurrences and for featurizing.

[Datseris]

for recurrences you can just literally map the endpoints to the attractors in the new parameter. so you can just use this as mapping...

[Datseris]

@KalelR at the code you pasted above what is atts_all supposed to be?

[KalelR]

for recurrences you can just literally map the endpoints to the attractors in the new parameter. so you can just use this as mapping...

sorry, how? I don't get it.

[KalelR]

@KalelR at the code you pasted above what is atts_all supposed to be?

The vector of dictionaries returned by ´continuation. I believe now we'd call it attractors_cont`.

[Datseris]

can you do a video call today? I can show you what I mean

[Datseris]

@KalelR I've pushed the basin enclosure continuation in src/continuation/continuation_basin_enclosure.jl in branch https://github.com/JuliaDynamics/Attractors.jl/blob/matcher_interface/src/continuation/continuation_basin_enclosure.jl

Do you mind branch off my branch now, to finish this, while I simplify all the "matching" new interface that I've created that is not really necessary?

[KalelR]

@KalelR I've pushed the basin enclosure continuation in src/continuation/continuation_basin_enclosure.jl in branch https://github.com/JuliaDynamics/Attractors.jl/blob/matcher_interface/src/continuation/continuation_basin_enclosure.jl

Do you mind branch off my branch now, to finish this, while I simplify all the "matching" new interface that I've created that is not really necessary?

Ok! Are you in a hurry about finishing this part?

[Datseris]

it would be nice to have it by tuesday, I am presenting attractors.jl then.

[Datseris]

@KalelR it's okay; this doesn't block anything. The initial source code is there but it is not exported or discussed in any docs, so you can take your time. I can go ahead and finish and merge #132 without requiring the basin enclosure to be complete.

We can also think of a better name than "basin enclosure".