Merge pull request #62 from JuliaDynamics/deprecate_reconstruct

Deprecate `reconstruct` in favor of `embed`.

CHANGELOG.md

@@ -1,3 +1,6 @@
+# v1.13.0
+* `reconstruct` is deprecated in favor of `embed`.
+
 # v1.12.0
 * Possible delay times in `optimal_traditional_de` are now `1:100` for increased accuracy.
 * New method for univariate non-unified delay embedding by Hegger, Kantz

Project.toml

@@ -1,7 +1,7 @@
 name = "DelayEmbeddings"
 uuid = "5732040d-69e3-5649-938a-b6b4f237613f"
 repo = "https://github.com/JuliaDynamics/DelayEmbeddings.jl.git"
-version = "1.12.1"
+version = "1.13.0"
 
 [deps]
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"

src/embeddings.jl

@@ -10,119 +10,97 @@ export WeightedDelayEmbedding, AbstractEmbedding
 abstract type AbstractEmbedding end
 
 """
-    DelayEmbedding(γ, τ) → `embedding`
-Return a delay coordinates embedding structure to be used as a functor,
+    DelayEmbedding(γ, τ, h = nothing) → `embedding`
+Return a delay coordinates embedding structure to be used as a function-like-object,
 given a timeseries and some index. Calling
 ```julia
 embedding(s, n)
 ```
 will create the `n`-th delay vector of the embedded space, which has `γ`
-temporal neighbors with delay(s) `τ`. See [`reconstruct`](@ref) for more.
+temporal neighbors with delay(s) `τ`.
+`γ` is the embedding dimension minus 1, `τ` is the delay time(s) while `h` are
+extra weights, as in [`embed`](@ref) for more.
 
 **Be very careful when choosing `n`, because `@inbounds` is used internally.**
 Use [`τrange`](@ref)!
 """
-struct DelayEmbedding{γ} <: AbstractEmbedding
+struct DelayEmbedding{γ, H} <: AbstractEmbedding
     delays::SVector{γ, Int}
+    h::SVector{γ, H}
 end
-
-@inline DelayEmbedding(γ, τ) = DelayEmbedding(Val{γ}(), τ)
-@inline function DelayEmbedding(::Val{γ}, τ::Int) where {γ}
+@inline DelayEmbedding(γ, τ, h = nothing) = DelayEmbedding(Val{γ}(), τ, h)
+@inline function DelayEmbedding(g::Val{γ}, τ::Int, h::H) where {γ, H}
     idxs = [k*τ for k in 1:γ]
-    return DelayEmbedding{γ}(SVector{γ, Int}(idxs...))
+    hs = hweights(g, h)
+    htype = H <: Union{Nothing, Real} ? H : eltype(H)
+    return DelayEmbedding{γ, htype}(SVector{γ, Int}(idxs...), hs)
 end
-@inline function DelayEmbedding(::Val{γ}, τ::AbstractVector) where {γ}
+@inline function DelayEmbedding(g::Val{γ}, τ::AbstractVector, h::H) where {γ, H}
     γ != length(τ) && throw(ArgumentError(
-    "Delay time vector length must equal the number of temporal neighbors."
+        "Delay time vector length must equal the embedding dimension minus 1."
     ))
-    return DelayEmbedding{γ}(SVector{γ, Int}(τ...))
+    hs = hweights(g, h)
+    htype = H <: Union{Nothing, Real} ? H : eltype(H)
+    return DelayEmbedding{γ, htype}(SVector{γ, Int}(τ...), hs)
 end
+hweights(::Val{γ}, h::Nothing) where {γ} = SVector{γ, Nothing}(fill(nothing, γ))
+hweights(::Val{γ}, h::T) where {γ, T<:Real} = SVector{γ, T}([h^b for b in 1:γ]...)
+hweights(::Val{γ}, h::AbstractVector) where {γ} = SVector{γ}(h)
 
-@generated function (r::DelayEmbedding{γ})(s::AbstractArray{T}, i) where {γ, T}
+@generated function (r::DelayEmbedding{γ, Nothing})(s::AbstractArray{T}, i) where {γ, T}
     gens = [:(s[i + r.delays[$k]]) for k=1:γ]
     quote
         @_inline_meta
         @inbounds return SVector{$γ+1,T}(s[i], $(gens...))
     end
 end
 
