DynamicAxisWarping.jl
Dynamic Time Warping (DTW), matrix profile and related algorithms in Julia.
Dynamic Time Warping is a method used to compare, or measure the "distance" between two signals. Contrary to, e.g., the Euclidean distance between two signals, which operates pointwise, DTW compares the signals by warping them in time to produce the closest possible match, with the constraint that the warping is monotonic. Underneath the hood, DTW solves an
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"
Simple usage
Inputs of dimension larger than 1 will be treated as sequences where time is in the last dimension. When using higherdimensional 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 equivalent to dtw_cost(a, b, dist, 5)
# The Distances.jl interface is supported
d = DTW(radius=5)
d(a,b)
Plotting
dtwplot(a, b, [dist=SqEuclidean()]; transportcost = 1)
matchplot(a, b, [dist=SqEuclidean()])
matchplot2(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)
Example, twodimension timeseries, simple example
using DynamicAxisWarping, Plots, Distances
# Create signals q(u) ∈ ℜ², y(u) ∈ ℜ²
fs = 70
u = collect(range(0,stop=1,step=1/fs))
# Create a template query signal
q, tq = [sin.(2pi .* u) cos.(2pi .* u)] .+ 0.01 .* randn.(), u
# Create a similar signal
y = [sin.(3pi .*u) cos.(3pi .* u)] .+ 0.01 .* randn.()
last_peak = findlast(isapprox.(y[:,2], maximum(y[:,2]),atol=0.05))
y, ty = y[1:last_peak, :], u[1:last_peak]
y[end10:end,:] .+= q[end10:end,:]
y[10:13] .+= 0.5
# Plot signals
kws = (;linewidth=3, zlabel="index", xlabel="signal comp. 1", ylabel="signal comp. 2",
xticks=1:1:1, yticks=1:1:1, asepct_ratio=1, legend=nothing)
cs, cq = theme_palette(:auto).colors.colors[1:2]
orig= plot(eachcol(q)...,1:size(q,1); c=cq, label="query", kws...)
plot!(eachcol(y)..., 1:size(y,1) ; c=cs, label="similar signal", kws...)
# Warp 2D time signals and visualize
cost, i1, i2 = dtw(y', q', SqEuclidean(); transportcost = 1)
kws=(;kws..., legend=nothing)
warped=plot(eachcol(q[i2,:])..., 1:length(i2); c=cq, label="query", kws...);
plot!(eachcol(y[i1,:])..., 1:length(i1); c=cs, label="signal", kws..., linewidth=1);
plot(orig, warped)
# Visualizing matched ℜ² signal points in ℜ²
kws=(;kws..., legend=nothing, ds=3, separation=0, xlabel="signal comp. 1", ylabel="signal comp. 2")
mp1 = matchplot2(y', q'; kws...)
# Or in ℜ³ with a signal index axis
mp2 = matchplot2(y', q'; showindex=true, zlabel="warped signal index", kws...)
plot(mp1, mp2)
Find a short pattern in a long time series
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, e.g., ZNormalizer
as the fifth argument.
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.loc1)]])

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.  Available normalizers are:
ZNormalizer, DiagonalZNormalizer, NormNormalizer
Multithreaded search
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}
Optimizations
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, ZNormalizer, prune_endpoints = true, prune_envelope = true)
# 853.336 ms (25519 allocations: 5.00 MiB)
Differentiable SoftDTW
The SoftDTW 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 arraymutation 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
Generalized DTW
The gdtw
function implements the algorithm from A General Optimization Framework for Dynamic Time Warping, which takes two continuoustime signals x
and y
on the interval [0,1]
, and warps them by means of warping functions ϕ, ψ
, so that x ∘ ϕ ≈ y ∘ ψ
, where either
ψ(s) = 2s  ϕ(s)
(both signals warped symmetrically, the default), or ψ(s)=s
(only the x
signal is warped). The method allows regularization by imposing penalties on ϕ(t)  t
(the "cumulative warping") and on ϕ'(t)  1
(the "instantaneous warping").
using LinearAlgebra
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
Clustering and barycenter averaging
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 "SoftDTW: a Differentiable Loss Function for TimeSeries" or "SpatioTemporal Alignments: Optimal transport through space and time" for a method that produces better barycenters at the expense of a much higher computational cost.
Sparse distance matrix
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)
Matrix Profile
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
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 diagonal"
dtwplot(a,b, transportcost=1.1) # Should be almost completely diagonal
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.
See also the keyword argument postprocess
that allows you to pass a function D_modified = f(D)
that is used to filter/postprocess the cost matrix. Lowpass filtering can be used to effectively remove smallscale warping.
Combine with optimal transport
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.
Distances.jl interface
d = DTW(radius=radius, dist=SqEuclidean()) # Or FastDTW / SoftDTW
d(a,b)
Signal alignment
See SignalAlignment.jl which uses this package to align outofsync signals and signals with different sample rates.
Acknowledgements
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