CompatHelper: bump compat for IntervalArithmetic to 0.22, (keep existing compat) (#129)

* CompatHelper: bump compat for IntervalArithmetic to 0.22, (keep existing compat)

* fix intervals_to_box to work with IntervalArithmatics 0.22

* fix some of the tests

* replace CoupledSDEsBasin extension with AttractorsExt

* fixed import

* Update Project.toml

* format

* fix docs make

* updated docs

---------

Co-authored-by: CompatHelper Julia <compathelper_noreply@julialang.org>
Co-authored-by: Orjan Ameye <orjan.ameye@hotmail.com>
Co-authored-by: reykboerner <reyk.boerner@reading.ac.uk>

Author: 4 people

Project.toml

@@ -33,7 +33,7 @@ StateSpaceSets = "40b095a5-5852-4c12-98c7-d43bf788e795"
 
 [extensions]
 ChaosToolsExt = ["ChaosTools"]
-CoupledSDEsBasin = ["ChaosTools", "Attractors"]
+AttractorsExt = ["Attractors"]
 DistMesh2D = ["DelaunayTriangulation"]
 
 [compat]
@@ -45,7 +45,7 @@ DocStringExtensions = "0.9.3"
 DynamicalSystemsBase = "3.13"
 Format = "1"
 Interpolations = "0.15.1"
-IntervalArithmetic = "0.20"
+IntervalArithmetic = "0.22"
 LinearSolve = "2.37.0, 3"
 Optimization = "4"
 OptimizationOptimisers = "0.3"

docs/make.jl

@@ -48,7 +48,7 @@ makedocs(;
     modules=[
         CriticalTransitions,
         Base.get_extension(CriticalTransitions, :ChaosToolsExt),
-        Base.get_extension(CriticalTransitions, :CoupledSDEsBasin),
+        Base.get_extension(CriticalTransitions, :AttractorsExt),
         Base.get_extension(DynamicalSystemsBase, :StochasticSystemsBase),
         # DynamicalSystemsBase
     ],

docs/src/examples/gMAM_Maierstein.md

@@ -127,24 +127,6 @@ lines!(ax2, reduce(hcat, tr), linewidth=1, color=(:black, 0.2))
 fig
 ```
 
-## Basins of attraction
-
-Basins of attraction are the regions in the state space that lead to a particular attractor. We can find the basins of attraction using the function `basins`.
-
-```@example GMAM
-using Attractors
-
-ba = basins(sys, [0.0, 0], [0.0, 1], [1.0, 0], intervals_to_box([-2, -2], [2, 2]), bstep=[0.01, 0.01], ϵ_mapper=0.001, Ttr=100)
-Ur, Vr, atr, M = ba
-heatmap(Ur, Vr, M)
-```
-
-The basin boundaries can be quickly extracted using the function `basin_boundaries`.
-
-```@example GMAM
-bb = basinboundary(ba)
-```
-
 ## Transitions
 
 We can quickly find a path which computes a transition from one attractor to another using the function `transition.

docs/src/man/systemanalysis.md

@@ -15,11 +15,4 @@ introduced by [Battelino et al. (1988)](https://doi.org/10.1016/0167-2789(88)900
 ```@docs
 edgetracking
 bisect_to_edge
-```
-
-## Baisin of attraction
-```@docs
-basinboundary
-basboundary
-basins
 ```

docs/src/man/utils.md

@@ -4,15 +4,15 @@
 
 ```@docs
 solver
+drift
+div_drift
 StochasticSystemsBase.covariance_matrix
 StochasticSystemsBase.diffusion_matrix
 noise_process
-drift
-div_drift
 ```
 
 ## More
 
 ```@docs
-intervals_to_box
+CriticalTransitions.intervals_to_box
 ```

ext/AttractorsExt.jl

@@ -0,0 +1,58 @@
+module AttractorsExt
+
+using CriticalTransitions
+using Attractors: Attractors, AttractorsViaProximity, AttractorMapper
+import Attractors: edgetracking, bisect_to_edge
+using DynamicalSystemsBase: ParallelDynamicalSystem
+
+export edgetracking, bisect_to_edge
+
+"""
+    edgetracking(ds::CoupledSDEs, attractors::Dict; diffeq, kwars...)
+
+Runs the edge tracking algorithm for the deterministic part of the CoupledSDEs `ds`.
+
+## Keyword arguments
+- `diffeq=(;alg = Vern9(), reltol=1e-11)`: ODE solver settings
+- `kwargs...`: all keyword arguments of `Attractors.edgetracking`
+"""
+function edgetracking(
+    ds::CoupledSDEs, attractors::Dict; diffeq=(; alg=Vern9(), reltol=1e-11), kwargs...
+)
+    return Attractors.edgetracking(CoupledODEs(ds; diffeq), attractors; kwargs...)
+end
+
+"""
+    bisect_to_edge(ds::CoupledSDEs, attractors::Dict;
+        u1, u2, bisect_thresh, diffeq, verbose, kwargs...)
+
+Runs the `bisect_to_edge` function for the deterministic part of the CoupledSDEs `ds`.
+
+## Keyword arguments
+- `u1`: State 1 (default: first state in `attractors`)
+- `u2`: State 2 (default: second state in `attractors`)
+- `bisect_thresh=1e-6`: distance threshold
+- `diffeq=(;alg = Vern9(), reltol=1e-11)`: ODE solver settings
+- `verbose=false`: Verbosity of output
+- `ϵ_mapper=nothing`: ϵ argument of `Attractors.AttractorsViaProximity`
+- `kwargs...`: Keyword arguments passed to `Attractors.AttractorsViaProximity`
+"""
+function bisect_to_edge(
+    ds::CoupledSDEs,
+    attractors::Dict;
+    u1=collect(values(attractors))[1][1],
+    u2=collect(values(attractors))[2][1],
+    bisect_thresh=1e-6,
+    diffeq=(; alg=Vern9(), reltol=1e-11),
+    ϵ_mapper=nothing,
+    verbose=false,
+    kwargs...,
+)
+    odes = CoupledODEs(ds; diffeq)
+    pds = Attractors.ParallelDynamicalSystem(odes, [u1, u2])
+    mapper = Attractors.AttractorsViaProximity(odes, attractors, ϵ_mapper; kwargs...)
+
+    return Attractors.bisect_to_edge(pds, mapper; bisect_thresh, verbose)
+end
+
+end # module AttractorsExt

ext/ChaosToolsExt.jl

@@ -20,16 +20,10 @@ Generates a box from specifying the interval limits in each dimension.
 """
 function CriticalTransitions.intervals_to_box(bmin::Vector, bmax::Vector)
     # Generates a box from specifying the interval limits
-    intervals = []
-    dim = length(bmin)
-    for i in 1:dim
-        push!(intervals, interval(bmin[i], bmax[i]))
+    if length(bmin) != length(bmax)
+        @warn "bmin and bmax must have the same length."
     end
-    box = intervals[1]
-    for i in 2:dim
-        box = IntervalArithmetic.cross(box, intervals[i])
-    end
-    return box
+    return SVector{length(bmin)}([interval(bmin[i], bmax[i]) for i in 1:length(bmin)])
 end
 
 """
@@ -65,49 +59,4 @@ function ChaosTools.fixedpoints(sys::CoupledSDEs, box)
     return fixedpoints(CoupledODEs(sys), box, jacobian(ds))
 end
 
-# function saddles_idx(fps::Tuple)
-#     num = size(fps[1],1); # number of fixed points
-#     dim = size(fps[1],2); # dimension of the system
-#     eigenvalues = fps[2];
-#     idx = [false for i ∈ 1:num];
-#     for ii ∈ 1:num
-#         imag_parts = [imag(eigenvalues[ii][jj]) for jj ∈ 1:dim]
-#         if all(imag_parts.==0) # we have purely real eigenvalues
-#             real_parts = [real(eigenvalues[ii][jj]) for jj ∈ 1:dim];
-#             if prod(real_parts) < 0 # we have at least positive eigenvalue and at least one negative eigenvalue
-#                 idx[ii] = true;
-#             end
-#         end
-#     end
-#     idx
-# end
-
-# function repellers_idx(fps::Tuple)
-#     num = size(fps[1],1); # number of fixed points
-#     dim = size(fps[1],2); # dimension of the system
-#     eigenvalues = fps[2];
-#     idx = [false for i ∈ 1:num];
-#     for ii ∈ 1:num
-#         real_parts = [real(eigenvalues[ii][jj]) for jj ∈ 1:dim];
-#         if all(real_parts .> 0)
-#             idx[ii] = true;
-#         end
-#     end
-#     idx
-# end
-
-# function attractors_idx(fps::Tuple)
-#     num = size(fps[1],1); # number of fixed points
-#     dim = size(fps[1],2); # dimension of the system
-#     eigenvalues = fps[2];
-#     idx = [false for i ∈ 1:num];
-#     for ii ∈ 1:num
-#         real_parts = [real(eigenvalues[ii][jj]) for jj ∈ 1:dim];
-#         if all(real_parts .< 0)
-#             idx[ii] = true;
-#         end
-#     end
-#     idx
-# end
-
 end