Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Merged
Merged
Show file tree
Hide file tree
Changes from 6 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
@@ -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`.
Expand Down
6 changes: 3 additions & 3 deletions Project.toml
Original file line number Diff line number Diff line change
@@ -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"
Expand All @@ -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"
Expand Down
43 changes: 22 additions & 21 deletions src/stability/fixedpoints.jl
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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`.
Expand All @@ -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).

"""
Expand All @@ -60,22 +60,21 @@ function fixedpoints(ds::DynamicalSystem, box, J = nothing;
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!
Expand All @@ -85,11 +84,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
Expand Down Expand Up @@ -132,7 +133,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
Expand Down
14 changes: 7 additions & 7 deletions test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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")

Expand Down
16 changes: 8 additions & 8 deletions test/stability/fixedpoints.jl
Original file line number Diff line number Diff line change
@@ -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]
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The problem is here: the initial vector must have the same type as the output of the function.

Since henon_rule and henon_jacob use SArrays, this should be a SVector (and same for the two other testsets below).

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

thank you!!! That fixed it!!!


# henon fixed points analytically
hmfp(a, b) = (x = -(1-b)/2a - sqrt((1-b)^2 + 4a)/2a; return SVector(x, b*x))
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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(ρ, β) = [
Expand Down
Loading