Rename: `continuation` to `global_continuation`

[Datseris]

I don't remember well if I have shared the news, but @JonasKoziorek is a GSOC student currently working on a new submodule for DynamicalSystems.jl : https://github.com/JuliaDynamics/PeriodicOrbits.jl . The same project can also serve as a first step for a common interface between BifurcationKit.jl and Attractors.jl.

After spending more time re-reading BifurcationKit.jl, and re-doing the comparison we have done in our 2023 CHAOS paper, I am convinced that it is possible to re-write / interface BifurcationKit.jl to be based on DynamicalSystems.jl and integrate with the ecosystem more, and also share the continuation interface. That is the bright unified future I have planned in my head, but alas it would take too much effort from me to make all the changes in the BifurcationKit.jl code.

In any case, as we know, the continuation in each software is fundamentally different. That is why I propose to rename our continuation to global_continuation. "global" because it views the whole state space (within a specified box). BifurcationKit.jl may perhaps rename its continuation to local_continuation, "local" because it tracks a specific point in the state space (or in the periodic orbit), and perhaps more importantly because it computes a local stability indicator (in particular the Jacobian). I don't know whether @rveltz is interested in renaming the continuation function, but in the bright unified future I have dreamed of having these two different names would be better for the user.

cc @awage @KalelR

[Datseris]

cc @rveltz I spelt your GitHub tag wrong in the original post.

[rveltz]

Another way would be to define a new continuation algorithm GLOBALCONT (subtype of abstract type of BK.AbstractContinuationAlgorithm) and define in DS a method continuation(prob, GLOBALCONT(),...). But that would mean to depend on BK. Perhaps, I (we?) should define a ContinuationBase.jl which regroups those types for easier/smaller dependency?

There is also the "global" continuation BK.DefCont in BifurcationKit

[Datseris]

Another way would be to define a new continuation algorithm GLOBALCONT (subtype of abstract type of BK.AbstractContinuationAlgorithm) and define in DS a method continuation(prob, GLOBALCONT(),...)

I am not sure how much sense this would make from an interface point of view. Different algorithms should be part of the same interface, here continuation, if the output is of the same type. The output of Attractors.continuation and BifurcationKit.continuation is different at a very fundamental level:

While for BifurcationKit.jl I found: https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/library/#BifurcationKit.continuation

In fact, BifurcationKit.jl does not appear to define a continuation interface. I found at least 5 docstrings for continuation, as it appears that different algorithms have different inputs to the continuation function itself via varying keyword arguments (correct me if I am wrong). There is also the "resume continuation" that is also called continuation, which does a branch switching from a previously found bifurcation point. In Attractor.continuation the input and output is always the same irrespectively of the algorithm (all options are given to the algorithm type itself).

That is why I don't see any obvious benefit, for neither the user nor the developer, for this type of output to be in the same interface continuation, hence the suggestion of global_ or local_.

---

There is also the "global" continuation BK.DefCont in BifurcationKit

Oh, I didn't know about this. This looks very promising, yet one more reason to unify the two approaches. What does continuation output in with this algorithm? It can't be the typical output that is "one branch", right? Or is the same struct that represents one branch also used to represent several branches? I also wonder what are the capabilities of this method. Is it limited only to fixed points (stable and unstable)? How robust is it? If there is coexistence of various fixed points and limit cycles, typical in Kuramoto networks, would this method at least find the fixed points?

[rveltz]

output in with this algorithm?

DCResult

In fact, BifurcationKit.jl does not appear to define a continuation interface. I found at least 5 docstrings for continuation, as it appears that different algorithms have different inputs to the continuation function itself via varying keyword arguments (correct me if I am wrong)

Exact. You are right and this API is probably not the end of the story. I thought the current interface is easier to the user. Compare:

