Merge pull request #28 from JuliaDynamics/gener_embed

Generalized delay embedding

CHANGELOG.md

@@ -1,3 +1,6 @@
+# v1.4
+* New generalized embeddings: `genembed`, `GeneralizedEmbedding` that can do any combination possible from a dataset.
+* Using `SizedArrays` in multidimensional embeddings is deprecated and will be removed in later versions.
 # v1.3
 * Allow indexing datasets with boolean vectors.
 # v1.2.0

Project.toml

@@ -1,7 +1,7 @@
 name = "DelayEmbeddings"
 uuid = "5732040d-69e3-5649-938a-b6b4f237613f"
 repo = "https://github.com/JuliaDynamics/DelayEmbeddings.jl.git"
-version = "1.3.1"
+version = "1.4.0"
 
 [deps]
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
@@ -21,9 +21,8 @@ StatsBase = "0.24, 0.32"
 julia = "1"
 
 [extras]
-ChaosTools = "608a59af-f2a3-5ad4-90b4-758bdf3122a7"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+DynamicalSystemsBase = "6e36e845-645a-534a-86f2-f5d4aa5a06b4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["ChaosTools", "Test", "OrdinaryDiffEq"]
+test = ["DynamicalSystemsBase", "Test"]

src/DelayEmbeddings.jl

@@ -5,7 +5,7 @@ Basic package used in the ecosystem of DynamicalSystems.jl.
 module DelayEmbeddings
 
 include("dataset.jl")
-include("reconstruction.jl")
+include("embeddings.jl")
 include("various.jl")
 include("neighborhoods.jl")
 

src/reconstruction.jl → src/embeddings.jl

@@ -1,29 +1,31 @@
 using StaticArrays
 using Base: @_inline_meta
-export reconstruct, DelayEmbedding, AbstractEmbedding, MTDelayEmbedding, embed
+export reconstruct, DelayEmbedding, MTDelayEmbedding, embed, τrange
 export WeightedDelayEmbedding
 
 #####################################################################################
-#                        Delay Embedding Reconstruction                             #
+# Univariate Delay Coordinates
 #####################################################################################
 """
     AbstractEmbedding
-Super-type of embedding methods. Use `subtypes(AbstractEmbedding)` for available
-methods.
+Super-type of embedding methods.
 """
-abstract type AbstractEmbedding <: Function end
+abstract type AbstractEmbedding end
 
 """
-    DelayEmbedding(γ, τ) -> `embedding`
+    DelayEmbedding(γ, τ) → `embedding`
 Return a delay coordinates embedding structure to be used as a functor,
 given a timeseries and some index. Calling
 ```julia
 embedding(s, n)
 ```
-will create the `n`-th reconstructed vector of the embedded space, which has `γ`
+will create the `n`-th delay vector of the embedded space, which has `γ`
 temporal neighbors with delay(s) `τ`. See [`reconstruct`](@ref) for more.
 
-**Be very careful when choosing `n`, because `@inbounds` is used internally.**
+**Be very careful when choosing `n`, because `@inbounds` is used internally.
+It must be that `n ≤ length(s) - maximum(τ)`.**
+
+Convience function [`τrange`](@ref) gives all valid `n` indices.
 """
 struct DelayEmbedding{γ} <: AbstractEmbedding
     delays::SVector{γ, Int}
@@ -52,12 +54,15 @@ end
 # Weighted version
 export WeightedDelayEmbedding
 """
-    WeightedDelayEmbedding(γ, τ, w) -> `embedding`
+    WeightedDelayEmbedding(γ, τ, w) → `embedding`
 Similar with [`DelayEmbedding`](@ref), but the entries of the
 embedded vector are further weighted with `w^γ`.
 See [`reconstruct`](@ref) for more.
 
-**Be very careful when choosing `n`, because `@inbounds` is used internally.**
+**Be very careful when choosing `n`, because `@inbounds` is used internally.
+It must be that `n ≤ length(s) - maximum(τ)`.**
+
+Convience function [`τrange`](@ref) gives all valid `n` indices.
 """
 struct WeightedDelayEmbedding{γ, T<:Real} <: AbstractEmbedding
     delays::SVector{γ, Int}
