Return & transit time statistics

src/ChaosTools.jl

@@ -9,7 +9,8 @@ using Reexport
 using DynamicalSystemsBase: DS, DDS, CDS
 using DynamicalSystemsBase: MDI, TDI
 using DynamicalSystemsBase: stateeltype
-using DiffEqBase: AbstractODEIntegrator, u_modified!
+using DiffEqBase: AbstractODEIntegrator, u_modified!, DEIntegrator
+DEI = DEIntegrator
 
 include("orbitdiagrams/discrete_diagram.jl")
 include("orbitdiagrams/poincare.jl")
@@ -22,6 +23,7 @@ include("nlts.jl")
 
 include("period_return/periodic_points.jl")
 include("period_return/period.jl")
+include("period_return/transit_times_statistics.jl")
 
 include("chaosdetection/lyapunovs.jl")
 include("chaosdetection/gali.jl")

src/nlts.jl

@@ -38,7 +38,7 @@ Typically `R` is the result of delay coordinates of a single timeseries.
 If the dataset/reconstruction
 exhibits exponential divergence of nearby states, then it should clearly hold
 ```math
-E(k) \\approx \\lambda\\Delta t k + E(0)
+E(k) \\approx \\lambda\\cdot k \\cdot \\Delta t + E(0)
 ```
 for a *well defined region* in the `k` axis, where ``\\lambda`` is
 the approximated

src/orbitdiagrams/poincare.jl

@@ -1,5 +1,5 @@
 using DynamicalSystemsBase: DEFAULT_DIFFEQ_KWARGS, _get_solver
-using Roots: find_zero, A42, AlefeldPotraShi, Brent
+using Roots: find_zero, A42
 export poincaresos, produce_orbitdiagram, PlaneCrossing
 
 const ROOTS_ALG = A42()
@@ -135,7 +135,7 @@ function poincaresos(integ, planecrossing, tfinal, Ttr, j, rootkw)
         integ.t > tfinal + Ttr && break
 
         # I am now guaranteed to have `t` in negative and `tprev` in positive
-        tcross = find_zero(f, (integ.tprev, integ.t), ROOTS_ALG; rootkw...)
+        tcross = Roots.find_zero(f, (integ.tprev, integ.t), ROOTS_ALG; rootkw...)
         ucross = integ(tcross)
 
         push!(data, ucross[j])

src/period_return/period.jl

