Dynamic Time Warping (DTW), matrix profile and related algorithms in Julia.
This package supports arbitrary metrics and arbitrary "spaces", i.e., as long as you are passing a vector or higher dimensional array of something that your distance can operate on, you're good to go. Time is always considered to be the last dimension.
This package is registered and can be installed with:
using Pkg
pkg"add DynamicAxisWarping"
Inputs of dimension larger than 1 will be treated as sequences where time is in the last dimension. When using higher-dimensional series, make sure the provided distance accepts them.
Any distance implementing the Distances.jl interface works, as well as functions on the form dist(x,y) -> ℝ
.
using DynamicAxisWarping, Distances, Plots
cost, i1, i2 = dtw(a,b, [dist=SqEuclidean()]; transportcost = 1)
cost, i1, i2 = fastdtw(a,b, dist, radius) # note https://arxiv.org/abs/2003.11246
cost = dtw_cost(a, b, dist, radius) # Optimized method that only returns cost. Supports early stopping, see docstring. Can be made completely allocation free.
# dtw supports arbitrary upper and lower bound vectors constraining the warping path.
imin,imax = radiuslimits(5,20,20), plot([imin imax])
dtw(a, b, dist, imin, imax) # Cost eqivalent to dtw_cost(a, b, dist, 5)
dtwplot(a, b, [dist=SqEuclidean()]; transportcost = 1)
matchplot(a, b, [dist=SqEuclidean()])
Example:
using DynamicAxisWarping, Plots
fs = 70
t = range(0,stop=1,step=1/fs)
y0 = sin.(2pi .*t)
y1 = sin.(3pi .*t)
y = [y0;y1[2:end]] .+ 0.01 .* randn.()
q = [y0;y0[2:end]] .+ 0.01 .* randn.()
y[10:15] .+= 0.5
q[13:25] .+= 0.5
f1 = plot([q y])
f2 = dtwplot(q,y,lc=:green, lw=1)
f3 = matchplot(q,y,ds=3,separation=1)
plot(f1,f2,f3, legend=false, layout=3, grid=false)
The function dtwnn
searches for a pattern in a long time series. By default, it does not normalize the data over each window, to do this, pass normalizer = ZNormalizer
(this only works for 1D and 2D data).
using DynamicAxisWarping, Distances
radius = 5
a = sin.(0.1 .* (1:100)) .+ 0.1 .* randn.()
b = sin.(0.1 .* (1:100_000)) .+ 0.1 .* randn.()
res = dtwnn(a, b, SqEuclidean(), radius, saveall=false, bsf_multiplier=1) # takes < 0.1s # DynamicAxisWarping.DTWSearchResult(0.4625287975222824, 73452, (prune_end = 79108, prune_env = 0))
plot([a b[eachindex(a) .+ (res.loc-1)]])
saveall
causes the entire distance profile to be computed. This will take longer time to compute. It is stored inres.dists
.bsf_multiplier = 1
: If > 1, require lower bound to exceedbsf_multiplier*best_so_far
. This allows you to find several nearby points without having to compute the entire distance profile.
Below is an example of how several long series y ∈ Y
can be searched for the occurance of query q
in a multithreaded fashion, using tmap
from ThreadTools.jl. In this example, we first create a unique workspace object for each thread to save on allocations
using ThreadTools
const workspaces = map(1:Threads.nthreads()) do i
DTWWorkspace(q, dist, radius)
end
@time results = tmap(Y) do y
dtwnn(workspaces[Threads.threadid()], y, showprogress = false)
end
mincost, minind = findmin(results) # special method for Vector{DTWSearchResult}
The following optimizations are implemented.
- Endpoint lower bound pruning
- Envelope lower bound pruning
- DTW early termination
- Online normalization (see
ZNormalizer
) - Sorting of query series
- All algorithms operate on arbitrary precision numbers. If you pass them
Float32
instead ofFloat64
, they can become up to twice as fast.
dtwnn
is fairly performant, below is a small benchmark performed on a 2014 laptop
a = sin.(0.1f0 .* (1:100)) .+ 0.1f0 .* randn.(Float32)
b = sin.(0.1f0 .* (1:1000_000)) .+ 0.1f0 .* randn.(Float32)
@btime dtwnn($a, $b, SqEuclidean(), 5, prune_endpoints = true, prune_envelope = true, normalizer=Val(ZNormalizer))
# 853.336 ms (25519 allocations: 5.00 MiB)
The Soft-DTW algorithm is provided through the function
soft_dtw_cost(a, b, [SqEuclidean()]; γ = 1, transportcost = 1)
γ
is the smoothing parameters and a smaller value of γ
makes the distance closer to the standard DTW distance.
To differentiate w.r.t. the first argument, try
using ReverseDiff
da = ReverseDiff.gradient(a->soft_dtw_cost(a,b; γ=1), a)
Zygote.jl will not work due to the array-mutation limitation.
See also function soft_dtw_cost_matrix
.
The following example illustrates how to calculate a barycenter (generalized average) using Soft DTW and Optim.jl, the result is shown below, together with three instances of the input series
The gdtw
function implements the algorithm from A General Optimization Framework for Dynamic Time Warping, which takes two continuous-time signals x
and y
on the interval [0,1]
, and warps the first into the second by means of a warping function ϕ
, so that x ∘ ϕ ≈ y
. The method allows regularization by imposing penalties on ϕ(t) - t
(the "cumulative warping") and on ϕ'(t)
(the "instantaneous warping").
ts = range(0, stop=4π, length=128)
x = LinearInterpolation(sin.(ts) .+ 0.1 .* randn.())
y = LinearInterpolation(sin.(1.1 .* ts))
norm(x.(ts) - y.(ts)) # 1.7184237220575787
cost, ϕ = gdtw(x,y)
norm(x.(ϕ.(ts)) - y.(ts)) # 0.9266090849096682
barycenter = dba(vector_of_arrays)
result = dbaclust(data, nclust, DTW())
Note that dba
is known to not always produce the best barycenters. See, e.g., soft_dtw_cost
above and "Soft-DTW: a Differentiable Loss Function for Time-Series" or "Spatio-Temporal Alignments: Optimal transport through space and time" for a method that produces better barycenters at the expense of a much higher computational cost.
The early termination and some of the stopping heuristics can be used to efficiently calculate a sparse distance matrix where only the k nearest neighbors are fully caluclated and stored. To this end, we have the function
dists, inds = sparse_distmat(y::Vector{Vector}, k, dist, radius)
This package defines specialized methods for MatrixProfile
.matrix_profile
, making use of early stopping to accelerate the computation of the matrix profile. The interface is
profile = matrix_profile(y, m, DTW(radius, [transportcost]))
transportcost
adds an additional penalty multiplier for "transporting", i.e., deviations from the Euclidean matching. The standard DTW distance does not consider this added cost and the default is 1. A value greater than 1 multiplies the cost of moving horizontally or vertically in the coupling matrix, promoting a diagonal move, corresponding to the standard Euclidean matching. The influence of the transport cost can be visualized with
a = sin.(1:100); b = sin.(1:100) .+ randn.();
dtwplot(a,b, transportcost=1) # Default
dtwplot(a,b, transportcost=1.01) # Should be "more diagnoal"
dtwplot(a,b, transportcost=1.1) # Should be almost completely diagnoal
You can try a transportcost < 1
as well, but then it is preferable to make weird alignments and I'm not sure how much sense that would make.
The distance between two datapoints can be any distance supporting the Distances.jl interface.
See the file frequency_warping.jl
(notebook) for an example combining dynamic time warping with optimal transport along the frequency axis for spectrograms. This example makes use of SpectralDistances.jl.
d = DTW(radius=radius, dist=SqEuclidean()) # Or FastDTW / SoftDTW
d(a,b)
This package is a fork of https://github.com/ahwillia/TimeWarp.jl which is no longer maintained.
Special thanks to Joseph Fowler (@joefowler) who contributed a substantial portion of the initial code for TimeWarp.jl