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.