Time Warp Edit Distance

Time Warp Edit Distance (TWED) is a measure of similarity (or dissimilarity) for discrete time series matching with time 'elasticity'. In comparison to other distance measures, (e.g. DTW (dynamic time warping) or LCS (longest common subsequence problem)), TWED is a metric. Its computational time complexity is ${\displaystyle O(n^{2})}$, but can be drastically reduced in some specific situations by using a corridor to reduce the search space. Its memory space complexity can be reduced to ${\displaystyle O(n)}$. It was first proposed in 2009 by P.-F. Marteau.

Definition

${\displaystyle \delta _{\lambda ,\nu }(A_{1}^{p},B_{1}^{q})=Min{\begin{cases}\delta _{\lambda ,\nu }(A_{1}^{p-1},B_{1}^{q})+\Gamma (a_{p}^{'}\to \Lambda )&{\rm {delete\ in\ A}}\\\delta _{\lambda ,\nu }(A_{1}^{p-1},B_{1}^{q-1})+\Gamma (a_{p}^{'}\to b_{q}^{'})&{\rm {match\ or\ substitution}}\\\delta _{\lambda ,\nu }(A_{1}^{p},B_{1}^{q-1})+\Gamma (\Lambda \to b_{q}^{'})&{\rm {delete\ in\ B}}\end{cases}}}$
whereas

