Time Warp Edit Distance

${\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}$

import numpy as np


def dlp(A, B, p=2):
    cost = np.sum(np.power(np.abs(A - B), p))
    return np.power(cost, 1 / p)


def twed(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 arguments
    if len(A) != len(timeSA):
        print("The length of A is not equal length of timeSA")
        return None, None

    if len(B) != len(timeSB):
        print("The length of B is not equal length of timeSB")
        return None, None

    if nu < 0:
        print("nu is negative")
        return None, None

    # Add padding
    A = 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 programming
    DP = np.zeros((n, m))

    # Initialize DP Matrix and set first row and column to infinity
    DP[0, :] = np.inf
    DP[:, 0] = np.inf
    DP[0, 0] = 0

    # Compute minimal cost
    for i in range(1, n):
        for j in range(1, m):
            # Calculate and save cost of various operations
            C = np.ones((3, 1)) * np.inf
            # Deletion in A
            C[0] = (
                DP[i - 1, j]
                + dlp(A[i - 1], A[i])
                + nu * (timeSA[i] - timeSA[i - 1])
                + _lambda
            )
            # Deletion in B
            C[1] = (
                DP[i, j - 1]
                + dlp(B[j - 1], B[j])
                + nu * (timeSB[j] - timeSB[j - 1])
                + _lambda
            )
            # Keep data points in both time series
            C[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 Matrix
            DP[i, j] = np.min(C)
    distance = DP[n - 1, m - 1]
    return distance, DP

def backtracking(DP):
    # [ best_path ] = BACKTRACKING ( DP )
    # Compute the most cost-efficient path
    # DP := DP matrix of the TWED function

    x = np.shape(DP)
    i = x[0] - 1
    j = x[1] - 1

    # The indices of the paths are save in opposite direction
    # path = np.ones((i + j, 2 )) * np.inf;
    best_path = []

    steps = 0
    while i != 0 or j != 0:
        best_path.append((i - 1, j - 1))

        C = np.ones((3, 1)) * np.inf

        # Keep data points in both time series
        C[0] = DP[i - 1, j - 1]
        # Deletion in A
        C[1] = DP[i - 1, j]
        # Deletion in B
        C[2] = DP[i, j - 1]

        # Find the index for the lowest cost
        idx = np.argmin(C)

        if idx == 0:
            # Keep data points in both time series
            i = i - 1
            j = j - 1
        elif idx == 1:
            # Deletion in A
            i = i - 1
            j = j
        else:
            # Deletion in B
            i = i
            j = j - 1
        steps = steps + 1

    best_path.append((i - 1, j - 1))

    best_path.reverse()
    return best_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 arguments
    if length(A) ~= length(timeSA)
        warning('The length of A is not equal length of timeSA')
        return
    end
 
    if length(B) ~= length(timeSB)
        warning('The length of B is not equal length of timeSB')
        return
    end
 
    if nu < 0
        warning('nu is negative')
        return
    end
    % Add padding
    A = [0 A];
    timeSA = [0 timeSA];
    B = [0 B];
    timeSB = [0 timeSB];
 
    % Dynamical programming
    DP = zeros(length(A), length(B));
 
    % Initialize DP Matrix and set first row and column to infinity
    DP(1, :) = inf;
    DP(:, 1) = inf;
    DP(1, 1) = 0;
 
    n = length(timeSA);
    m = length(timeSB);
    % Compute minimal cost
    for i = 2:n
        for j = 2:m
            cost = Dlp(A(i), B(j));
         
            % Calculate and save cost of various operations
            C = ones(3, 1) * inf;
         
            % Deletion in A
            C(1) = DP(i - 1, j) + Dlp(A(i - 1), A(i)) + nu * (timeSA(i) - timeSA(i - 1)) + lambda;
            % Deletion in B
            C(2) = DP(i, j - 1) + Dlp(B(j - 1), B(j)) + nu * (timeSB(j) - timeSB(j - 1)) + lambda;
            % Keep data points in both time series
            C(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 Matrix
            DP(i, j) = min(C);
        end
    end
 
    distance = DP(n, m);
 
    % Function to calculate euclidean distance
    function [cost] = Dlp(A, B)
        cost = sqrt(sum((A - B) .^ 2, 2));
    end
 
end

function [path] = backtracking(DP)
    % [ path ] = BACKTRACKING ( DP )
    % Compute the most cost-efficient path
    % DP := DP matrix of the TWED function
 
    x = size(DP);
    i = x(1);
    j = x(2);
 
    % The indices of the paths are save in opposite direction
    path = 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 series
        C(1) = DP(i - 1, j - 1);
        % Deletion in A
        C(2) = DP(i - 1, j);
        % Deletion in B
        C(3) = DP(i, j - 1);
     
        % Find the index for the lowest cost
        [~, idx] = min(C);
     
        switch idx
            case 1
                % Keep data points in both time series
                i = i - 1;
                j = j - 1;
            case 2
                % Deletion in A
                i = i - 1;
                j = j;
            case 3
                % Deletion in B
                i = i;
                j = j - 1;
        end
        steps = steps + 1;
    end
    path(steps, :) = [i j];
 
    % Path was calculated in reversed direction.
    path = path(1:steps, :);
    path = path(end: - 1:1, :);
 
end

• 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.