Delay Coordinates Embedding

A timeseries recorded in some manner from a dynamical system can be used to gain information about the dynamics of the entire phase-space of the system. This can be done by reconstructing a new phase-space from the timeseries. One method that can do this is what is known as delay coordinates embedding or delay coordinates reconstruction.

Reconstruct/Embed

Delay embedding reconstructions are done through reconstruct or embed:

DelayEmbeddings.reconstruct — Function

reconstruct(s, γ, τ [, w])

Reconstruct s using the delay coordinates embedding with γ temporal neighbors and delay τ and return the result as a Dataset.

Use embed for the version that accepts the embedding dimension D = γ+1 instead.

Description

Single Timeseries

If τ is an integer, then the $n$-th entry of the embedded space is

\[(s(n), s(n+\tau), s(n+2\tau), \dots, s(n+γ\tau))\]

If instead τ is a vector of integers, so that length(τ) == γ, then the $n$-th entry is

\[(s(n), s(n+\tau[1]), s(n+\tau[2]), \dots, s(n+\tau[γ]))\]

The reconstructed dataset can have same invariant quantities (like e.g. lyapunov exponents) with the original system that the timeseries were recorded from, for proper γ and τ. This is known as the Takens embedding theorem [1, 2]. The case of different delay times allows reconstructing systems with many time scales, see [3].

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$, like so