@@ -81,10 +86,10 @@ end
 """
     reconstruct(s, γ, τ [, w])
 Reconstruct `s` using the delay coordinates embedding with `γ` temporal neighbors
-and delay `τ` and return the result as a [`Dataset`](@ref).
+and delay `τ` and return the result as a [`Dataset`](@ref). Optionally use weight `w`.
 
 Use [`embed`](@ref) for the version that accepts the embedding dimension `D = γ+1`
-instead.
+instead. Here `τ ≥ 0`, use [`genembed`](@ref) for a generalized version.
 
 ## Description
 ### Single Timeseries
@@ -115,10 +120,7 @@ of the embedded vector are further weighted with ``w^\\gamma``, like so
 
 ### 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` can be either:
-
-* `s::AbstractDataset{B}`
-* `s::SizedAray{A, B}`
+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
@@ -157,14 +159,20 @@ function reconstruct(s::AbstractVector{T}, γ, τ, w) where {T}
 end
 @inline function reconstruct(s::AbstractVector{T},
     de::Union{WeightedDelayEmbedding{γ}, DelayEmbedding{γ}}) where {T, γ}
-    L = length(s) - maximum(de.delays)
-    data = Vector{SVector{γ+1, T}}(undef, L)
-    @inbounds for i in 1:L
+    r = τrange(s, de)
+    data = Vector{SVector{γ+1, T}}(undef, length(r))
+    @inbounds for i in r
         data[i] = de(s, i)
     end
     return Dataset{γ+1, T}(data)
 end
 
+"""
+    τrange(s, de::AbstractEmbedding)
+Return the range `r` of valid indices `n` to create delay vectors out of `s` using `de`.
+"""
+τrange(s, de::AbstractEmbedding) = 1:(length(s) - maximum(de.delays))
+
 """
     embed(s, D, τ)
 Perform a delay coordinates embedding on signal `s` with embedding dimension `D`
@@ -178,21 +186,24 @@ embed(s, D, τ) = reconstruct(s, D-1, τ)
 
 
 #####################################################################################
-#                              MultiDimensional R                                   #
+# Multiple timeseries
 #####################################################################################
 """
     MTDelayEmbedding(γ, τ, B) -> `embedding`
 Return a delay coordinates embedding structure to be used as a functor,
-given multiple timeseries (`B` in total), either as a [`Dataset`](@ref) or a
-`SizedArray`), and some index.
+that embeds multiple timeseries (`B` in total) given in the form of a [`Dataset`](@ref).
+
 Calling
 ```julia
 embedding(s, n)
 ```
-will create the `n`-th reconstructed vector of the embedded space, which has `γ`
-temporal neighbors with delay(s) `τ`. See [`reconstruct`](@ref) for more.
+where `s` is a `Dataset` will create the `n`-th delay vector of the embedded space,
+which has `γ` temporal neighbors with delay(s) `τ`. See [`reconstruct`](@ref) for more.
+
+**Be very careful when choosing `n`, because `@inbounds` is used internally.
+It must be that `n ≤ length(s) - maximum(τ)`.**
 
-**Be very careful when choosing `n`, because `@inbounds` is used internally.**
+Convience function [`τrange`](@ref) gives all valid `n` indices.
 """
 struct MTDelayEmbedding{γ, B, X} <: AbstractEmbedding
     delays::SMatrix{γ, B, Int, X} # X = γ*B = total dimension number
@@ -223,6 +234,7 @@ 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}
+    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]
     quote
@@ -231,6 +243,8 @@ end
     end
 end
 
+τrange(s, de::MTDelayEmbedding) = 1:(size(s)[1] - maximum(de.delays))
+
 @inline function reconstruct(
     s::Union{AbstractDataset{B, T}, SizedArray{Tuple{A, B}, T, 2, M}},
     γ, τ) where {A, B, T, M}
@@ -245,13 +259,107 @@ end
     s::Union{AbstractDataset{B, T}, SizedArray{Tuple{A, B}, T, 2, M}},
     de::MTDelayEmbedding{γ, B, F}) where {A, B, T, M, γ, F}
 
-    L = size(s)[1] - maximum(de.delays)
+    r = τrange(s, de)
     X = (γ+1)*B
-    data = Vector{SVector{X, T}}(undef, L)
-    @inbounds for i in 1:L
+    data = Vector{SVector{X, T}}(undef, length(r))
+    @inbounds for i in r
         data[i] = de(s, i)
     end
     return Dataset{X, T}(data)
 end
 
 reconstruct(s::AbstractMatrix, args...) = reconstruct(Dataset(s), args...)
