Merge pull request #86 from JuliaDynamics/partially

re-open PR #83

.gitignore

@@ -2,3 +2,4 @@
 *.jl.*.cov
 *.jl.mem
 docs/build
+Manifest.toml

CHANGELOG.md

@@ -1,3 +1,7 @@
+# 1.6
+* Implementation of  Wernecke, H., Sándor, B. & Gros, C.
+      *How to test for partially predictable chaos*. [Scientific Reports **7**, (2017)](https://www.nature.com/articles/s41598-017-01083-x). Thanks and welcome to our new contributor @efosong ! Implemented with the function `predictability`.
+
 # 1.5
 * Updated everything to new DiffEqBase 5.0 and the new default integrator (`SimpleATsit5`)
 * Simplified low-level call signature for `poincaresos`.

Project.toml

@@ -1,12 +1,13 @@
 name = "ChaosTools"
 uuid = "608a59af-f2a3-5ad4-90b4-758bdf3122a7"
 repo = "https://github.com/JuliaDynamics/ChaosTools.jl.git"
-version = "1.5.1"
+version = "1.6.0"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 DynamicalSystemsBase = "6e36e845-645a-534a-86f2-f5d4aa5a06b4"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/partially_predictable.jl

@@ -2,31 +2,88 @@ using LinearAlgebra, Distributions
 using Statistics: mean
 export predictability
 
-# TODO - write docstring
-#
-# [1] H. Wernecke, B. Sándor, and C. Gros, ‘How to test for partially
-#     predictable chaos’, Scientific Reports, vol. 7, no. 1, Dec. 2017.
+"""
+   predictability(ds::DynamicalSystem; kwargs...) -> chaos_type, ν, C
 
+Determine whether `ds` displays strongly chaotic, partially-predictable chaotic
+or regular behaviour, using the method by Wernecke et al. described in [1].
+
+Return the type of the behavior, the cross-distance scaling coefficient `ν`
+and the correlation coefficient `C`.
+Typical values for `ν`, `C` and `chaos_type` are given in Table 2 of [1]:
+
+| `chaos_type` | `ν` | `C` |
+|--------------|-----|-----|
+|    `:SC`     |  0  |  0  |
+|    `:PPC`    |  0  |  1  |
+|    `:REG`    |  1  |  1  |
+
+
+## Keyword Arguments
+* `Ttr = 200` : Extra "transient" time to evolve the system before sampling from
+   the trajectory. Should be `Int` for discrete systems.
+* `T_sample = 1e4` : Time to evolve the system for taking samples. Should be
+  `Int` for discrete systems.
+* `n_samples = 500` : Number of samples to take for use in calculating statistics.
+* `λ_max = lyapunov(ds, 5000)` : Value to use for largest Lyapunov exponent
+  for finding the Lyapunov prediction time. If it is less than zero a regular
+  result is returned immediatelly.
+* `d_tol = 1e-3` : tolerance distance to use for calculating Lyapunov prediction time.
+* `T_multiplier = 10` : Multiplier from the Lyapunov prediction time to the evaluation time.
+* `T_max = Inf` : Maximum time at which to evaluate trajectory distance. If the internally
+   computed evaluation time is larger than `T_max`, stop at `T_max` instead.
+* `δ_range = 10.0 .^ (-9:-6)` : Range of initial condition perturbation distances
+   to use to determine scaling `ν`.
+* `diffeq...` : Keyword arguments propagated into `init` of DifferentialEquations.jl.
+  See [`trajectory`](@ref) for examples. Only valid for continuous systems.
+
+## Description
+Samples points from a trajectory of the system to be used as initial conditions. Each of
+these initial conditions is randomly perturbed by a distance `δ`, and the trajectories for
+both the original and perturbed initial conditions are computed to the 'evaluation time'
+`T`.
+
+The average (over the samples) distance and cross-correlation coefficient
+of the state at time `T` is
+computed. This is repeated for a range of `δ` (defined by `δ_range`), and linear
+regression is used to determine how the distance and cross-correlation scale with `δ`,
+allowing for identification of chaos type.
+
+The evaluation time `T` is calculated as `T = T_multiplier*Tλ`, where the Lyapunov
+prediction time `Tλ = log(d_tol/δ)/λ_max`. This may be very large if the `λ_max` is small,
+e.g. when the system is regular, so this internally computed time `T` can be overridden by
+a smaller `T_max` set by the user.
+
+## Performance Notes
+For continuous systems, it is likely that the `maxiters` used by the integrators needs to
+be increased, e.g. to 1e9. This is part of the `diffeq` kwargs.
+In addition, be aware that this function does a *lot* of internal computations.
+It is operating in a different speed than e.g. [`lyapunov`](@ref).
+
+## References
+
+[1] : Wernecke, H., Sándor, B. & Gros, C.
+      *How to test for partially predictable chaos*. [Scientific Reports **7**, (2017)](https://www.nature.com/articles/s41598-017-01083-x).
+"""
 function predictability(ds::DynamicalSystem;
-                        T_transient::Real = 200,
-                        T_sample::Real = 1e5,
-                        n_samples::Integer = 1000,
-                        λ_max::Real = abs(lyapunov(ds, 5000)),
+                        Ttr::Real = 200,
+                        T_sample::Real = 1e4,
+                        n_samples::Integer = 500,
+                        λ_max::Real = lyapunov(ds, 5000),
                         d_tol::Real = 1e-3,
                         T_multiplier::Real = 10,
                         T_max::Real = Inf,
-                        δ_range::AbstractArray{T,1} = 10.0 .^ (-9:-6),
+                        δ_range::AbstractArray = 10.0 .^ (-9:-6),
                         diffeq...
-                       ) where T <: Real
+                        )
 
