From afb795240f46c86497f87592a4b8443a448728aa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kai=20Hauke=20Kr=C3=A4mer?= <36292849+hkraemer@users.noreply.github.com> Date: Tue, 31 Aug 2021 14:57:18 +0200 Subject: [PATCH] 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 | 4 +- src/unified_de/pecuzal.jl | 4 +- src/unified_de/uzal_cost.jl | 4 +- test/runtests.jl | 16 ++- test/unified/test_pecuzal_embedding.jl | 65 +++++---- test/unified/uzal_cost_test.jl | 192 +++++++++---------------- 6 files changed, 124 insertions(+), 161 deletions(-) diff --git a/src/unified_de/pecora.jl b/src/unified_de/pecora.jl index 51a2906..ad25446 100644 --- a/src/unified_de/pecora.jl +++ b/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 diff --git a/src/unified_de/pecuzal.jl b/src/unified_de/pecuzal.jl index d5037fb..80a170a 100644 --- a/src/unified_de/pecuzal.jl +++ b/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 diff --git a/src/unified_de/uzal_cost.jl b/src/unified_de/uzal_cost.jl index 9e9a033..303903b 100644 --- a/src/unified_de/uzal_cost.jl +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index 32263e5..e5f2b91 100644 --- a/test/runtests.jl +++ b/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 diff --git a/test/unified/test_pecuzal_embedding.jl b/test/unified/test_pecuzal_embedding.jl index 9619e07..582751a 100644 --- a/test/unified/test_pecuzal_embedding.jl +++ b/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 diff --git a/test/unified/uzal_cost_test.jl b/test/unified/uzal_cost_test.jl index 791e433..06fef89 100644 --- a/test/unified/uzal_cost_test.jl +++ b/test/unified/uzal_cost_test.jl @@ -1,4 +1,3 @@ -using DynamicalSystemsBase using DelayEmbeddings using Test using Random @@ -9,15 +8,23 @@ println("\nTesting uzal_cost.jl...") @testset "Random vectors" begin ## Check on random vector Random.seed!(1516578735) - tr = randn(10000) - tr = Dataset(tr) - L = uzal_cost(tr; + tr1 = randn(10000) + tr2 = randn(5000) + tr1 = Dataset(tr1) + tr2 = Dataset(tr2) + L = uzal_cost(tr1; + Tw = 60, K= 3, w = 1, samplesize = 1.0, + metric = Euclidean() + ) + L2 = uzal_cost(tr2; Tw = 60, K= 3, w = 1, samplesize = 1.0, metric = Euclidean() ) - L_max = -2.059 - @test L < L_max + L_max = -0.06 + L_min = -0.08 + @test L_min < L < L_max + @test L_min < L2 < L_max ## check on clean step function (0-distances) @@ -41,11 +48,20 @@ println("\nTesting uzal_cost.jl...") Tw = 60, K= 3, w = 1, samplesize = 1.0, metric = Euclidean() ) - @test L<-3 + @test -1.6 < L < -1.5 end @testset "Uzal local cost (Lorenz)" begin - tr = readdlm(joinpath(tsfolder, "3.csv")) + ## 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) # check local cost function output @@ -58,25 +74,29 @@ end @test length(L_local) == length(tr)-Tw @test maximum(L_local)>L @test minimum(L_local)