New tests for pecuzal, uzal_cost; corrected normalization in uzal cost (

#97)

* fixed tests w.r.t. deprecations in Peaks

* corrected normalization factor; revised tests; included time series

* transfer time series into JuliaDynamics-repo and download them when testing

* final revisions

* for loop in time series download

src/unified_de/pecora.jl

@@ -206,9 +206,7 @@ function pecora(
         p::T = 0.5) where {D, T<:Real}
 
     @assert K ≥ 8 "You must provide a δ-neighborhood size consisting of at least 8 neighbors."
-    @assert all(x -> x ≥ 0, τs) "τ's and j's for generalized embedding must be positive integers"
-    @assert all(x -> x ≥ 0, js) "τ's and j's for generalized embedding must be positive integers"
-    @assert all(x -> x ≥ 0, delays) "considered delay values must be positive integers"
+    @assert all(x -> x ≥ 0, js) "j's for generalized embedding must be positive integers"
     @assert 0 < samplesize ≤ 1 "`samplesize` must be ∈ (0,1]"
 
     N = floor(Int,samplesize*length(s)) #number of fiducial points

src/unified_de/pecuzal.jl

@@ -545,14 +545,14 @@ function compute_L_decrease(E²::Array{P, 2}, E²_trial::Array{P, 2}, ϵ²::Vect
     E²_avrg = mean(E²[1:NN,1:T], dims=2)                   # Eq. 15
     σ² = E²_avrg ./ ϵ²[1:NN] # noise amplification σ², Eq. 17
     σ²_avrg = mean(σ²) # averaged value of the noise amplification, Eq. 18
-    α² = 1 / sum(ϵ²[1:NN].^(-1)) # for normalization, Eq. 21
+    α² = 1 / mean(ϵ²[1:NN].^(-1)) # for normalization, Eq. 21
     L = log10(sqrt(σ²_avrg)*sqrt(α²))
     # 2nd dataset
     # Average E²[T] over all prediction horizons
     E²_avrg_trial = mean(E²_trial[1:NN,1:T], dims=2)                   # Eq. 15
     σ²_trial = E²_avrg_trial ./ ϵ²_trial[1:NN] # noise amplification σ², Eq. 17
     σ²_avrg_trial = mean(σ²_trial) # averaged value of the noise amplification, Eq. 18
-    α²_trial = 1 / sum(ϵ²_trial[1:NN].^(-1)) # for normalization, Eq. 21
+    α²_trial = 1 / mean(ϵ²_trial[1:NN].^(-1)) # for normalization, Eq. 21
     L_trial = log10(sqrt(σ²_avrg_trial)*sqrt(α²_trial))
 
     return L_trial - L

src/unified_de/uzal_cost.jl

@@ -93,7 +93,7 @@ function uzal_cost(Y::AbstractDataset{D, ET};
     end
     σ² = E²_avrg ./ ϵ² # noise amplification σ², Eq. 17
     σ²_avrg = mean(σ²) # averaged value of the noise amplification, Eq. 18
-    α² = 1 / sum(ϵ².^(-1)) # for normalization, Eq. 21
+    α² = 1 / mean(ϵ².^(-1)) # for normalization, Eq. 21
     L = log10(sqrt(σ²_avrg)*sqrt(α²))
 end
 
@@ -184,7 +184,7 @@ function uzal_cost_local(Y::AbstractDataset{D, ET};
         E²_avrg[i] = mean(E²)                   # Eq. 15
     end
     σ² = E²_avrg ./ ϵ² # noise amplification σ², Eq. 17
-    α² = 1 / sum(ϵ².^(-1)) # for normalization, Eq. 21
+    α² = 1 / mean(ϵ².^(-1)) # for normalization, Eq. 21
     L_local = log10.(sqrt.(σ²).*sqrt(α²))
     return L_local
 end

test/runtests.jl

@@ -3,6 +3,8 @@ using StaticArrays
 using Test
 
 # Download some test timeseries
+
+# Old:
 repo = "https://raw.githubusercontent.com/JuliaDynamics/NonlinearDynamicsTextbook/master/exercise_data"
 tsfolder = joinpath(@__DIR__, "timeseries")
 todownload = ["$n.csv" for n in 1:4]
@@ -12,6 +14,13 @@ for a in todownload
     download(repo*"/"*a, joinpath(tsfolder, a))
 end
 
+#New:
+todownload = ["test_time_series_lorenz_standard_N_10000_multivariate.csv", "test_time_series_roessler_N_10000_multivariate.csv"]
+repo = "https://raw.githubusercontent.com/JuliaDynamics/JuliaDynamics/master/timeseries"
+for a in todownload
+    download(repo*"/"*a, joinpath(tsfolder, a))
+end
+
 ti = time()
 
 diffeq = (atol = 1e-9, rtol = 1e-9, maxiters = typemax(Int))
@@ -23,14 +32,15 @@ diffeq = (atol = 1e-9, rtol = 1e-9, maxiters = typemax(Int))
     include("traditional/embedding_dimension_test.jl")
     include("unified/test_pecora.jl")
     include("unified/uzal_cost_test.jl")
+    include("unified/test_pecuzal_embedding.jl")
 
     # The following methods have been tested throughly and also published in research.
-    # We know that they work, but unfortunately the tests we have written 
-    # about them are not good. 
+    # We know that they work, but unfortunately the tests we have written
+    # about them are not good.
     # See https://github.com/JuliaDynamics/DelayEmbeddings.jl/issues/95
     # include("unified/mdop_tests.jl")
     # include("unified/test_garcia.jl")
-    # include("unified/test_pecuzal_embedding.jl")
+
 end
 
 ti = time() - ti

test/unified/test_pecuzal_embedding.jl

@@ -1,11 +1,20 @@
-using DelayEmbeddings, DynamicalSystemsBase
-using Test, Random
+using DelayEmbeddings
+using Test, Random, DelimitedFiles
 
 println("\nTesting pecuzal_method.jl...")
 @testset "PECUZAL" begin
 
-lo = Systems.lorenz([1.0, 1.0, 50.0])
-tr = trajectory(lo, 100; Δt = 0.01, Ttr = 10)
+## Test Lorenz example
+# For comparison reasons using Travis CI we carry out the integration on a UNIX
+# OS and save the resulting time series
+# See https://github.com/JuliaDynamics/JuliaDynamics for the storage of
+# the time series used for testing
+#
+# u0 = [0, 10.0, 0.0]
+# lo = Systems.lorenz(u0; σ=10, ρ=28, β=8/3)
+# tr = trajectory(lo, 100; Δt = 0.01, Ttr = 100)
+tr = readdlm(joinpath(tsfolder, "test_time_series_lorenz_standard_N_10000_multivariate.csv"))
+tr = Dataset(tr)
 
 @testset "Univariate example" begin
 
@@ -14,19 +23,21 @@ tr = trajectory(lo, 100; Δt = 0.01, Ttr = 10)
     Tmax = 100
 
     @time Y, τ_vals, ts_vals, Ls , εs = pecuzal_embedding(s[1:5000];
-                                        τs = 0:Tmax , w = w)
+                                        τs = 0:Tmax , w = w, L_threshold = 0.05, econ=true)
 
-    @test -0.728 < Ls[1] < -0.727
-    @test -0.324 < Ls[2] < -0.323
+    @test -1.08 < Ls[1] < -1.04
+    @test -0.51 < Ls[2] < -0.47
+    @test -0.18 < Ls[3] < -0.14
 
-    @test τ_vals[2] == 18
-    @test τ_vals[3] == 9
-    @test length(ts_vals) == 3
+    @test τ_vals[2] == 19
+    @test τ_vals[3] == 96
+    @test τ_vals[4] == 80
+    @test length(ts_vals) == 4
 
-    @time Y, τ_vals, ts_vals, Ls , εs = pecuzal_embedding(s[1:5000];
-                                        τs = 0:Tmax , w = w, econ = true)
-    @test -0.728 < Ls[1] < -0.726
-    @test -0.322 < Ls[2] < -0.321
+    @time Y, τ_vals, ts_vals, Ls , εs = pecuzal_embedding(s;
+                                        τs = 0:Tmax , w = w, L_threshold = 0.05)
+    @test -0.95 < Ls[1] < -0.91
+    @test -0.48 < Ls[2] < -0.44
 
     @test τ_vals[2] == 18
     @test τ_vals[3] == 9
@@ -41,7 +52,7 @@ end
 
 
 @testset "Multivariate example" begin
-## Test of the proposed Pecora-Uzal-embedding-method (multivariate case)
+    ## Test of the proposed PECUZAL-embedding-method (multivariate case)
 
     w1 = estimate_delay(tr[:,1], "mi_min")
     w2 = estimate_delay(tr[:,2], "mi_min")
@@ -52,27 +63,23 @@ end
     @time Y, τ_vals, ts_vals, Ls , ε★ = pecuzal_embedding(tr[1:5000,:];
                                         τs = 0:Tmax , w = w, econ = true)
 
-    @test length(ts_vals) == 5
-    @test ts_vals[3] == ts_vals[4] == ts_vals[5] == 1
+    @test length(ts_vals) == 4
+    @test ts_vals[2] == ts_vals[3] == ts_vals[4] == 1
     @test ts_vals[1] == 3
-    @test ts_vals[2] == 2
     @test τ_vals[1] == 0
-    @test τ_vals[2] == 0
-    @test τ_vals[3] == 62
-    @test τ_vals[4] == 48
-    @test τ_vals[5] == 0
-    @test -0.9338 < Ls[1] < -0.9337
-    @test -0.356 < Ls[2] < -0.355
-    @test -0.1279 < Ls[3] < -0.1278
-    @test -0.015 < Ls[4] < -0.014
+    @test τ_vals[2] == 9
+    @test τ_vals[3] == 64
+    @test τ_vals[4] == 53
+    @test -1.40 < Ls[1] < -1.36
+    @test -0.76 < Ls[2] < -0.72
+    @test -0.1 < Ls[3] < -0.06
 
     @time Y, τ_vals, ts_vals, Ls , ε★ = pecuzal_embedding(tr[1:5000,:];
                                         τs = 0:Tmax , w = w, econ = true, L_threshold = 0.2)
-    @test -0.9338 < Ls[1] < -0.9337
-    @test -0.356 < Ls[2] < -0.355
+    @test -1.40 < Ls[1] < -1.36
+    @test -0.76 < Ls[2] < -0.72
     @test size(Y,2) == 3
 
-
 end
 
 @testset "Dummy example" begin