update predictability to DynamicalSystems.jl v3.0. (#312)

* update predictability to DynamicalSystems.jl v3.0.

* update  and its tests to v3.0: add correct reinitialization and others

* update  and its tests to v3.0: enable tests

* update predictability to DS v3.0: minor bugfixes

Project.toml

@@ -1,7 +1,7 @@
 name = "ChaosTools"
 uuid = "608a59af-f2a3-5ad4-90b4-758bdf3122a7"
 repo = "https://github.com/JuliaDynamics/ChaosTools.jl.git"
-version = "3.1.0"
+version = "3.1.1"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"

src/chaosdetection/partially_predictable.jl

@@ -1,9 +1,12 @@
 using LinearAlgebra
 using Random: Xoshiro
-using Distributions: Exponential, Geometric, UnivariateDistribution
+using Distributions: Exponential, Geometric, Normal, UnivariateDistribution
 using Statistics: mean
 export predictability
 
+slope(x, y) = linreg(x, y)[2]
+
+
 """
     predictability(ds::CoreDynamicalSystem; kwargs...) -> chaos_type, ν, C
 
@@ -74,15 +77,13 @@ function predictability(ds::CoreDynamicalSystem;
         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 = 10.0 .^ (-9:-6),
-        ν_threshold = 0.5, C_threshold = 0.5,
+        ν_threshold = 0.5, C_threshold = 0.5
     )
 
-    error("This function has not yet been updated to DynamicalSystems.jl v3.0.
-    Please consider a Pull Request :) (very easy!)")
-
+	reinit!(ds)
     rng = Xoshiro()
     λ_max < 0 && return :REG, 1.0, 1.0
-    samples = sample_trajectory(ds, Ttr, T_sample, n_samples; diffeq)
+    samples = sample_trajectory(ds, Ttr, T_sample, n_samples)
     # Calculate initial mean position and variance of the trajectory. (eq. [1] pg. 5)
     # using random samples approach instead of direct integration
     μ = mean(samples)
@@ -91,9 +92,7 @@ function predictability(ds::CoreDynamicalSystem;
     # Calculate cross-distance scaling and correlation scaling
     distances = Float64[] # Mean distances at time T for different δ
     correlations = Float64[] # Cross-correlation at time T for different δ
-    # TODO: p_integ is type unstable, the rest of the code should be moved
-    # into a separate function. It's also good todo for more readable code...
-    p_integ = parallel_integrator(ds, samples[1:2]; diffeq)
+    p_sys = ParallelDynamicalSystem(ds, samples[1:2])
     for δ in δ_range
         # TODO: some kind of warning should be thrown for very large Tλ
         Tλ = log(d_tol/δ)/λ_max
@@ -106,11 +105,11 @@ function predictability(ds::CoreDynamicalSystem;
             n /= norm(n)
             û = u + δ*n
             # re-integrate to time T
-            reinit!(p_integ, [u, û])
-            step!(p_integ, T)
+            reinit!(p_sys, [u, û])
+            step!(p_sys, T)
             # Accumulate distance and square-distance
             # TODO: This integrator access is not using the API!
-            d = norm(p_integ.u[1]-p_integ.u[2], 2)
+            d = norm(current_states(p_sys)[1]-current_states(p_sys)[2], 2)
             Σd  += d
             Σd² += d^2
         end
@@ -144,36 +143,33 @@ end
 # Samples *approximately* `n_samples` points via exponential distribution sampling times
 function sample_trajectory(ds::ContinuousDynamicalSystem,
                            Ttr::Real, T_sample::Real,
-                           n_samples::Real;
-                           diffeq = NamedTuple())
+                           n_samples::Real)
     β = T_sample/n_samples
     D_sample = Exponential(β)
-    sample_trajectory(ds, Ttr, T_sample, D_sample; diffeq)
+    sample_trajectory(ds, Ttr, T_sample, D_sample)
 end
 # Same as above but geometric distribution due to discrete nature of time
 function sample_trajectory(ds::DiscreteDynamicalSystem,
                            Ttr::Real, T_sample::Real,
-                           n_samples::Real;
-                           diffeq = NamedTuple())
+                           n_samples::Real) 
+
     @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, Ttr, T_sample, D_sample; diffeq)
+    sample_trajectory(ds, Ttr, T_sample, D_sample)
 end
 
 function sample_trajectory(ds::DynamicalSystem,
                            Ttr::Real, T_sample::Real,
-                           D_sample::UnivariateDistribution;
-                           diffeq = NamedTuple())
+                           D_sample::UnivariateDistribution)
     # Simulate initial transient
-    integ = integrator(ds; diffeq)
-    step!(integ, Ttr)
+    step!(ds, Ttr)
     # Time to the next sample is sampled from the distribution D_sample
-    samples = typeof(integ.u)[]
-    while integ.t < Ttr + T_sample
-        step!(integ, rand(D_sample))
-        push!(samples, integ.u)
+    samples = typeof(current_state(ds))[]
+    while current_time(ds) < Ttr + T_sample
+        step!(ds, rand(D_sample))
+        push!(samples, deepcopy(current_state(ds)))
     end
     samples
 end

test/chaosdetection/partially_predictable.jl

@@ -2,17 +2,30 @@ using ChaosTools
 using Test
 using Random
 
+@inbounds function lorenz_rule!(du, u, p, t)
+    σ = p[1]; ρ = p[2]; β = p[3]
+    du[1] = σ*(u[2]-u[1])
+    du[2] = u[1]*(ρ-u[3]) - u[2]
+    du[3] = u[1]*u[2] - β*u[3]
+    return nothing
+end
+
+u0 = [0, 10.0, 0]
+p0 = [10, 28, 8/3]
+diffeq = (maxiters = 1e9,abstol = 1.0e-6, reltol = 1.0e-6)
+
+lz = CoupledODEs(lorenz_rule!, u0, p0; diffeq)
+
 ν_thresh_lower, ν_thresh_upper = 0.1, 0.9
 C_thresh_lower, C_thresh_upper = 0.15, 0.9
-diffeqmaxit = (maxiters = 1e9,)
 println("\nTesting Partially predictable chaos...")
 
 @testset "Predictability Lorenz" begin
     @testset "strongly chaotic" begin
         Random.seed!(12)
-        lz = Systems.lorenz(ρ=180.70)
+        set_parameter!(lz,2,180.70)
         @time chaos_type, ν, C = predictability(lz;
-            λ_max = 1.22, diffeq=diffeqmaxit, T_max = 1e3, n_samples = 100,
+            λ_max = 1.22, T_max = 1e3, n_samples = 100,
         )
         @test chaos_type == :SC
         @test ν < ν_thresh_lower
@@ -21,9 +34,9 @@ println("\nTesting Partially predictable chaos...")
     end
     @testset "PPC 1" begin
         Random.seed!(12)
-        lz = Systems.lorenz(ρ=180.78)
+        set_parameter!(lz,2,180.78)
         @time chaos_type, ν, C = predictability(lz;
-            λ_max = 0.4, diffeq=diffeqmaxit, n_samples = 100, T_max = 400
+            λ_max = 0.4, n_samples = 100, T_max = 400
         )
         @test chaos_type == :PPC
         @test ν < ν_thresh_lower
@@ -32,9 +45,9 @@ println("\nTesting Partially predictable chaos...")
     end
     @testset "PPC 2" begin
         Random.seed!(12)
-        lz = Systems.lorenz(ρ=180.95)
+        set_parameter!(lz,2,180.95)
         @time chaos_type, ν, C = predictability(lz;
-            λ_max = 0.1, diffeq=diffeqmaxit, n_samples = 100, T_max = 400
+            λ_max = 0.1, n_samples = 100, T_max = 400
         )
         @test chaos_type == :PPC
         @test ν < ν_thresh_lower
@@ -43,10 +56,9 @@ println("\nTesting Partially predictable chaos...")
     end
     @testset "laminar" begin
         Random.seed!(12)
-        lz = Systems.lorenz(ρ=181.10)
+        set_parameter!(lz,2,181.10)
         @time chaos_type, ν, C = predictability(lz;
-            λ_max = 0.01, T_max = 400, n_samples = 100, diffeq=diffeqmaxit
-        )
+            λ_max = 0.01, T_max = 400, n_samples = 100)
         @test chaos_type == :REG
         @test ν > ν_thresh_upper
         @test C > C_thresh_upper
@@ -56,13 +68,17 @@ end
 
 @testset "Predictability Discrete" begin
     @testset "Henon map" begin
-        ds = Systems.henon()
+        henon_rule(x, p, n) = SVector(1.0 - p[1]*x[1]^2 + x[2], p[2]*x[1])
+		u0 = zeros(2)
+		p0 = [1.4, 0.3]
+		hn = DeterministicIteratedMap(henon_rule, u0, p0)
+
         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 = 100000)
+            set_parameter!(hn, 1, a)
+            chaos_type, ν, C = predictability(hn; T_max = 100000)
             @test chaos_type == res[i]
         end
     end

test/runtests.jl

@@ -22,7 +22,7 @@ include("timeevolution/orbitdiagram.jl")
 
 testfile("chaosdetection/lyapunovs.jl")
 testfile("chaosdetection/gali.jl")
-# include("chaosdetection/partially_predictable_tests.jl")
+include("chaosdetection/partially_predictable.jl")
 include("chaosdetection/01test.jl")
 testfile("chaosdetection/expansionentropy.jl")