Time Warp Edit Distance

Last updated

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 , 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 . It was first proposed in 2009 by P.-F. Marteau.




Whereas the recursion is initialized as:



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.


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:]


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

Related Research Articles

Reed–Solomon codes are a group of error-correcting codes that were introduced by Irving S. Reed and Gustave Solomon in 1960. They have many applications, the most prominent of which include consumer technologies such as MiniDiscs, CDs, DVDs, Blu-ray discs, QR codes, data transmission technologies such as DSL and WiMAX, broadcast systems such as satellite communications, DVB and ATSC, and storage systems such as RAID 6.

Wavenumber Spatial frequency of a wave

In the physical sciences, the wavenumber is the spatial frequency of a wave, measured in cycles per unit distance or radians per unit distance. It is analogous to temporal frequency, which is defined as the number of wave cycles per unit time or radians per unit time.

In the special theory of relativity, four-force is a four-vector that replaces the classical force.

Rayleigh–Jeans law Approximation of a black bodys spectral radiance

In physics, the Rayleigh–Jeans law is an approximation to the spectral radiance of electromagnetic radiation as a function of wavelength from a black body at a given temperature through classical arguments. For wavelength , it is:

The Ising model, named after the physicists Ernst Ising and Wilhelm Lenz, is a mathematical model of ferromagnetism in statistical mechanics. The model consists of discrete variables that represent magnetic dipole moments of atomic "spins" that can be in one of two states. The spins are arranged in a graph, usually a lattice, allowing each spin to interact with its neighbors. Neighboring spins that agree have a lower energy than those that disagree; the system tends to the lowest energy but heat disturbs this tendency, thus creating the possibility of different structural phases. The model allows the identification of phase transitions as a simplified model of reality. The two-dimensional square-lattice Ising model is one of the simplest statistical models to show a phase transition.

In mathematical optimization theory, duality or the duality principle is the principle that optimization problems may be viewed from either of two perspectives, the primal problem or the dual problem. If the primal is a minimization problem then the dual is a maximization problem. Any feasible solution to the primal (minimization) problem is at least as large as any feasible solution to the dual (maximization) problem. Therefore, the solution to the primal is an upper bound to the solution of the dual, and the solution of the dual is a lower bound to the solution of the primal. This fact is called weak duality.

Covariant formulation of classical electromagnetism

The covariant formulation of classical electromagnetism refers to ways of writing the laws of classical electromagnetism in a form that is manifestly invariant under Lorentz transformations, in the formalism of special relativity using rectilinear inertial coordinate systems. These expressions both make it simple to prove that the laws of classical electromagnetism take the same form in any inertial coordinate system, and also provide a way to translate the fields and forces from one frame to another. However, this is not as general as Maxwell's equations in curved spacetime or non-rectilinear coordinate systems.

Maxwells equations in curved spacetime Electromagnetism in general relativity

In physics, Maxwell's equations in curved spacetime govern the dynamics of the electromagnetic field in curved spacetime or where one uses an arbitrary coordinate system. These equations can be viewed as a generalization of the vacuum Maxwell's equations which are normally formulated in the local coordinates of flat spacetime. But because general relativity dictates that the presence of electromagnetic fields induce curvature in spacetime, Maxwell's equations in flat spacetime should be viewed as a convenient approximation.

Ellipsoidal coordinates are a three-dimensional orthogonal coordinate system that generalizes the two-dimensional elliptic coordinate system. Unlike most three-dimensional orthogonal coordinate systems that feature quadratic coordinate surfaces, the ellipsoidal coordinate system is based on confocal quadrics.

Inverse Gaussian distribution Family of continuous probability distributions

In probability theory, the inverse Gaussian distribution is a two-parameter family of continuous probability distributions with support on (0,∞).

Newman–Penrose formalism Notation in general relativity