-    # ======================== Internal Constants ========================= #
+    λ_max < 0 && return :REG, 1.0, 1.0
+    # Internal Constants
     ν_threshold = 0.5
     C_threshold = 0.5
-    # ===================================================================== #
-
 
     # Sample points from a single trajectory of the system
-    samples = sample_trajectory(ds, T_transient, T_sample, n_samples; diffeq...)
+    samples = sample_trajectory(ds, Ttr, T_sample, n_samples; diffeq...)
 
     # Calculate the mean position and variance of the trajectory. ([1] pg. 5)
     # Using samples 'Monte Carlo' approach instead of direct integration
@@ -38,6 +95,8 @@ function predictability(ds::DynamicalSystem;
     correlations = Float64[] # Cross-correlation at time T for different δ
     p_integ = parallel_integrator(ds, samples[1:2]; diffeq...)
     for δ in δ_range
+        # some kind of warning should be thrown for very large
+        # Tλ.
         Tλ = log(d_tol/δ)/λ_max
         T = min(T_multiplier * Tλ, T_max)
         Σd = 0
@@ -71,13 +130,13 @@ function predictability(ds::DynamicalSystem;
     # Perform regression to check cross-distance scaling
     ν = slope(log.(δ_range), log.(distances))
     C = mean(correlations)
-    
+
     # Determine chaotic nature of the system
     if ν > ν_threshold && C > C_threshold
-        chaos_type = :LAM
+        chaos_type = :REG
     elseif ν <= ν_threshold && C > C_threshold
         chaos_type = :PPC
-    elseif ν <= ν_threshold && C <= C_threshold
+    elseif ν <= ν_threshold && C ≤ C_threshold
         chaos_type = :SC
     else
         # Covers the case when ν > ν_threshold but C <= C_threshold
@@ -88,42 +147,42 @@ function predictability(ds::DynamicalSystem;
 end
 
 
-function sample_trajectory(ds::ContinuousDynamicalSystem, 
-                           T_transient::Real, T_sample::Real, 
-                           n_samples::Real; 
+function sample_trajectory(ds::ContinuousDynamicalSystem,
+                           Ttr::Real, T_sample::Real,
+                           n_samples::Real;
                            diffeq...)
     # Samples *approximately* `n_samples` points.
     β = T_sample/n_samples
     D_sample = Exponential(β)
-    sample_trajectory(ds, T_transient, T_sample, D_sample; diffeq...)
+    sample_trajectory(ds, Ttr, T_sample, D_sample; diffeq...)
 end
 
-function sample_trajectory(ds::DiscreteDynamicalSystem, 
-                           T_transient::Real, T_sample::Real, 
+function sample_trajectory(ds::DiscreteDynamicalSystem,
+                           Ttr::Real, T_sample::Real,
                            n_samples::Real;
                            diffeq...)
-    @assert n_samples < T_sample "DiscreteDynamicalSystems must satisfy n_samples < T_sample"
+    @assert n_samples < T_sample "discrete systems must satisfy n_samples < T_sample"
     # Samples *approximately* `n_samples` points.
     p = n_samples/T_sample
     D_sample = Geometric(p)
-    sample_trajectory(ds, T_transient, T_sample, D_sample; diffeq...)
+    sample_trajectory(ds, Ttr, T_sample, D_sample; diffeq...)
 end
 
-function sample_trajectory(ds::DynamicalSystem, 
-                           T_transient::Real, T_sample::Real,
+function sample_trajectory(ds::DynamicalSystem,
+                           Ttr::Real, T_sample::Real,
                            D_sample::UnivariateDistribution;
                            diffeq...)
     # Simulate initial transient
     integ = integrator(ds; diffeq...)
-    while integ.t < T_transient
+    while integ.t < Ttr
         step!(integ)
     end
 
     # Time to the next sample is sampled from the distribution D_sample
     # e.g. Continuous systems: D_sample is Exponential distribution
     samples = typeof(integ.u)[]
-    while integ.t < T_transient + T_sample
-        step!(integ, rand(D_sample), true)
+    while integ.t < Ttr + T_sample
+        step!(integ, rand(D_sample))
         push!(samples, integ.u)
     end
     samples

test/chaos_detection_tests.jl → test/gali_tests.jl

@@ -5,7 +5,7 @@ using OrdinaryDiffEq: Tsit5
 
 test_value = (val, vmin, vmax) -> @test vmin <= val <= vmax
 
-println("\nTesting chaos detection algorithms...")
+println("\nTesting GALI...")
 
 @testset "GALI discrete" begin
 

test/partially_predictable_tests.jl

@@ -3,10 +3,10 @@ using Test, StaticArrays
 using OrdinaryDiffEq: Vern9
 using DynamicalSystemsBase.Systems: lorenz
 
-ν_thresh_lower, ν_thresh_upper = 0.05, 0.95
-C_thresh_lower, C_thresh_upper = 0.05, 0.95
+ν_thresh_lower, ν_thresh_upper = 0.1, 0.9
+C_thresh_lower, C_thresh_upper = 0.1, 0.9
 
-println("\nTesting predictability...")
+println("\nTesting Partially predictable chaos...")
 
 @testset "Predictability Continuous" begin
     @testset "Lorenz map" begin
@@ -34,9 +34,23 @@ println("\nTesting predictability...")
         @testset "Lorenz map - laminar" begin
             lz = lorenz(ρ=181.10)
             chaos_type, ν, C = predictability(lz; T_max = 400, alg=Vern9(), maxiters=1e9)
-            @test chaos_type == :LAM
+            @test chaos_type == :REG
             @test ν > ν_thresh_upper
             @test C > C_thresh_upper
         end
     end
 end
+
+@testset "Predictability Discrete" begin
+    @testset "Henon map" begin
+        ds = Systems.henon()
+        a_reg = 0.8; a_ppc = 1.11; a_reg2 = 1.0; a_cha = 1.4
+        res = [:REG, :REG, :PPC, :SC]
+
+        for (i, a) in enumerate((a_reg, a_reg2, a_ppc, a_cha))
+            set_parameter!(ds, 1, a)
+            chaos_type, ν, C = predictability(ds; T_max = 400000)
+            @test chaos_type == res[i]
+        end
+    end
+end

test/runtests.jl

@@ -1,8 +1,11 @@
 ti = time()
 
-# Continuous
+# Continuous and discrete
 include("lyapunov_exponents.jl")
-include("chaos_detection_tests.jl")
+include("gali_tests.jl")
+include("partially_predictable_tests.jl")
+
+# Continuous
 include("poincare_tests.jl")
 
 # Discrete