Update of StabilityMeasuresAccumulator (#184)

* Make possible that nls parameters are functions

* Check epsilon distance minimum before saving ct

* Corrected median conv time/pace and optimized critical shocks

* notes on how to simplify code

* Comprehensive refactoring to record only needed information

Has been tested to show that the results did not change

* Small fix for PR tests

* Adjusted saving of ics and corrected expected test values

* More type parametrization

* Fix of max amplification

* Handle Inf values in plotting of continuation curves

* Updated docs

* add some dummy entries to bib for now

* use `distace`

* Continuation with variable distance parameter

* Bump version

* fix name specification

* add instructions how to extent the accumulator

---------

Co-authored-by: Datseris <[email protected]>

Author: andreasmorr

CHANGELOG.md

@@ -1,3 +1,9 @@
+# v1.29
+
+- Keyword `metric` of `StabilityMeasuresAccumulator` is no longer used as it was incorrect.
+- New keyword `distance` for `StabilityMeasuresAccumulator`.
+- Source code of `StabilityMeasuresAccumulator` has been drastically simplified, and instructions for how to extent its functionality were added.
+
 # v1.28
 
 - `extract_attractors` will now automatically name the dimensions of the attractors for MTK-derived systems.

Project.toml

@@ -2,7 +2,7 @@ name = "Attractors"
 uuid = "f3fd9213-ca85-4dba-9dfd-7fc91308fec7"
 authors = ["George Datseris <[email protected]>", "Kalel Rossi", "Alexandre Wagemakers"]
 repo = "https://github.com/JuliaDynamics/Attractors.jl.git"
-version = "1.28.0"
+version = "1.29.0"
 
 
 [deps]

docs/refs.bib

@@ -296,12 +296,12 @@ @article{Krakovska2024ResilienceDynamicalSystems
 }
 
 @article{Morr2025,
-  title = {To be published},
-  author = {Andreas Morr, George Datseris},
-  journal = {},
-  volume = {},
-  number = {},
-  pages = {},
+  title = {Numerically Estimating Resilience Measures in Dynamical Systems},
+  author = {Andreas Morr and Christian Kuehn and George Datseris},
+  journal = {to be published},
+  volume = {1},
+  number = {1},
+  pages = {1},
   year = {2025},
 }
 

ext/plotting.jl

@@ -301,6 +301,10 @@ function Attractors.plot_continuation_curves!(ax, continuation_info, prange = 1:
     )
 
     series = continuation_series(continuation_info, NaN, ukeys)
+    for (k, v) in series
+        v[isinf.(v)] .= NaN
+    end
+    
     for k in ukeys
         scatterlines!(ax, prange, series[k];
             color = colors[k], label = "$(labels[k])", marker = markers[k],

src/continuation/continuation_stability_measures.jl

@@ -108,7 +108,7 @@ continuation steps.
 
 - `ε = nothing`: given to [`AttractorsViaProximity`](@ref).
 - `proximity_mapper_options = NamedTuple()`: extra keywords for `AttractorsViaProximity`.
-- `metric, finite_time, weighting_distribution`: given to [`StabilityMeasuresAccumulator`](@ref).
+- `distance, finite_time, weighting_distribution`: given to [`StabilityMeasuresAccumulator`](@ref).
 - `samples_per_parameter = 1000`: how many samples to use when estimating stability measures
   via [`StabilityMeasuresAccumulator`](@ref). Ignored when `ics` is not a function.
 """
@@ -121,7 +121,7 @@ function stability_measures_along_continuation(
         weighting_distribution = EverywhereUniform(),
         finite_time = 1.0,
         samples_per_parameter = 1000,
-        metric = Euclidean(),
+        distance = Centroid(),
         proximity_mapper_options = NamedTuple(),
         show_progress=true
     )
@@ -130,16 +130,41 @@ function stability_measures_along_continuation(
     )
     measures_cont = []
     for (i, p) in enumerate(pcurve)
-        ε_ = ε isa AbstractVector ? ε[i] : ε # if its a vector, get i-th entry
-        weighting_distribution_ = weighting_distribution isa AbstractVector ?
-                                weighting_distribution[i] : weighting_distribution
-        finite_time_ = finite_time isa AbstractVector ? finite_time[i] : finite_time
         set_parameters!(ds, p)
         attractors = attractors_cont[i]
+        if ε isa AbstractVector
+            ε_ = ε[i]
+        elseif ε isa Function
+            ε_ = ε(p, attractors)
+        else
+            ε_ = ε
+        end
+        if weighting_distribution isa AbstractVector
+            wd = weighting_distribution[i]
+        elseif weighting_distribution isa Function
+            wd = weighting_distribution(p, attractors)
+        else
+            wd = weighting_distribution
+        end
+        if finite_time isa AbstractVector
+            ft = finite_time[i]
+        elseif finite_time isa Function
+            ft = finite_time(p, attractors)
+        else
+            ft = finite_time
+        end
+        if distance isa AbstractVector
+            d = distance[i]
+        elseif distance isa Function
+            d = distance(p, attractors)
+        else
+            d = distance
+        end
+
         accumulator = StabilityMeasuresAccumulator(
             AttractorsViaProximity(ds, attractors; ε = ε_, proximity_mapper_options...);
-            weighting_distribution = weighting_distribution_, finite_time = finite_time_,
-            metric = metric
+            weighting_distribution=wd, finite_time=ft,
+            distance=d
         )
         N = ics isa Function ? samples_per_parameter : length(ics)
         for i ∈ 1:N

src/mapping/attractor_mapping_proximity.jl

@@ -69,6 +69,9 @@ function AttractorsViaProximity(ds::DynamicalSystem, attractors::Dict;
     if !(valtype(attractors) <: AbstractStateSpaceSet)
         error("The input attractors must be a dictionary with values of `StateSpaceSet`s.")
     end
+    if isempty(attractors)
+        error("The input attractors cannot be an empty dictionary.")
+    end
     if dimension(ds) ≠ dimension(first(attractors)[2])
         error("Dimension of the dynamical system and candidate attractors must match.")
     end