Interface redesign proposal: matchers and a generic continuation type

[Datseris]

This proposal is motivated by two things:

1. The currently open PR #126
2. Some offline discussion I had with Kalel about a new way to mattch attractors, that is fundamentally different from the way we match currently. It is so fundamentally different, that one cannot just adjust the "distance function" in our current match_statespacesets! function so that the same result is achieved.

I want to propose a new interface, that will naturally make Attractors.jl easier to extend both in continuation algorithms, and in matching algorithms.

Matchers

The first step is a complete separation of continuation and matching. In fact, this is already achieved in the source code, but not at a formal interface level.

I propose a new abstract type SSSetMatcher (state space set, a.k.a. attractor, matcher). Subtypes of this interface extend the function

function match_continuation!(attractors_info, matcher::SSSetMatcher)
    # do your thing, store the replacement maps
    return rmaps::Vector{Dict}
end

This function operates given the "attractors" (vector of dictionaries). The extra convenience function

function match_continuation!(continuation_quantity::Vector{Dict}, rmaps::Vector{Dict})
    for (quantity, rmap) in zip(continuation_quantity, rmaps)
        swap_dict_keys!(quantity, rmap)
    end
    return
end

can be applied instantly to the basins fractions curves, or any other measure of attractors that we track, such as those in the example of PR #129 .

This new interface allows us to have the current matcher as one type, and the new matcher that Kalel has designed (but is not yet pushed anywhere online here) as a second matcher.

New generic continuation type

This orthogonality of matchers and mappers means that we can generalize the current RecurrencesFindAndMatch to a more general AttractorsContinueAndMatch. An alternative name could be AttractorsSeedAndMatch or AttractorsSeedContinueMatch. let me know which one you prefer.

Here is its docstring:

AttractorsContinueAndMatch(mapper, matcher [, seeding])

A continuation method for continuation. mapper is any subtype of AttractorMapper which implements extract_attractors, while matcher is any subtype of SSSetMatcher.

Description

This global continuation method has a simple basis. At the first parameter slice of the continuation process, attractors and their fractions are found using the given mapper. At each subsequent parameter slice, new attractors are found by seeding initial conditions from the previously found attractors and then running these initial conditions through the mapper. Seeding initial conditions close to previous attractors increases the probability that if an attractor continues to exist in the new parameter, it is found. Seeding is controlled by the seeding input. It is a function that given a state space set (an attractor) it returns an iterator of initial conditions sampled from the attractor. This generates the seeds.

After the special seeded initial conditions are mapped to attractors, attractor basin fractions are computed by sampling additional initial conditions using the provided ics in continuation). I.e., exactly as in basins_fractions. Naturally, during this step new attractors may be found, besides those found using the "seeding from previous attractors". Once the basins fractions are computed, the parameter is incremented again and we perform the steps as before.

This process continues for all parameter values. After all parameters are exhausted, the found attractors (and their fractions) are "matched" to the previous ones. This means: their IDs are changed, so that attractors that are "similar" to those at a previous parameter get assigned the same ID. Matching is done by the provided matcher. In fact, matching is a rather trivial call to match_continuation! at the end of the continuation function. This call is almost free. If you don't like the final matching output, you may use arbitrarily different matcher in match_continuation! without having to recompute the whole continuation. That is also why how matching works is described in the docstrings of each matcher.

and here is its source code:

function continuation(acam::AttractorsContinueAndMatch, prange, pidx, ics;
        samples_per_parameter = 100, show_progress = true,
    )
    progress = ProgressMeter.Progress(length(prange);
        desc = "Continuing attractors and basins:", enabled=show_progress
    )
    mapper = acam.mapper
    reset_mapper!(mapper)
    # first parameter is run in isolation, as it has no prior to seed from
    set_parameter!(referenced_dynamical_system(mapper), pidx, prange[1])
    fs = basins_fractions(mapper, ics; show_progress = false, N = samples_per_parameter)
    # At each parmaeter `p`, a dictionary mapping attractor ID to fraction is created.
    fractions_curves = [fs]
    # The attractors are also stored (and are the primary output)
    prev_attractors = deepcopy(extract_attractors(mapper))
    attractors_info = [prev_attractors]
    ProgressMeter.next!(progress; showvalues = [("previous parameter", prange[1]),])
    # Continue loop over all remaining parameters
    for p in prange[2:end]
        set_parameter!(referenced_dynamical_system(mapper), pidx, p)
        reset_mapper!(mapper)
        # Seed initial conditions from previous attractors
        # Notice that one of the things that happens here is some attractors have
        # really small basins. We find them with the seeding process here, but the
        # subsequent random sampling in `basins_fractions` doesn't. This leads to
        # having keys in `mapper.bsn_nfo.attractors` that do not exist in the computed
        # fractions. The fix is easy: we add the initial conditions mapped from
        # seeding to the fractions using an internal argument.
        seeded_fs = Dict{Int, Int}()
        for att in values(prev_attractors)
            for u0 in acam.seeding(att)
                # We map the initial condition to an attractor, but we don't care
                # about which attractor we go to. This is just so that the internal
                # array of `AttractorsViaRecurrences` registers the attractors
                label = mapper(u0; show_progress = false)
                seeded_fs[label] = get(seeded_fs, label, 0) + 1
            end
        end
        # Now perform basin fractions estimation as normal, utilizing found attractors
        # (the function comes from attractor_mapping.jl)
        fs = basins_fractions(mapper, ics;
            additional_fs = seeded_fs, show_progress = false, N = samples_per_parameter
        )
        # We do not match attractors here; the matching is independent step done at the end
        current_attractors = deepcopy(extract_attractors(mapper))
        push!(fractions_curves, fs)
        push!(attractors_info, current_attractors)
        overwrite_dict!(prev_attractors, current_attractors)
        ProgressMeter.next!(progress; showvalues = [("previous parameter", p),])
    end
    # Match attractors (and basins)
    rmaps = match_continuation!(attractors, matcher)
    match_continuation!(rmaps, fractions_curves)
    return fractions_curves, attractors_info
end

---

As you can imagine, i have actually already implemented this proposal. Not only it works immediatelly, but infact, it makes PR #126 obsolete. The existing AttractorsViaFeaturizing already directly works with this code, we just need to ensure that extract_attractors throws an error if ics is not a state space set but a generating function.

---

Thoughts? If you approve I can put in the PR and we can merge this, so that Kalel can go ahead and do the new matching PR.

The continuation function remains exactly as is. So is the current "grouping continuation".

[Datseris]

cc @awage @KalelR .