+
+#####################################################################################
+# Generalized embedding (arbitrary combination of timeseries and delays)
+#####################################################################################
+export GeneralizedEmbedding, genembed
+
+"""
+    GeneralizedEmbedding(τs, js) -> `embedding`
+Return a delay coordinates embedding structure to be used as a functor.
+Given a timeseries *or* trajectory (i.e. `Dataset`) `s` and calling
+```julia
+embedding(s, n)
+```
+will create the `n`-th delay vector of `s` in the embedded space using
+`generalized` embedding (see [`genembed`](@ref).
+
+`js` is ignored for timeseries input `s` (since all entries of `js` must be `1` in
+this case).
+
+**Be very careful when choosing `n`, because `@inbounds` is used internally.
+It must be that `minimum(τs) + 1 ≤ n ≤ length(s) - maximum(τs)`.
+In addition please ensure that all entries of `js` are valid dimensions of `s`.**
+
+Convience function [`τrange`](@ref) gives all valid `n` indices.
+"""
+struct GeneralizedEmbedding{D} <: AbstractEmbedding
+    τs::NTuple{D, Int}
+    js::NTuple{D, Int}
+end
+
+function Base.show(io::IO, g::GeneralizedEmbedding{D}) where {D}
+    print(io, "$D-dimensional generalized embedding\n")
+    print(io, "  τs: $(g.τs)\n")
+    print(io, "  js: $(g.js)")
+end
+
+# timeseries input
+@generated function (g::GeneralizedEmbedding{D})(s::AbstractArray{T}, i::Int) where {D, T}
+    gens = [:(s[i + g.τs[$k]]) for k=1:D]
+    quote
+        @_inline_meta
+        @inbounds return SVector{$D,T}($(gens...))
+    end
+end
+
+# dataset input
+@generated function (g::GeneralizedEmbedding{D})(s::Dataset{X, T}, i::Int) where {D, X, T}
+    gens = [:(s[i + g.τs[$k], g.js[$k]]) for k=1:D]
+    quote
+        @_inline_meta
+        @inbounds return SVector{$D,T}($(gens...))
+    end
+end
+
+τrange(s, ge::GeneralizedEmbedding) =
+max(1, (-minimum(ge.τs) + 1)):min(length(s), length(s) - maximum(ge.τs))
+
+
+"""
+    genembed(s, τs, js = ones(...)) → dataset
+Create a generalized embedding of `s` which can be a timeseries or arbitrary `Dataset`
+and return the result as a new `dataset`.
+
+The generalized embedding works as follows:
+- `τs::NTuple{D, Int}` denotes what delay times will be used for each of the entries
+  of the delay vector. It is strongly recommended that `τs[1] = 0`.
+  `τs` is allowed to have *negative entries* as well.
+- `js::NTuple{D, Int}` denotes which of the timeseries contained in `s`
+  will be used for the entries of the delay vector. `js` can contain duplicate indices.
+
+For example, imagine input trajectory ``s = [x, y, z]`` where ``x, y, z`` are timeseries
+(the columns of the `Dataset`).
+If `js = (1, 3, 2)` and `τs = (0, 2, -7)` the created delay vector at
+each step ``n`` will be
+```math
+(x(n), z(n+2), y(n-7))
+```
+
+`js` can be skipped, defaulting to index 1 (first timeseries) for all delay entries.
+
+See also [`reconstruct`](@ref). Internally uses [`GeneralizedEmbedding`](@ref).
+"""
+function genembed(s, τs::NTuple{D, Int}, js::NTuple{D, Int}) where {D}
+    ge::GeneralizedEmbedding{D} = GeneralizedEmbedding(τs, js)
+    r = τrange(s, ge)
+    T = eltype(s)
+    data = Vector{SVector{D, T}}(undef, length(r))
+    @inbounds for (i, n) in enumerate(r)
+        data[i] = ge(s, n)
+    end
+    return Dataset{D, T}(data)
+end
+
+genembed(s, τs::NTuple{D, Int}) where {D} = genembed(s, τs, NTuple{D, Int}(ones(D)))

src/estimate_delay.jl

@@ -160,7 +160,7 @@ export mutualinformation
 
 Calculate the mutual information between the time series `s` and its images
 delayed by `τ` points for `τ` ∈ `τs`, using an _improvement_ of the method
-outlined by Fraser & Swinney in [1].
+outlined by Fraser & Swinney in[^Fraser1986].
 
 ## Description
 
@@ -188,9 +188,7 @@ For performance and stability reasons, the automatic partition method implemente
 in this function is only used to divide the axes of the grid, using the marginal
 frequencies of `s`.
 
-## References
-[1]: Fraser A.M. & Swinney H.L. "Independent coordinates for strange attractors
-from mutual information" *Phys. Rev. A 33*(2), 1986, 1134:1140.
+[^Fraser1986]: Fraser A.M. & Swinney H.L. "Independent coordinates for strange attractors from mutual information" *Phys. Rev. A 33*(2), 1986, 1134:1140.
 """
 function mutualinformation(s::AbstractVector{T}, τs::AbstractVector{Int};
     kwargs...) where {T}