CompatHelper: bump compat for IntervalRootFinding to 0.6, (keep existing compat) (#338)

* CompatHelper: bump compat for IntervalRootFinding to 0.6, (keep existing compat)

* address brealing changes

* force the new version

* update more interval arithmetic breakages

* also update DSP compat

* organize code better, still breaking

* fixed fixed points!!!

---------

Co-authored-by: CompatHelper Julia <compathelper_noreply@julialang.org>
Co-authored-by: Datseris <datseris.george@gmail.com>

Author: 3 people

CHANGELOG.md

@@ -1,5 +1,10 @@
 # main
 
+# v3.3
+
+- Updated IntervalRootFinding.jl to v0.6 which [has breaking changes](https://github.com/JuliaIntervals/IntervalRootFinding.jl/releases/tag/v0.6.0) regarding
+  the `fixedpoints` function (no more `IntervalBox` type).
+
 # v3.2
  - `periodicorbits` uses new internal datastructure to avoid duplicates in the output. New keyword argument `abstol` is introduced and keyword argument `roundtol` is deprecated.
  - `fixedpoints` can now use `order` keyword argument to compute fixed points of n-th order iterations of `DeterministicIteratedMap`.

Project.toml

@@ -1,7 +1,7 @@
 name = "ChaosTools"
 uuid = "608a59af-f2a3-5ad4-90b4-758bdf3122a7"
 repo = "https://github.com/JuliaDynamics/ChaosTools.jl.git"
-version = "3.2.1"
+version = "3.3.0"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
@@ -25,12 +25,12 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 
 [compat]
 Combinatorics = "0.7, 1"
-DSP = "0.6, 0.7"
+DSP = "0.6, 0.7, 0.8"
 DataStructures = "0.18"
 Distances = "0.7, 0.8, 0.9, 0.10"
 Distributions = "0.21, 0.22, 0.23, 0.24, 0.25"
 DynamicalSystemsBase = "3"
-IntervalRootFinding = "0.5.9"
+IntervalRootFinding = "0.6"
 LombScargle = "1"
 Neighborhood = "0.2.2"
 Optim = "1.7"

src/stability/fixedpoints.jl

@@ -1,6 +1,6 @@
 import IntervalRootFinding, LinearAlgebra
-using IntervalRootFinding: (..), (×), IntervalBox, interval
-export fixedpoints, .., ×, IntervalBox, interval
+using IntervalRootFinding: (..), (×), interval
+export fixedpoints, .., ×, interval
 
 """
     fixedpoints(ds::CoreDynamicalSystem, box, J = nothing; kwargs...) → fp, eigs, stable
@@ -12,11 +12,11 @@ Fixed points are returned as a [`StateSpaceSet`](@ref).
 For convenience, a vector of the Jacobian eigenvalues of each fixed point, and whether
 the fixed points are stable or not, are also returned.
 
-`box` is an appropriate `IntervalBox` from IntervalRootFinding.jl.
+`box` is an appropriate interval box from IntervalRootFinding.jl (a vector of intervals).
 E.g. for a 3D system it would be something like
 ```julia
-v, z = -5..5, -2..2   # 1D intervals, can use `interval(-5, 5)` instead
-box = v × v × z       # `\\times = ×`, or use `IntervalBox(v, v, z)` instead
+v, z = interval(-5, 5), interval(-2, 2)  # 1D intervals
+box = [v, v, z]
 ```
 
 `J` is the Jacobian of the dynamic rule of `ds`.
@@ -37,14 +37,14 @@ as a start of a continuation process. See also [`periodicorbits`](@ref).
   see the docs of IntervalRootFinding.jl for all possibilities.
 - `tol = 1e-15` is the root-finding tolerance.
 - `warn = true` throw a warning if no fixed points are found.
-- `order = nothing` search for fixed points of the n-th iterate of 
+- `order = nothing` search for fixed points of the n-th iterate of
   [`DeterministicIteratedMap`](@ref). Must be a positive integer or `nothing`.
   Select `nothing` or 1 to search for the fixed points of the original map.
 
 ## Performance notes
 
-Setting `order` to a value greater than 5 can be very slow. Consider using 
-more suitable algorithms for periodic orbit detection, such as 
+Setting `order` to a value greater than 5 can be very slow. Consider using
+more suitable algorithms for periodic orbit detection, such as
 [`periodicorbits`](@ref).
 
 """
@@ -55,27 +55,30 @@ function fixedpoints(ds::DynamicalSystem, box, J = nothing;
     if isinplace(ds)
         error("`fixedpoints` currently works only for out-of-place dynamical systems.")
     end
+    if box isa Vector
+        # since the code only works for out of place we can always convert
+        box = SVector{length(box)}(box)
+    end
 
     if !(isnothing(order) || (isa(order, Int) && order > 0))
         error("`order` must be a positive integer or `nothing`.")
     end
 
-    # Jacobian: copy code from `DynamicalSystemsBase`
-    f = dynamic_rule(ds)
     p = current_parameters(ds)
-    if isnothing(J)
-        Jf(u, p, t) = DynamicalSystemsBase.ForwardDiff.jacobian(x -> f(x, p, 0.0), u)
-    else
-        Jf = J
-    end
     # Find roots via IntervalRootFinding.jl
     fun = to_root_f(ds, p, order)
-    jac = to_root_J(Jf, ds, p, order)
-    r = IntervalRootFinding.roots(fun, jac, box, method, tol)
+    jac = to_root_J(J, ds, p, order)
+    r = IntervalRootFinding.roots(fun, box; abstol = tol, derivative = jac, contractor = method)
     D = dimension(ds)
     fp = roots_to_dataset(r, D, warn)
-    # Find eigenvalues and stability
+    # extract eigenvalues and stability
     eigs = Vector{Vector{Complex{Float64}}}(undef, length(fp))
+    if isnothing(J)
+        f = dynamic_rule(ds)
+        Jf(u, p, t = 0) = DynamicalSystemsBase.ForwardDiff.jacobian(x -> f(x, p, t), u)
+    else
+        Jf = J
+    end
     Jm = zeros(dimension(ds), dimension(ds)) # `eigvals` doesn't work with `SMatrix`
     for (i, u) in enumerate(fp)
         Jm .= Jf(u, p, 0) # notice that we use the "pure" jacobian, no -u!
@@ -85,11 +88,13 @@ function fixedpoints(ds::DynamicalSystem, box, J = nothing;
     return fp, eigs, stable
 end
 
+to_root_J(Jf::Nothing, ::DynamicalSystem, args...) = nothing
+
 to_root_f(ds::CoupledODEs, p, ::Nothing) = u -> dynamic_rule(ds)(u, p, 0.0)
-to_root_J(Jf, ::CoupledODEs, p, ::Nothing) = u -> Jf(u, p, 0.0)
+to_root_J(Jf::Function, ::CoupledODEs, p, ::Nothing) = u -> Jf(u, p, 0)
 
 to_root_f(ds::DeterministicIteratedMap, p, ::Nothing) = u -> dynamic_rule(ds)(u, p, 0) - u
-function to_root_J(Jf, ds::DeterministicIteratedMap, p, ::Nothing)
+function to_root_J(Jf::Function, ds::DeterministicIteratedMap, p, ::Nothing)
     c = Diagonal(ones(typeof(current_state(ds))))
     return u -> Jf(u, p, 0) - c
 end
@@ -132,7 +137,7 @@ function roots_to_dataset(r, D, warn)
     end
     F = zeros(length(r), D)
     for (j, root) in enumerate(r)
-        F[j, :] .= map(i -> (i.hi + i.lo)/2, root.interval)
+        F[j, :] .= map(i -> (i.bareinterval.hi + i.bareinterval.lo)/2, root.region)
     end
     return StateSpaceSet(F; warn = false)
 end

test/runtests.jl

@@ -18,18 +18,18 @@ defaultname(file) = uppercasefirst(replace(splitext(basename(file))[1], '_' => '
 
 @testset "ChaosTools" begin
 
-include("timeevolution/orbitdiagram.jl")
+testfile("timeevolution/orbitdiagram.jl")
 
 testfile("chaosdetection/lyapunovs.jl")
 testfile("chaosdetection/gali.jl")
-include("chaosdetection/partially_predictable.jl")
-include("chaosdetection/01test.jl")
+testfile("chaosdetection/partially_predictable.jl")
+testfile("chaosdetection/01test.jl")
 testfile("chaosdetection/expansionentropy.jl")
 
-include("stability/fixedpoints.jl")
-include("periodicity/periodicorbits.jl")
-include("periodicity/davidchacklai.jl")
-include("periodicity/period.jl")
+testfile("stability/fixedpoints.jl")
+testfile("periodicity/periodicorbits.jl")
+testfile("periodicity/davidchacklai.jl")
+testfile("periodicity/period.jl")
 
 # testfile("rareevents/return_time_tests.jl", "Return times")
 

test/stability/fixedpoints.jl

@@ -1,12 +1,12 @@
 using ChaosTools, Test
 
 @testset "Henon map" begin
-    henon_rule(x, p, n) = SVector{2}(1.0 - p[1]*x[1]^2 + x[2], p[2]*x[1])
-    henon_jacob(x, p, n) = SMatrix{2,2}(-2*p[1]*x[1], p[2], 1.0, 0.0)
+    henon_rule(x, p, n = 0) = SVector{2}(1.0 - p[1]*x[1]^2 + x[2], p[2]*x[1])
+    henon_jacob(x, p, n = 0) = SMatrix{2,2}(-2*p[1]*x[1], p[2], 1.0, 0.0)
     ds = DeterministicIteratedMap(henon_rule, zeros(2), [1.4, 0.3])
     x = interval(-1.5, 1.5)
     y = interval(-0.5, 0.5)
-    box = x × y
+    box = [x, y]
 
     # henon fixed points analytically
     hmfp(a, b) = (x = -(1-b)/2a - sqrt((1-b)^2 + 4a)/2a; return SVector(x, b*x))
@@ -45,7 +45,7 @@ end
     ds = DeterministicIteratedMap(henon_rule, zeros(2), [1.4, 0.3])
     x = interval(-3.0, 3.0)
     y = interval(-10.0, 10.0)
-    box = x × y
+    box = [x, y]
     o = 4
     fp, eigs, stable = fixedpoints(ds, box, henon_jacob; order=o)
     tol = 1e-8
@@ -74,10 +74,10 @@ end
         end
     end
 
-    x = -20..20
-    y = -20..20
-    z = 0.0..40
-    box = x × y × z
+    x = interval(-20, 20)
+    y = interval(-20, 20)
+    z = interval(0, 40)
+    box = [x, y, z]
     σ = 10.0
     β = 8/3
     lorenzfp(ρ, β) = [