The Newman–Penrose (NP) formalism is a set of notation developed by Ezra T. Newman and Roger Penrose for general relativity (GR). Their notation is an effort to treat general relativity in terms of spinor notation, which introduces complex forms of the usual variables used in GR. The NP formalism is itself a special case of the tetrad formalism, where the tensors of the theory are projected onto a complete vector basis at each point in spacetime. Usually this vector basis is chosen to reflect some symmetry of the spacetime, leading to simplified expressions for physical observables. In the case of the NP formalism, the vector basis chosen is a null tetrad: a set of four null vectors—two real, and a complex-conjugate pair. The two real members asymptotically point radially inward and radially outward, and the formalism is well adapted to treatment of the propagation of radiation in curved spacetime. The Weyl scalars, derived from the Weyl tensor, are often used. In particular, it can be shown that one of these scalars— in the appropriate frame—encodes the outgoing gravitational radiation of an asymptotically flat system.

The linear attenuation coefficient, attenuation coefficient, or narrow-beam attenuation coefficient characterizes how easily a volume of material can be penetrated by a beam of light, sound, particles, or other energy or matter. A coefficient value that is large represents a beam becoming 'attenuated' as it passes through a given medium, while a small value represents that the medium had little effect on loss. The SI unit of attenuation coefficient is the reciprocal metre (m−1). Extinction coefficient is another term for this quantity, often used in meteorology and climatology. Most commonly, the quantity measures the exponential decay of intensity, that is, the value of downward e-folding distance of the original intensity as the energy of the intensity passes through a unit thickness of material, so that an attenuation coefficient of 1 m−1 means that after passing through 1 metre, the radiation will be reduced by a factor of e, and for material with a coefficient of 2 m−1, it will be reduced twice by e, or e2. Other measures may use a different factor than e, such as the decadic attenuation coefficient below. The broad-beam attenuation coefficient counts forward-scattered radiation as transmitted rather than attenuated, and is more applicable to radiation shielding.

Conway–Maxwell–Poisson distribution Probability distribution

In probability theory and statistics, the Conway–Maxwell–Poisson distribution is a discrete probability distribution named after Richard W. Conway, William L. Maxwell, and Siméon Denis Poisson that generalizes the Poisson distribution by adding a parameter to model overdispersion and underdispersion. It is a member of the exponential family, has the Poisson distribution and geometric distribution as special cases and the Bernoulli distribution as a limiting case.

In mathematics, the Littlewood–Richardson rule is a combinatorial description of the coefficients that arise when decomposing a product of two Schur functions as a linear combination of other Schur functions. These coefficients are natural numbers, which the Littlewood–Richardson rule describes as counting certain skew tableaux. They occur in many other mathematical contexts, for instance as multiplicity in the decomposition of tensor products of finite-dimensional representations of general linear groups, or in the decomposition of certain induced representations in the representation theory of the symmetric group, or in the area of algebraic combinatorics dealing with Young tableaux and symmetric polynomials.

In physics, quantum beats are simple examples of phenomena that cannot be described by semiclassical theory, but can be described by fully quantized calculation, especially quantum electrodynamics. In semiclassical theory (SCT), there is an interference or beat note term for both V-type and -type atoms. However, in the quantum electrodynamic (QED) calculation, V-type atoms have a beat term but -types do not. This is strong evidence in support of quantum electrodynamics.

Laser linewidth is the spectral linewidth of a laser beam.

Lomax distribution

The Lomax distribution, conditionally also called the Pareto Type II distribution, is a heavy-tail probability distribution used in business, economics, actuarial science, queueing theory and Internet traffic modeling. It is named after K. S. Lomax. It is essentially a Pareto distribution that has been shifted so that its support begins at zero.

In the Newman–Penrose (NP) formalism of general relativity, independent components of the Ricci tensors of a four-dimensional spacetime are encoded into seven Ricci scalars which consist of three real scalars , three complex scalars and the NP curvature scalar . Physically, Ricci-NP scalars are related with the energy–momentum distribution of the spacetime due to Einstein's field equation.

In statistics, the complex Wishart distribution is a complex version of the Wishart distribution. It is the distribution of times the sample Hermitian covariance matrix of zero-mean independent Gaussian random variables. It has support for Hermitian positive definite matrices.

Birth process

In probability theory, a birth process or a pure birth process is a special case of a continuous-time Markov process and a generalisation of a Poisson process. It defines a continuous process which takes values in the natural numbers and can only increase by one or remain unchanged. This is a type of birth–death process with no deaths. The rate at which births occur is given by an exponential random variable whose parameter depends only on the current value of the process


  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.