\[(s(n), w*s(n+\tau), w^2*s(n+2\tau), \dots,w^\gamma * 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 can be either:

s::AbstractDataset{B}
s::SizedAray{A, B}

If the trajectory is for example $(x, y)$ and τ is integer, then the $n$-th entry of the embedded space is

\[(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

\[(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!

References

[1] : F. Takens, Detecting Strange Attractors in Turbulence — Dynamical Systems and Turbulence, Lecture Notes in Mathematics 366, Springer (1981)

[2] : T. Sauer et al., J. Stat. Phys. 65, pp 579 (1991)

[3] : K. Judd & A. Mees, Physica D 120, pp 273 (1998)

DelayEmbeddings.embed — Function

embed(s, D, τ)

Perform a delay coordinates embedding on signal s with embedding dimension D and delay time τ. The result is returned as a Dataset, which is a vector of static vectors.

See reconstruct for an advanced version that supports multiple delay times and can reconstruct multiple timeseries efficiently.

Here are some examples of reconstructing a 3D continuous chaotic system:

using DynamicalSystems, PyPlot

ds = Systems.gissinger(ones(3))
data = trajectory(ds, 1000.0, dt = 0.05)

xyz = columns(data)

figure(figsize = (12,10))
k = 1
for i in 1:3
    for τ in [5, 30, 100]
        R = reconstruct(xyz[i], 1, τ)
        ax = subplot(3,3,k)
        plot(R[:, 1], R[:, 2], color = "C$(k-1)", lw = 0.8)
        title("var = $i, τ = $τ")
        global k+=1
    end
end

tight_layout()
suptitle("2D reconstructed space")
subplots_adjust(top=0.9)

`τ` and `dt`

Keep in mind that whether a value of τ is "reasonable" for continuous systems depends on dt. In the above example the value τ=30 is good, only for the case of using dt = 0.05. For shorter/longer dt one has to adjust properly τ so that their product τ*dt is the same.

You can also reconstruct multidimensional timeseries. For this to be possible, the number of timeseries must be known by Type:

using StaticArrays: Size
a = rand(1000, 3) # my trajectory

A = Size(1000, 3)(a) # create array with the size as Type information
R = reconstruct(A, 2, 2) #aaaall good

9-dimensional Dataset{Float64} with 996 points
 0.256484   0.790223   0.174453   …  0.343874    0.628909   0.863272 
 0.120142   0.373044   0.16791       0.121285    0.985374   0.661154 
 0.799597   0.818934   0.993801      0.0280622   0.0179787  0.924353 
 0.306764   0.262371   0.332753      0.134357    0.189774   0.0180224
 0.343874   0.628909   0.863272      0.114968    0.51378    0.0216869
 0.121285   0.985374   0.661154   …  0.92977     0.661606   0.335877 
 0.0280622  0.0179787  0.924353      0.394555    0.956662   0.986223 
 0.134357   0.189774   0.0180224     0.657785    0.577648   0.142756 
 0.114968   0.51378    0.0216869     0.974613    0.0650369  0.992449 
 0.92977    0.661606   0.335877      0.457956    0.379947   0.0588448
 ⋮                                ⋱                                  
 0.905557   0.226827   0.212771      0.252239    0.0484241  0.105439 
 0.44893    0.662725   0.138321      0.311222    0.383862   0.74507  
 0.996653   0.316605   0.218122      0.555507    0.209989   0.455309 
 0.548062   0.386795   0.032705      0.0450945   0.120385   0.734084 
 0.252239   0.0484241  0.105439   …  0.206956    0.497334   0.81929  
 0.311222   0.383862   0.74507       0.593817    0.110002   0.333277 
 0.555507   0.209989   0.455309      0.783529    0.347379   0.706067 
 0.0450945  0.120385   0.734084      0.722248    0.592311   0.44646  
 0.206956   0.497334   0.81929       0.00678783  0.412704   0.850435

ds = Systems.towel(); tr = trajectory(ds, 10000)
R = reconstruct(tr, 2, 2) # Dataset size is also known by Type!

9-dimensional Dataset{Float64} with 9997 points
 0.085     -0.121       0.075     …  0.837347   0.0372633   0.555269
 0.285813  -0.0675286   0.238038     0.51969    0.0616256   0.940906
 0.76827   -0.038933    0.672094     0.966676  -0.00171595  0.2225  
 0.681871   0.0508933   0.825263     0.112748   0.0674955   0.653573
 0.837347   0.0372633   0.555269     0.386547  -0.0886542   0.869349
 0.51969    0.0616256   0.940906  …  0.910741  -0.0316828   0.411607
 0.966676  -0.00171595  0.2225       0.306095   0.0689305   0.909129
 0.112748   0.0674955   0.653573     0.824263  -0.056185    0.326064
 0.386547  -0.0886542   0.869349     0.545332   0.0508239   0.819404
 0.910741  -0.0316828   0.411607     0.954994   0.00453815  0.569534
 ⋮                                ⋱                                 
 0.914702  -0.0315439   0.294266     0.90246    0.0242141   0.539502
 0.289932   0.0641239   0.778698     0.335976   0.0735803   0.943945
 0.793854  -0.0552801   0.664223     0.86657   -0.0497658   0.214728
 0.62671    0.0557527   0.832001     0.430816   0.0535742   0.62743 
 0.90246    0.0242141   0.539502  …  0.936955  -0.0200112   0.894333
 0.335976   0.0735803   0.943945     0.237481   0.0983265   0.353212
 0.86657   -0.0497658   0.214728     0.681538  -0.0476555   0.883219
 0.430816   0.0535742   0.62743      0.836353   0.0363264   0.380351
 0.936955  -0.0200112   0.894333     0.515471   0.0534613   0.898152

Embedding Functors

The high level functions embed, reconstruct utilize a low-level interface for creating embedded vectors on-the-fly. The high level interface simply loops over the low level interface. The low level interface is composed of the following two structures:

DelayEmbeddings.DelayEmbedding — Type

DelayEmbedding(γ, τ) -> `embedding`

Return a delay coordinates embedding structure to be used as a functor, given a timeseries and some index. Calling

embedding(s, n)

will create the n-th reconstructed vector of the embedded space, which has γ temporal neighbors with delay(s) τ. See reconstruct for more.

Be very careful when choosing n, because @inbounds is used internally.

DelayEmbeddings.MTDelayEmbedding — Type

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 or a SizedArray), and some index. Calling

embedding(s, n)

will create the n-th reconstructed vector of the embedded space, which has γ temporal neighbors with delay(s) τ. See reconstruct for more.

Be very careful when choosing n, because @inbounds is used internally.