-export WeightedDelayEmbedding
-"""
-    WeightedDelayEmbedding(γ, τ, w) → `embedding`
-Similar with [`DelayEmbedding`](@ref), but the entries of the
-embedded vector are further weighted with `w^γ`.
-See [`reconstruct`](@ref) for more.
-"""
-struct WeightedDelayEmbedding{γ, T<:Real} <: AbstractEmbedding
-    delays::SVector{γ, Int}
-    w::T
-end
-
-@inline WeightedDelayEmbedding(γ, τ, w) = WeightedDelayEmbedding(Val{γ}(), τ, w)
-@inline function WeightedDelayEmbedding(::Val{γ}, τ::Int, w::T) where {γ, T}
-    idxs = [k*τ for k in 1:γ]
-    return WeightedDelayEmbedding{γ, T}(SVector{γ, Int}(idxs...), w)
-end
-
-@generated function (r::WeightedDelayEmbedding{γ, T})(s::AbstractArray{X}, i) where {γ, T, X}
-    gens = [:(r.w^($k) * s[i + r.delays[$k]]) for k=1:γ]
+# This is the version with weights
+@generated function (r::DelayEmbedding{γ})(s::AbstractArray{T}, i) where {γ, T, R<:Real}
+    gens = [:(r.h[$k]*(s[i + r.delays[$k]])) for k=1:γ]
     quote
         @_inline_meta
-        @inbounds return SVector{$γ+1,X}(s[i], $(gens...))
+        @inbounds return SVector{$γ+1,T}(s[i], $(gens...))
     end
 end
 
+function WeightedDelayEmbedding(γ, τ, w)
+    @warn "`WeightedDelayEmbedding` is deprecated in favor of `DelayEmbedding`."
+    return DelayEmbedding(γ, τ, w)
+end
+
 """
-    reconstruct(s, γ, τ [, w])
-Reconstruct `s` using the delay coordinates embedding with `γ` temporal neighbors
-and delay `τ` and return the result as a [`Dataset`](@ref). Optionally use weight `w`.
+    embed(s, d, τ [, h])
+Embed `s` using delay coordinates with embedding dimension `d` and delay time `τ`
+and return the result as a [`Dataset`](@ref). Optionally use weight `h`, see below.
 
-Use [`embed`](@ref) for the version that accepts the embedding dimension `D = γ+1`
-instead. Here `τ ≥ 0`, use [`genembed`](@ref) for a generalized version.
+Here `τ > 0`, use [`genembed`](@ref) for a generalized version.
 
 ## Description
-### Single Timeseries
 If `τ` is an integer, then the ``n``-th entry of the embedded space is
 ```math
-(s(n), s(n+\\tau), s(n+2\\tau), \\dots, s(n+γ\\tau))
+(s(n), s(n+\\tau), s(n+2\\tau), \\dots, s(n+(d-1)\\tau))
 ```
-If instead `τ` is a vector of integers, so that `length(τ) == γ`,
+If instead `τ` is a vector of integers, so that `length(τ) == d-1`,
 then the ``n``-th entry is
 ```math
-(s(n), s(n+\\tau[1]), s(n+\\tau[2]), \\dots, s(n+\\tau[γ]))
+(s(n), s(n+\\tau[1]), s(n+\\tau[2]), \\dots, s(n+\\tau[d-1]))
 ```
 
-The reconstructed dataset can have same
+The resulting set can have same
 invariant quantities (like e.g. lyapunov exponents) with the original system
-that the timeseries were recorded from, for proper `γ` and `τ`.
+that the timeseries were recorded from, for proper `d` and `τ`.
 This is known as the Takens embedding theorem [^Takens1981] [^Sauer1991].
-The case of different delay times allows reconstructing systems with many time scales,
+The case of different delay times allows embedding systems with many time scales,
 see[^Judd1998].
 
-*Notice* - The dimension of the returned dataset (i.e. embedding dimension) is `γ+1`!
-
-If `w` (a "weight") is provided as an extra argument, then the entries
-of the embedded vector are further weighted with ``w^\\gamma``[^Farmer1988]
+If provided, `h` can be weights to multiply the entries of the embedded space.
+If `h isa Real` then the embedding is
 ```math
-(s(n), w*s(n+\\tau), w^2*s(n+2\\tau), \\dots,w^\\gamma * s(n+γ\\tau))
+(s(n), h \\cdot s(n+\\tau), w^2 \\cdot s(n+2\\tau), \\dots,w^{d-1} \\cdot s(n+γ\\tau))
 ```
-
-### Multiple Timeseries
-To make a reconstruction out of a multiple timeseries (i.e. trajectory) the number
-of timeseries must be known by type, so `s` must be a `Dataset`.
-
-If the trajectory is for example ``(x, y)`` and `τ` is integer, then the ``n``-th
-entry of the embedded space is
-```math
-(x(n), y(n), x(n+\\tau), y(n+\\tau), \\dots, x(n+γ\\tau), y(n+γ\\tau))
-```
-If `τ` is an `AbstractMatrix{Int}`, so that `size(τ) == (γ, B)`,
-then we have
-```math
-(x(n), y(n), x(n+\\tau[1, 1]), y(n+\\tau[1, 2]), \\dots, x(n+\\tau[γ, 1]), y(n+\\tau[γ, 2]))
-```
-
-*Notice* - The dimension of the returned dataset is `(γ+1)*B`!
+Otherwise `h` can be a vector of length `d-1`, which the decides the weights of each
+entry directly.
 
 ## References
 [^Takens1981] : F. Takens, *Detecting Strange Attractors in Turbulence — Dynamical
@@ -135,22 +113,16 @@ Systems and Turbulence*, Lecture Notes in Mathematics **366**, Springer (1981)
 [^Farmer1988]: Farmer & Sidorowich, [Exploiting Chaos to Predict the Future and Reduce Noise"](http://www.nzdl.org/gsdlmod?e=d-00000-00---off-0cltbibZz-e--00-1----0-10-0---0---0direct-10---4-------0-1l--11-en-50---20-home---00-3-1-00-0--4--0--0-0-11-10-0utfZz-8-00&a=d&cl=CL3.16&d=HASH013b29ffe107dba1e52f1a0c_1245)
 
 """