@@ -5,7 +5,7 @@ import LombScargle
 export estimate_period
 
 """
-    estimate_period(v, method, t=0:length(v)-1;, method_specific_kwargs...)
+    estimate_period(v::Vector, method, t=0:length(v)-1; kwargs...)
 Estimate the period of the signal `v`, with accompanying time vector `t`,
 using the given `method`.
 

src/period_return/transit_times_statistics.jl

@@ -0,0 +1,159 @@
+using LinearAlgebra
+using Roots
+using Distances
+
+# Shortcuts (already defined)
+# MDI = DynamicalSystemsBase.MinimalDiscreteIntegrator
+# DEI = DynamicalSystemsBase.DiffEqBase.DEIntegrator
+
+# TODO: Add example with rectangle, perhaps not clear
+
+export transit_time_statistics, transit_return
+
+"""
+    transit_time_statistics(ds::DynamicalSystem, u₀, εs, T; diffeq...) → exits, entries
+Collect transit time statistics for a ball/box centered at `u₀` with radii `εs` (see below),
+in the state space of the given dynamical system.
+Return the exit and (re-)entry return times to the set(s), where each of these is a vector
+containing all collected times for the respective `ε`-radius set, for `ε ∈ εs`.
+
+Use `transit_return(exits, entries)` to transform the result to transit and return time
+statistics instead.
+
+## Keywords
+* `diffeq...`: All keywords are propagated to the integrators of DifferentialEquations.jl
+  for continuous systems (discrete have no keywords).
+
+## Description
+Transit time statistics are important for the transport properties of dynamical systems[^Meiss1997]
+and can even be connected with the fractal dimension of chaotic sets[^Boev2014].
+
+The current algorithm collects exit and re-entry times to given sets in the state space,
+which are centered at `u₀` (**algorithm always starts at `u₀`** and the initial state of
+`ds` is irrelevant). `εs` is always a `Vector`.
+
+The sets around `u₀` are either nested hyper-spheres of radius `ε ∈ εs`, if each entry of
+`εs` is a real number. The sets can also be
+hyper-rectangles (boxes), if each entry of `εs` is a vector itself.
+Then, the `i`-th box is defined by the space covered by `u0 .± εs[i]` (thus the actual
+box size is `2εs[i]`!).
+
+The reason to input multiple `εs` at once is purely for performance.
+
+For discrete systems, exit time is recorded immediatelly after exitting of the set, and
+re-entry is recorded immediatelly on re-entry. This means that if an orbit needs
+1 step to leave the set and then it re-enters immediatelly on the next step, the return time
+is 1. For continuous systems high-order
+interpolation is done to accurately record the time of exactly crossing the `ε`-ball/box.
+
+[^Meiss1997]: Meiss, J. D. *Average exit time for volume-preserving maps*, Chaos (1997)](https://doi.org/10.1063/1.166245)
+
+[^Boev2014]: Boev, Vadivasova, & Anishchenko, *Poincaré recurrence statistics as an indicator of chaos synchronization*, Chaos (2014)](https://doi.org/10.1063/1.4873721)
+"""
+function transit_time_statistics(ds::DynamicalSystem, u0, εs, T; diffeq...)
+    check_εs_sorting(εs, length(u0))
+    integ = integrator(ds, u0; diffeq...)
+    transit_time_statistics(integ, u0, εs, T)
+end
+
+"""
+    transit_return(exits, entries) → transit, return
+Convert the output of [`return_time_statistics`](@ref) to the transit and return times.
+"""
+function transit_return(exits, entries)
+    returns = [en .- view(ex, 1:length(en)) for (en, ex) in zip(entries, exits)]
+    transits = similar(entries)
+    for (j, (en, ex)) in enumerate(zip(entries, exits))
+        M, N = length(en), length(ex)
+        enr, exr = M == N ? (1:N-1, 2:N) : (1:M, 2:N)
+        transits[j] = view(ex, exr) .- view(en, enr)
+    end
+    return transits, returns
+end
+
+##########################################################################################
+# ε-distances
+##########################################################################################
+function check_εs_sorting(εs, L)
+    correct = if εs[1] isa Real
+        issorted(εs; rev = true)
+    elseif εs[1] isa AbstractVector
+        @assert all(e -> length(e) == L, εs) "Boxes must have same dimension as state space!"
+        issorted(εs; by = maximum, rev = true)
+    end
+    if !correct
+        throw(ArgumentError("`εs` must be sorted from largest to smallest ball/box size."))
+    end
+    return correct
+end
+
+# Support both types of sets: balls and boxes
+"Return `true` if state is outside ε-ball"
+function isoutside(u, u0, ε::AbstractVector)
+    @inbounds for i in 1:length(u)
+        abs(u[i] - u0[i]) > ε[i] && return true
+    end
+    return false
+end
+isoutside(u, u0, ε::Real) = euclidean(u, u0) > ε
+
+"Return the **signed** distance of state to ε-ball (negative means inside ball)"
+function εdistance(u, u0, ε::AbstractVector)
+    m = eltype(u)(-Inf)
+    @inbounds for i in 1:length(u)
+        m2 = abs(u[i] - u0[i]) - ε[i]
+        m2 > m && (m = m2)
+    end
+    return m
+end
+εdistance(u, u0, ε::Real) = euclidean(u, u0) - ε
+
+##########################################################################################
+# Discrete systems
+##########################################################################################
+function transit_time_statistics(integ::MDI, u0, εs, T)
+    E = length(εs); check_εs_sorting(εs, length(u0))
+    pre_outside = fill(false, length(εs)) # `true` if outside the ball. Previous step
+    cur_outside = copy(pre_outside)       # current step.
+    exits = [typeof(integ.t)[] for _ in 1:length(εs)]
+    entries = [typeof(integ.t)[] for _ in 1:length(εs)]
+
+    while integ.t < T
+        step!(integ)
+
+        # here i gives the index of the largest ε-ball that the trajectory is out of.
+        # It is guaranteed that the trajectory is thus outside all other boxes
+        i = first_outside_index(integ, u0, εs, E) # TODO: Continuous version
+        cur_outside[i:end] .= true
+        cur_outside[1:i-1] .= false
+
+        update_exit_times!(exits, i, pre_outside, cur_outside, integ)
+        update_entry_times!(entries, i, pre_outside, cur_outside, integ)
+        pre_outside .= cur_outside
+    end
+    return_times = [en .- view(ex, 1:length(en)) for (en, ex) in zip(entries, exits)]
+    return exits, entries
+end
+
+function first_outside_index(integ::MDI, u0, εs, E)::Int
+    i = findfirst(e -> isoutside(integ.u, u0, e), εs)
+    return isnothing(i) ? E+1 : i
+end
+
+function update_exit_times!(exits, i, pre_outside, cur_outside, integ::MDI)
+    @inbounds for j in i:length(pre_outside)
+        cur_outside[j] && !pre_outside[j] && push!(exits[j], integ.t)
+    end
+end
+
+function update_entry_times!(entries, i, pre_outside, cur_outside, integ::MDI)
+    # TODO: Can I use `i` here?
+    @inbounds for j in 1:length(pre_outside)
+        pre_outside[j] && !cur_outside[j] && push!(entries[j], integ.t)
+    end
+end
+
+
+##########################################################################################
+# Continuous
+##########################################################################################

test/period_return/transit_time_tests.jl

@@ -0,0 +1,106 @@
+using ChaosTools, Test, Statistics
+
+println("\nTesting transit time statistics...")
+@testset "Transit time statistics" begin
+
+@testset "Standard map (exact)" begin
+    # INPUT
+    ds = Systems.standardmap()
+    T = 10000 # maximum time
+
+    # period 3 of standard map
+    u0 = SVector(0.8121, 1.6243)
+    εs = sort!([4.0, 2.0, 0.01]; rev=true)
+    # radius of 4.0 covers the first 2 points of the period 3 while 2.0 covers only
+    # the first point. Therefore:
+    # all first return times must be 1, and all second return times must be 2,
+    # and 3rd must be same as second (because it doesn't matter how close are)
+    exits, entries = transit_time_statistics(ds, u0, εs, T)
+    transits, returns = transit_return(exits, entries)
+
+    @test all(x -> length(x) > 5, exits)
+    @test all(x -> length(x) > 5, entries)
+    @test all(issorted, exits)
+    @test all(issorted, entries)
+    @test all(isequal(2), transits[1])
+    @test all(isequal(1), transits[2])
+    @test transits[2] == transits[3]
+    @test all(isequal(1), returns[1])
+    @test all(isequal(2), returns[2])
+    @test returns[2] == returns[3]
+
+    # quasiperiodic around period 3:
+    u0 = SVector(0.877, 1.565)
+    εs = sort!([4.0, 0.5, 0.1]; rev=true)
+    exits, entries = transit_time_statistics(ds, u0, εs, T)
+    transits, returns = transit_return(exits, entries)
+
+    @test all(issorted, exits)
+    @test all(issorted, entries)
+    @test all(x -> length(x) > 5, exits)
+    @test all(x -> length(x) > 5, entries)
+
+    # For ε=4.0, nothing changes with the before
+    @test all(isequal(1), returns[1])
+    @test all(isequal(2), transits[1])
+
+    # Similarly, 0.5 should be the same as before
+    @test all(isequal(1), transits[2])
+    @test all(isequal(2), returns[2])
+
+    # But now, the third entry is different, because it has the size of the quasiperiodic
+    # stability island torous
+    @test returns[3] ≠ returns[2]
+    @test transits[3] ≠ transits[2]
+    @test any(>(3), returns[3])
+    @test all(isequal(1), transits[3]) # still need only one step to exit
+end
+
+@testset "Towel map (boxes)" begin
+    to = Systems.towel()
+    tr = trajectory(to, 5000; Ttr = 10)
+    u0 = tr[3000]
+
+    # With these boxes, in the first 5 steps, the trajectory enters the y and z range
+    # but not the x range. Therefore it should NOT enter the box. See figure!
+    εs = [
+        SVector(0.05, 0.05, 0.125),
+        SVector(0.005, 0.005, 0.025),
+    ]
+
+    # Visual guidance
+    # using PyPlot
+    # tr0 = trajectory(to, 5, u0)
+    # fig, axs = subplots(1,3)
+    # comb = ((1, 2), (1, 3), (2, 3))
+    # for i in 1:3
+    #     j, k = comb[i]
+    #     ax = axs[i]
+    #     ax.scatter(tr[:, j], tr[:, k], s = 2, color = "C$(i-1)")
+    #     ax.scatter([u0[j]], [u0[k]], s = 20, color = "k")
+    #     ax.plot(tr0[:, j], tr0[:, k], color = "k")
+    #     for l in 1:length(εs)
+    #         rect = matplotlib.patches.Rectangle(
+    #         u0[[j, k]] .- εs[l][[j, k]], 2εs[l][j], 2εs[l][k],
+    #         alpha = 0.25, color = "k"
+    #         )
+    #         ax.add_artist(rect)
+    #     end
+    # end
+
+    exits, entries = transit_time_statistics(to, u0, εs, 10000)
+    transits, returns = transit_return(exits, entries)
+
+    @test all(issorted, exits)
+    @test all(issorted, entries)
+    @test length(exits[1]) > length(exits[2])
+    @test returns[1][1] > 5
+    @test mean(returns[1]) < mean(returns[2])
+end
+# ds = Systems.roessler()
+# T = 1000.0 # maximum time
+# u0 = trajectory(ds, 10; Ttr = 100)[end] # return center
+# εs = sort!([0.1, 0.01, 0.001]; rev=true)
+# diffeq = (;)
+
+end

test/runtests.jl

@@ -13,6 +13,7 @@ include("chaosdetection/expansionentropy_tests.jl")
 
 include("period_return/periodicity_tests.jl")
 include("period_return/period_tests.jl")
+include("period_return/transit_time_tests.jl")
 
 include("dimensions/entropy_dimension.jl")
 include("nlts_tests.jl")