[WIP] nonlocal stability measure estimator

docs/src/api.md

@@ -20,11 +20,13 @@ subdivision_based_grid
 ```
 
 ### Proximity
+
 ```@docs
 AttractorsViaProximity
 ```
 
 ### Featurizing
+
 ```@docs
 AttractorsViaFeaturizing
 ```
@@ -122,6 +124,13 @@ MFSBlackBoxOptim
 MFSBruteForce
 ```
 
+## Nonlocal Stability
+
+```@docs
+StabilityMeasuresAccumulator
+finalize_accumulator
+```
+
 ## Global continuation
 
 ```@docs

ext/plotting.jl

@@ -295,15 +295,15 @@ function Attractors.plot_continuation_curves!(ax, continuation_info, prange = 1:
         labels = Dict(ukeys .=> ukeys),
         add_legend = length(ukeys) < 7,
         markers = markers_from_keys(ukeys),
-        slines_kwargs = (;linewidth = 3,),
+        series_kwargs = NamedTuple(),
         axislegend_kwargs = (position = :lt,)
     )
 
     series = continuation_series(continuation_info, NaN, ukeys)
     for k in ukeys
         scatterlines!(ax, prange, series[k];
             color = colors[k], label = "$(labels[k])", marker = markers[k],
-            markersize = 10, linewidth = 3, slines_kwargs...
+            markersize = 10, linewidth = 3, series_kwargs...
         )
     end
     xlims!(ax, minimum(prange), maximum(prange))
@@ -381,7 +381,6 @@ function Attractors.plot_basins_attractors_curves!(axb, axa, fractions_cont, att
     return
 end
 
-
 ##########################################################################################
 # Videos
 ##########################################################################################

src/continuation/continuation_ascm_generic.jl

@@ -83,7 +83,7 @@ without having to recompute the whole global continuation!
 
 ### Step 5 - Finish
 
-After matching the parameter is incremented.
+After matching the parameter(s) is incremented.
 Steps 1-4 repeat until all parameter values are exhausted.
 
 ### Further note

src/continuation/continuation_stability_measures.jl

@@ -0,0 +1,113 @@
+# TODO: Can we make this function not duplicate the source code
+# of the original `ASCM` implementation?
+
+function global_continuation(
+    acam::AttractorSeedContinueMatch{<:StabilityMeasuresAccumulator}, pcurve, ics;
+        samples_per_parameter = 100, show_progress = true,
+    )
+    N = samples_per_parameter
+    progress = ProgressMeter.Progress(length(pcurve); desc = "Continuing attractors and stability:", enabled=show_progress)
+    mapper = acam.mapper
+    reset_mapper!(mapper)
+
+    set_parameters!(referenced_dynamical_system(mapper), pcurve[1])
+    if ics isa Function
+        fs = basins_fractions(mapper, ics; show_progress = false, N = samples_per_parameter)
+    else # we ignore labels in this continuation algorithm
+        fs, = basins_fractions(mapper, ics; show_progress = false)
+    end
+
+    # this is the first difference with the standard continuation: we have to call finalize.
+    measures = finalize_accumulator(mapper)
+    measures_cont = [measures]
+
+    prev_attractors = deepcopy(extract_attractors(mapper))
+    attractors_cont = [deepcopy(prev_attractors)] # we need the copy
+    ProgressMeter.next!(progress)
+
+    # Continue loop over all remaining parameters
+    ### For @Datseris: Is this what you had in mind when we emailed about 
+    ### the global continuation algorithm with AttractorsViaProximity?
+    ### Currently, line 44 has an error because we cant pass ics to the AttractorsViaRecurrences.
+    for p in @view(pcurve[2:end])
+        set_parameters!(referenced_dynamical_system(mapper), p)
+        reset_mapper!(mapper)
+        if typeof(mapper.mapper) <: AttractorsViaRecurrences
+            fs = if allows_mapper_u0(mapper)
+                seed_attractors_to_fractions_individual(mapper, prev_attractors, ics, N, acam.seeding)
+            else
+                seed_attractors_to_fractions_grouped(mapper, prev_attractors, ics, N, acam.seeding)
+            end
+            current_attractors = deepcopy(extract_attractors(mapper))
+        else
+            error("Unsupported mapper type: $(typeof(mapper.mapper))")
+        end
+        push!(attractors_cont, current_attractors)
+        overwrite_dict!(prev_attractors, current_attractors)
+        ProgressMeter.next!(progress)
+
+        # this is the second difference in global continuation
+        measures = finalize_accumulator(mapper)
+        push!(measures_cont, measures)
+    end
+
+    rmaps = match_sequentially!(attractors_cont, acam.matcher; pcurve, ds = referenced_dynamical_system(mapper))
+
+    # This is the third difference in global continuation
+    for i in 2:length(pcurve)
+        for dict in values(measures_cont[i])
+            swap_dict_keys!(dict, rmaps[i-1])
+        end
+    end
+    transposed = Dict{String, Vector{Dict{Int64, Float64}}}()
+    for measure in measures_cont[1]
+        measure_name = measure[1]
+        transposed[measure_name] = Vector{Dict{Int64, Float64}}()
+    end
+
+    for measures in measures_cont
+        for (measure_name, measure_dict) in measures
+            push!(transposed[measure_name], measure_dict)
+        end
+    end
+    measures_cont = transposed
+
+    return measures_cont, attractors_cont
+end
+
+
+
+# make sure to allow the possiblity that the proximity options can also be
+# vectors of same length as `pcurve`; Same for the distributions
+function stability_measures_along_continuation(ds::DynamicalSystem, attractors_cont, ics, pcurve, εs, distributions, Ts; metric=Euclidean(), proximity_mapper_options=[])
+    measures_cont = []
+    for (i, p) in enumerate(pcurve)
+        ε = εs isa AbstractVector ? εs[i] : εs # if its a vector, get i-th entry
+        d = distributions isa AbstractVector ? distributions[i] : distributions # if its a vector, get i-th entry
+        T = Ts isa AbstractVector ? Ts[i] : Ts # if its a vector, get i-th entry
+        set_parameters!(ds, p)
+        attractors = attractors_cont[i]
+        accumulator = StabilityMeasuresAccumulator(
+            AttractorsViaProximity(ds, attractors, ε, proximity_mapper_options...),
+            d=d, T=T, metric=metric
+        )
+        for u0 in ics
+            id = accumulator(u0)
+        end
+        push!(measures_cont, finalize_accumulator(accumulator))
+    end
+
+    transposed = Dict{String, Vector{Dict{Int64, Float64}}}()
+    for measure in measures_cont[1]
+        measure_name = measure[1]
+        transposed[measure_name] = Vector{Dict{Int64, Float64}}()
+    end
+
+    for measures in measures_cont
+        for (measure_name, measure_dict) in measures
+            push!(transposed[measure_name], measure_dict)
+        end
+    end
+    measures_cont = transposed
+    return measures_cont
+end

src/continuation/global_continuation_api.jl

@@ -1,4 +1,4 @@
-export global_continuation, GlobalContinuationAlgorithm, continuation_series
+export global_continuation, GlobalContinuationAlgorithm, continuation_series, stability_measures_along_continuation
 
 """
     GlobalContinuationAlgorithm
@@ -25,6 +25,9 @@ Return:
 1. `fractions_cont::Vector{Dict{Int, Float64}}`. The fractions of basins of attraction.
    `fractions_cont[i]` is a dictionary mapping attractor IDs to their basin fraction
    at the `i`-th parameter combination.
+   -  This output is different if you are using [`StabilityMeasuresAccumulator`](@ref)
+      in combination with [`AttractorSeedContinueMatch`](@ref). See the docstring
+      of [`StabilityMeasuresAccumulator`](@ref) for more details.
 2. `attractors_cont::Vector{Dict{Int, <:Any}}`. The continued attractors.
    `attractors_cont[i]` is a dictionary mapping attractor ID to the
    attractor set at the `i`-th parameter combination.
@@ -98,4 +101,5 @@ end
 include("continuation_ascm_generic.jl")
 include("continuation_recurrences.jl")
 include("continuation_grouping.jl")
-include("aggregate_attractor_fractions.jl")
+include("aggregate_attractor_fractions.jl")
+include("continuation_stability_measures.jl")

src/mapping/attractor_mapping.jl

@@ -15,7 +15,9 @@ export AttractorMapper,
     extract_attractors,
     subdivision_based_grid,
     SubdivisionBasedGrid,
-    reset_mapper!
+    reset_mapper!,
+    StabilityMeasuresAccumulator,
+    finalize_accumulator
 
 ValidICS = Union{AbstractVector, Function}
 
@@ -48,6 +50,9 @@ The mappers that can do this are:
 * [`AttractorsViaRecurrences`](@ref)
 * [`AttractorsViaFeaturizing`](@ref) with the [`GroupViaHistogram`](@ref) configuration.
 
+See also [`StabilityMeasuresAccumulator`](@ref) that extends this interface
+to accelerate estimation of stability measures.
+
 ## For developers
 
 `AttractorMapper` defines an extendable interface. A new type needs to subtype
@@ -183,3 +188,4 @@ function convergence_time end
 include("attractor_mapping_proximity.jl")
 include("recurrences/attractor_mapping_recurrences.jl")
 include("grouping/attractor_mapping_featurizing.jl")
+include("stability_measures_accumulator.jl")

src/mapping/attractor_mapping_proximity.jl

@@ -92,6 +92,8 @@ function AttractorsViaProximity(ds::DynamicalSystem, attractors::Dict, ε = noth
     return mapper
 end
 
+reset_mapper!(mapper::AttractorsViaProximity) = nothing # AM can we add this for consistency?
+
 function _deduce_ε_from_attractors(attractors, search_trees, verbose = false)
     if length(attractors) != 1
         verbose && @info("Computing minimum distance between attractors to deduce `ε`...")