-function reconstruct(s::AbstractVector{T}, γ, τ) where {T}
-    if γ == 0
-        return Dataset(s)
-    end
-    de::DelayEmbedding{γ} = DelayEmbedding(Val{γ}(), τ)
-    return reconstruct(s, de)
-end
-function reconstruct(s::AbstractVector{T}, γ, τ, w) where {T}
-    if γ == 0
+function embed(s::AbstractVector{T}, d, τ, h::H = nothing) where {T, H}
+    if d == 1
         return Dataset(s)
     end
-    de = WeightedDelayEmbedding(Val{γ}(), τ, w)
-    return reconstruct(s, de)
+    htype = H <: Union{Nothing, Real} ? H : eltype(H)
+    de::DelayEmbedding{d-1, htype} = DelayEmbedding(d-1, τ, h)
+    return embed(s, de)
 end
-@inline function reconstruct(s::AbstractVector{T},
-    de::Union{WeightedDelayEmbedding{γ}, DelayEmbedding{γ}}) where {T, γ}
+
+@inline function embed(s::AbstractVector{T}, de::DelayEmbedding{γ}) where {T, γ}
     r = τrange(s, de)
     data = Vector{SVector{γ+1, T}}(undef, length(r))
     @inbounds for i in r
@@ -176,6 +148,12 @@ times and can reconstruct multiple timeseries efficiently.
 """
 embed(s, D, τ) = reconstruct(s, D-1, τ)
 
+function reconstruct(s, γ, τ)
+    @warn "`reconstruct` is deprecated in favor of `embed`. In general, the old argument "*
+    "`γs` is being phased out throughout the library."
+    return embed(s, γ+1, τ)
+end
+
 
 #####################################################################################
 # Multiple timeseries
@@ -209,6 +187,8 @@ end
 end
 @inline function MTDelayEmbedding(
     ::Val{γ}, τ::AbstractMatrix{<:Integer}, ::Val{B}) where {γ, B}
+    @warn "Multi-timeseries delay embedding via `reconstruct` is deprecated in favor of "*
+    "using `genembed` directly."
     X = γ*B
     γ != size(τ)[1] && throw(ArgumentError(
     "`size(τ)[1]` must equal the number of spatial neighbors."
@@ -226,6 +206,8 @@ end
 @generated function (r::MTDelayEmbedding{γ, B, X})(
     s::Union{AbstractDataset{B, T}, SizedArray{Tuple{A, B}, T, 2, M}},
     i) where {γ, A, B, T, M, X}
+    @warn "Multi-timeseries delay embedding via `reconstruct` is deprecated in favor of "*
+    "using `genembed` directly."
     typeof(s) <: SizedArray && @warn "Using SizedArrays is deprecated. Use Dataset instead."
     gensprev = [:(s[i, $d]) for d=1:B]
     gens = [:(s[i + r.delays[$k, $d], $d]) for k=1:γ for d=1:B]

test/delaytime_test.jl

@@ -4,6 +4,7 @@ using DynamicalSystemsBase, Test
 println("\nTesting delay estimation...")
 
 testval = (val, vmin, vmax) -> @test vmin ≤ val ≤ vmax
+diffeq = (atol = 1e-9, rtol = 1e-9, maxiters = typemax(Int))
 
 @testset "Estimate Delay" begin
     # test exponential decay fit

test/embedding_tests.jl

@@ -9,41 +9,34 @@ println("\nTesting delay embeddings...")
 
     @testset "standard" begin
 
-    	@testset "D = $(D), τ = $(τ)" for D in [1,2], τ in [2,3]
-
-    		R = reconstruct(s, D, τ)
-
+    	@testset "D = $(D), τ = $(τ)" for D in [2,3], τ in [2,3]
+    		R = embed(s, D, τ)
     		@test R[(1+τ):end, 1] == R[1:end-τ, 2]
-    		@test size(R) == (length(s) - τ*D, D+1)
-            @test embed(s, D+1, τ) == R
+    		@test size(R) == (length(s) - τ*(D-1), D)
     	end
     end
 
     @testset "weighted" begin
-
-    	@testset "D = $(D), τ = $(τ)" for D in [1,2], τ in [2,3]
+    	@testset "D = $(D), τ = $(τ)" for D in [2,3], τ in [2,3]
 
             w = 0.1
-    		R1 = reconstruct(s, D, τ)
-    		R2 = reconstruct(s, D, τ, w)
-
-            for γ in 0:D
+    		R1 = embed(s, D, τ)
+    		R2 = embed(s, D, τ, w)
+            for γ in 0:D-1
                 @test  (w^γ) * R1[1, γ+1][1] == R2[1, γ+1]
                 @test  (w^γ) * R1[5, γ+1][1] == R2[5, γ+1]
             end
     	end
     end
 
-
     @testset "multi-time" begin
-
-        D = 2
+        D = 3
         τ1 = [2, 4]
         τ2 = [4, 8]
 
-        R0 = reconstruct(s, D, 2)
-        R1 = reconstruct(s, D, τ1)
-        R2 = reconstruct(s, D, τ2)
+        R0 = embed(s, D, 2)
+        R1 = embed(s, D, τ1)
+        R2 = embed(s, D, τ2)
 
         @test R1 == R0
 
@@ -52,50 +45,9 @@ println("\nTesting delay embeddings...")
         @test R2[:, 1] == R0[1:end-4, 1]
         @test size(R2) == (N-maximum(τ2), 3)
 
-        @test_throws ArgumentError reconstruct(s, 4, τ1)
-
-
+        @test_throws ArgumentError embed(s, 4, τ1)
     end
 
-    @testset "multidim " begin
-        @testset "D = $(D), B = $(B)" for  D in [2,3], B in [2,3]
-
-            τ = 3
-            si = Matrix(data[:,1:B])
-            sizedsi = SizedMatrix{N,B}(si)
-            R = reconstruct(sizedsi, D, τ)
-            tr = Dataset(si)
-            R2 = reconstruct(tr, D, τ)
-
-            @test R == R2
-
-            for dim in 1:B
-                @test R[(1+τ):end, dim] == R[1:end-τ, dim+B]
-                @test R2[(1+τ):end, dim] == R[1:end-τ, dim+B]
-            end
-            @test size(R) == (size(s,1) - τ*D, (D+1)*B)
-        end
-    end
-
-    @testset "multidim multi-time" begin
-
-        taus = [2 3; 4 6; 6 8]
-        data2 = data[:, 1:2]
-        data3 = SizedMatrix{N, 2}(Matrix(data2))
-        R1 = reconstruct(data2, 3, taus)
-        R2 = reconstruct(data3, 3, taus)
-
-        @test R1 == R2
-        @test R1[:, 1] == data2[1:end-8, 1]
-        @test R1[:, 2] == data2[1:end-8, 2]
-        @test R1[:, 3] == data2[3:end-6, 1]
-
-        # test error throws:
-        taus = [0 0 0; 2 3 0; 4 6 0; 6 8 0]
-        @test_throws ArgumentError reconstruct(data2, 5, taus)
-        @test_throws ArgumentError reconstruct(data2, 4, taus)
-
-    end
 end
 
 println("\nTesting generalized embedding...")