# compute branch of solutions with detection of bifurcations
br = continuation(prob, PALC(), options_cont)
# compute bifurcated branch from say the first bifurcation point which is a bp point
br_branched = continuation(br, 1)
# or continue the bifurcation point in two parameters, p2 being the additional parameter
br_br_codim2 = continuation(br, 2, (@lens _.p2)

with the other hypothetical idea that I have which is ala SciML, you find this API in the deprecated Bifurcation.jl

# compute branch of solutions with detection of bifurcations
br = continuation(prob, PALC(), options_cont)
# compute bifurcated branch from say the first bifurcation point which is a bp point
bp_prob = BranchFromBifurcation(br, 1)
br_branched = continuation(bp_prob, PALC())
# or continue the bifurcation point in two parameters, p2 being the additional parameter
br_codim2_prob = Codim2Problem(br, 1, (@lens _.p2))
br_br_codim2 = continuation(br_codim2_prob, PALC())

[Datseris]

Exact. You are right and this API is probably not the end of the story. I thought the current interface is easier to the user. Compare:

# compute branch of solutions with detection of bifurcations
br = continuation(prob, PALC(), options_cont)
# compute bifurcated branch from say the first bifurcation point which is a bp point
br_branched = continuation(br, 1)
# or continue the bifurcation point in two parameters, p2 being the additional parameter
br_br_codim2 = continuation(br, 2, (@lens _.p2)

Yes, this is nicer than the "SciML alternative" you quote below it. But I believe there is a version that is nicer than both. In particular:

# compute branch of solutions with detection of bifurcations
br = local_continuation(ds::DynamicalSystem, lca::LocalContinuationAlgorithm, pidx; cont_kwargs...)
# br references the `ds`. this gives the advantage that the interface becomes 
# decoupled from `ds` and would work for "any" dynamical system `lca` makes sense with.
# Then:
# compute bifurcated branch from say the first bifurcation point which is a bp point
br = continue_from_bfpoint!(br, 1) # adds information in-place to `br`
# or continue the bifurcation point in two parameters, p2 being the additional parameter
br_br_codim2 = continue_from_bf_point_2D(br, 2, p2idx)

br would reference ds and would store all branches in a Vector{ContinuationBranch}, and all bifurcation points in another Vector{BifurcationPoint}. A crucially important aspect of this approach is that cont_kwargs must not be related in any way to lca. They must be options that universaly apply to all local continuation, such as the parameter range, the maximum possible step to take in parameter axis, etc. All options that actually care about lca are given as inputs when you create lca. This eliminates the problem you have right now where continuation may take different keywords depending on the algorithm.

This approach is based on "each function does one thing, the thing the name suggests". Textbooks in good software design recommend this approach: don't use the same function name for different things (find branch, resume another branch, plotting). This also means that anything related to plotting is not given to local_continuation but to a new plotting function plot_continuation_branches or however we'd like to call it. THat is also why co-dim-2 continuation has a new name and returns a new type. Because the information returned is fundamentally dfferent, i.e., it can't be a Vector{ContinuationBranch}.

---

Off topic, but in DynamicalSystems.jl there is no reason for "lenses" to specify parameters, a construct which I find odd. You just give the index of the parameter container. Why did you use lenses at all? A beginner in Julia understands what is a mutable container and the index of the container to adjust its value. Lenses is an additional concept that so far I have not found a reason to use, instead of just having a mutable parameter container, which is also what SciML recommends for the ODE problem, and it is also what it creates if you have symbolic problems via MTK.

---

@KalelR and @awage , I am finishing my talk for JuliaCon today, and I would much rather refer to global_continuation rather than continuation if you are also okay with the name. I've thought about this extensively, and I am certain now that based on the new PeriodicOrbits.jl package, we can have both local and global continuations in harmony in DynamicalSystems.jl.

[awage]

I am finishing my talk for JuliaCon today, and I would much rather refer to global_continuation rather than continuation if you are also okay with the name.

Fine for me! local_continuation would be the interface to BK?

[rveltz]

Why did you use lenses at all? A beginner in Julia understands what is a mutable container and the index of the container to adjust its value. Lenses is an additional concept that so far I have not found a reason to use, instead of just having a mutable parameter container, which is also what SciML recommends for the ODE problem, and it is also what it creates if you have symbolic problems via MTK

It was a suggestion of @tkf. I think the best parameter container for bifurcation analysis is a named tuple but this is debatable. All old softwares auto, matcont, etc use index to specify the parameter but with the tooling of Julia, I think we deserve better. Look at here. You want to switch to continuation wrt U0? "hang on, what is the index". Additionally, it greatly simplifies plots where the parameter axis is correctly labelled. Even better: think about the case of 2 parameters continuation. If you dont use named tuples, you get generic p1, p2 as parameter axis.

I could add an interface like update_param!(pars, index, value) and use it all around in BK. Then it would be trivial to enable the array based parameter container (that you suggest) or the immutable one. (This is actually about to be pushed since I want to use Accessors instead of Setfield)

Why not a dict, or array then?

Because the compiler will remove the par.U0 or the change @set par.U0 = newpar so it costs nothing to be immutable (it does not allocate). I think same should be true for mutable (array) container used with with @lens _[1]. Also, if you use arrays, how do you give a name to the components?

which is also what SciML recommends for the ODE problem

where?

it is also what it creates if you have symbolic problems via MTK

Yes, I saw this. We discussed the same problem with @TorkelE but he manages to nicely propagate the bifurcation parameter in his interface. The interface MTK-BK is not fully tested yet. I like the MTK way to specify the continuation parameter but I can't expect all users to use MTK to specify their problem.

[rveltz]

On a side note, I think local_continuation is not a good name. I think what you have in mind is local_continuation -> connected_component and global_continuation -> all_connected_components

[Datseris]

I think what you have in mind is local_continuation -> connected_component and global_continuation -> all_connected_components

That is not entirely correct either unfortunately.

In our global_continuation you don't get just the connected components. You also get a measure of nonlocal stability (by default the relative basin fractions) which is achieved by globally sampling the state space. So its more like all_connected_components_and_nonlocal_stabilities.

In contrast, what you get in the "standard" output of traditional continuation software, with the exception of BK.DefCont, is single_connected_component_and_linear_stability.

So focusing only on whether you have one or multiple connected components is not an accurate reflection of the differences. That's why so far I prefer more local_ and global_ continuation, although I definitely would welcome alternative names that suite more the differences and are also simple.