${\displaystyle \Gamma (\alpha _{p}^{'}\to \Lambda )=d_{LP}(a_{p}^{'},a_{p-1}^{'})+\nu \cdot (t_{a_{p}}-t_{a_{p-1}})+\lambda }$
${\displaystyle \Gamma (\alpha _{p}^{'}\to b_{q}^{'})=d_{LP}(a_{p}^{'},b_{q}^{'})+d_{LP}(a_{p-1}^{'},b_{q-1}^{'})+\nu \cdot (|t_{a_{p}}-t_{b_{q}}|+|t_{a_{p-1}}-t_{b_{q-1}}|)}$
${\displaystyle \Gamma (\Lambda \to b_{q}^{'})=d_{LP}(b_{p}^{'},b_{p-1}^{'})+\nu \cdot (t_{b_{q}}-t_{b_{q-1}})+\lambda }$

Whereas the recursion ${\displaystyle \delta _{\lambda ,\nu }}$ is initialized as:
${\displaystyle \delta _{\lambda ,\nu }(A_{1}^{0},B_{1}^{0})=0,}$
${\displaystyle \delta _{\lambda ,\nu }(A_{1}^{0},B_{1}^{j})=\infty \ {\rm {{for\ }j\geq 1}}}$
${\displaystyle \delta _{\lambda ,\nu }(A_{1}^{i},B_{1}^{0})=\infty \ {\rm {{for\ }i\geq 1}}}$
with ${\displaystyle a'_{0}=b'_{0}=0}$

Implementations

An implementation of the TWED algorithm in C with a Python wrapper is available at ^[1]

TWED is also implemented into the Time Series Subsequence Search Python package (TSSEARCH for short) available at .

An R implementation of TWED has been integrated into the TraMineR, a R package for mining, describing and visualizing sequences of states or events, and more generally discrete sequence data. ^[2]

Additionally, cuTWED is a CUDA- accelerated implementation of TWED which uses an improved algorithm due to G. Wright (2020). This method is linear in memory and massively parallelized. cuTWED is written in CUDA C/C++, comes with Python bindings, and also includes Python bindings for Marteau's reference C implementation.

Python

importnumpyasnpdefdlp(A,B,p=2):cost=np.sum(np.power(np.abs(A-B),p))returnnp.power(cost,1/p)deftwed(A,timeSA,B,timeSB,nu,_lambda):# [distance, DP] = TWED( A, timeSA, B, timeSB, lambda, nu )# Compute Time Warp Edit Distance (TWED) for given time series A and B## A      := Time series A (e.g. [ 10 2 30 4])# timeSA := Time stamp of time series A (e.g. 1:4)# B      := Time series B# timeSB := Time stamp of time series B# lambda := Penalty for deletion operation# nu     := Elasticity parameter - nu >=0 needed for distance measure# Reference :#    Marteau, P.; F. (2009). "Time Warp Edit Distance with Stiffness Adjustment for Time Series Matching".#    IEEE Transactions on Pattern Analysis and Machine Intelligence. 31 (2): 306–318. arXiv:cs/0703033#    http://people.irisa.fr/Pierre-Francois.Marteau/# Check if input argumentsiflen(A)!=len(timeSA):print("The length of A is not equal length of timeSA")returnNone,Noneiflen(B)!=len(timeSB):print("The length of B is not equal length of timeSB")returnNone,Noneifnu<0:print("nu is negative")returnNone,None# Add paddingA=np.array([0]+list(A))timeSA=np.array([0]+list(timeSA))B=np.array([0]+list(B))timeSB=np.array([0]+list(timeSB))n=len(A)m=len(B)# Dynamical programmingDP=np.zeros((n,m))# Initialize DP Matrix and set first row and column to infinityDP[0,:]=np.infDP[:,0]=np.infDP[0,0]=0# Compute minimal costforiinrange(1,n):forjinrange(1,m):# Calculate and save cost of various operationsC=np.ones((3,1))*np.inf# Deletion in AC[0]=(DP[i-1,j]+dlp(A[i-1],A[i])+nu*(timeSA[i]-timeSA[i-1])+_lambda)# Deletion in BC[1]=(DP[i,j-1]+dlp(B[j-1],B[j])+nu*(timeSB[j]-timeSB[j-1])+_lambda)# Keep data points in both time seriesC[2]=(DP[i-1,j-1]+dlp(A[i],B[j])+dlp(A[i-1],B[j-1])+nu*(abs(timeSA[i]-timeSB[j])+abs(timeSA[i-1]-timeSB[j-1])))# Choose the operation with the minimal cost and update DP MatrixDP[i,j]=np.min(C)distance=DP[n-1,m-1]returndistance,DP

Backtracking, to find the most cost-efficient path:

defbacktracking(DP):# [ best_path ] = BACKTRACKING ( DP )# Compute the most cost-efficient path# DP := DP matrix of the TWED functionx=np.shape(DP)i=x[0]-1j=x[1]-1# The indices of the paths are save in opposite direction# path = np.ones((i + j, 2 )) * np.inf;best_path=[]steps=0whilei!=0orj!=0:best_path.append((i-1,j-1))C=np.ones((3,1))*np.inf# Keep data points in both time seriesC[0]=DP[i-1,j-1]# Deletion in AC[1]=DP[i-1,j]# Deletion in BC[2]=DP[i,j-1]# Find the index for the lowest costidx=np.argmin(C)ifidx==0:# Keep data points in both time seriesi=i-1j=j-1elifidx==1:# Deletion in Ai=i-1j=jelse:# Deletion in Bi=ij=j-1steps=steps+1best_path.append((i-1,j-1))best_path.reverse()returnbest_path[1:]

MATLAB

function[distance, DP] = twed(A, timeSA, B, timeSB, lambda, nu)% [distance, DP] = TWED( A, timeSA, B, timeSB, lambda, nu )% Compute Time Warp Edit Distance (TWED) for given time series A and B%% A      := Time series A (e.g. [ 10 2 30 4])% timeSA := Time stamp of time series A (e.g. 1:4)% B      := Time series B% timeSB := Time stamp of time series B% lambda := Penalty for deletion operation% nu     := Elasticity parameter - nu >=0 needed for distance measure%% Code by: P.-F. Marteau - http://people.irisa.fr/Pierre-Francois.Marteau/% Check if input argumentsiflength(A)~=length(timeSA)warning('The length of A is not equal length of timeSA')returnendiflength(B)~=length(timeSB)warning('The length of B is not equal length of timeSB')returnendifnu<0warning('nu is negative')returnend% Add paddingA=[0A];timeSA=[0timeSA];B=[0B];timeSB=[0timeSB];% Dynamical programmingDP=zeros(length(A),length(B));% Initialize DP Matrix and set first row and column to infinityDP(1,:)=inf;DP(:,1)=inf;DP(1,1)=0;n=length(timeSA);m=length(timeSB);% Compute minimal costfori=2:nforj=2:mcost=Dlp(A(i),B(j));% Calculate and save cost of various operationsC=ones(3,1)*inf;% Deletion in AC(1)=DP(i-1,j)+Dlp(A(i-1),A(i))+nu*(timeSA(i)-timeSA(i-1))+lambda;% Deletion in BC(2)=DP(i,j-1)+Dlp(B(j-1),B(j))+nu*(timeSB(j)-timeSB(j-1))+lambda;% Keep data points in both time seriesC(3)=DP(i-1,j-1)+Dlp(A(i),B(j))+Dlp(A(i-1),B(j-1))+...nu*(abs(timeSA(i)-timeSB(j))+abs(timeSA(i-1)-timeSB(j-1)));% Choose the operation with the minimal cost and update DP MatrixDP(i,j)=min(C);endenddistance=DP(n,m);% Function to calculate euclidean distancefunction[cost]=Dlp(A, B)cost=sqrt(sum((A-B).^2,2));endend

Backtracking, to find the most cost-efficient path:

function[path]=backtracking(DP)% [ path ] = BACKTRACKING ( DP )% Compute the most cost-efficient path% DP := DP matrix of the TWED functionx=size(DP);i=x(1);j=x(2);% The indices of the paths are save in opposite directionpath=ones(i+j,2)*Inf;steps=1;while(i~=1||j~=1)path(steps,:)=[i;j];C=ones(3,1)*inf;% Keep data points in both time seriesC(1)=DP(i-1,j-1);% Deletion in AC(2)=DP(i-1,j);% Deletion in BC(3)=DP(i,j-1);% Find the index for the lowest cost[~,idx]=min(C);switchidxcase1% Keep data points in both time seriesi=i-1;j=j-1;case2% Deletion in Ai=i-1;j=j;case3% Deletion in Bi=i;j=j-1;endsteps=steps+1;endpath(steps,:)=[ij];% Path was calculated in reversed direction.path=path(1:steps,:);path=path(end:-1:1,:);end

References

1. Marcus-Voß and Jeremie Zumer, pytwed. "Github repository". GitHub . Retrieved 2020-09-11.
2. TraMineR. "Website on the servers of the Geneva University, CH" . Retrieved 2016-09-11.

• Marteau, P.; F. (2009). "Time Warp Edit Distance with Stiffness Adjustment for Time Series Matching". IEEE Transactions on Pattern Analysis and Machine Intelligence. 31 (2): 306–318. arXiv: cs/0703033 . doi:10.1109/TPAMI.2008.76. PMID   19110495. S2CID   10049446.