Multidimensional Dynamic Time Warping Implementation in Python - confirm?

Question

I believe that I implemented MDTW in python here but I don't know if I did it correctly. The results seem intuitive. Can someone look at this code and tell me if you see anything wrong?

A lot of the papers I read on the subject kinda go over my head, but this just takes the mse of all column vectors in matrix s1 and s2 at each time point inside a given window w instead of single dimensional dynamic time warping which just takes the mse error of each value between two vectors and finds the smallest one.

I use panda's Panel, which is essentially a 3D version of a DataFrame.

Also, after I have this, I want to do clustering with some multidimensional time series. Thoughts on which clustering algorithm to run? Kmeans? Hierarchical? I'll start by building a dendrogram at least.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

def MDTWDistance(s1, s2, window=10, num_columns=1):
    DTW={}

    w = max(window, abs(len(s1)-len(s2)))

    for i in range(-1,len(s1)):
        for j in range(-1,len(s2)):
            DTW[(i, j)] = float('inf')
    DTW[(-1, -1)] = 0

    for i in range(len(s1)):
        for j in range(max(0, i-w), min(len(s2), i+w)):
            #print "Finding Distance of", s1.loc[i], s2.loc[j]
            dist= mdist(s1.loc[i], s2.loc[j], num_columns)
            #print "Dist", dist
            #print i, j, dist
            DTW[(i, j)] = dist + min(DTW[(i-1, j)],DTW[(i, j-1)], DTW[(i-1, j-1)])

    return np.sqrt(DTW[len(s1)-1, len(s2)-1])

def mdist(a, b, num_col):
    dist = 0
    for col in range(num_col):
        #print "Finding Distance of", a[col], b[col]
        dist = dist + (a[col]-b[col])**2
    return dist

x=np.linspace(0,50,100)
ts1=pd.Series(3.1*np.sin(x/1.5)+3.5)
ts2=pd.Series(2.2*np.sin(x/3.5+2.4)+3.2)
ts3=pd.Series(0.04*x+8.0)
ts4=pd.Series(0.048*x+8.6)
ts5=pd.Series(-0.17*x+4.1)
ts6=pd.Series(-0.14*x+4.5)

ts1.plot()
ts2.plot()
ts3.plot()
ts4.plot()
ts5.plot()
ts6.plot()

plt.ylim(-4,12)
plt.legend(['ts1','ts2','ts3','ts4','ts5','ts6'])
plt.show()

timeSeries = pd.Panel({0:pd.DataFrame(np.transpose([ts1, ts2])),
                 1:pd.DataFrame(np.transpose([ts3, ts4])),
                   2:pd.DataFrame(np.transpose([ts5, ts6]))
                  })

print "0 and 1:",MDTWDistance(timeSeries[0], timeSeries[1],window=10, num_columns=2)
print "0 and 2:",MDTWDistance(timeSeries[0], timeSeries[2],window=10, num_columns=2)
print "1 and 2:",MDTWDistance(timeSeries[1], timeSeries[2],window=10, num_columns=2)

The output is as follow:

0 and 1: 89.354619036
0 and 2: 58.8268328591
1 and 2: 133.434513377

With the graph:

EDIT: I found the distance using an MDTW package in R with the following code:

import rpy2.robjects.numpy2ri
from rpy2.robjects.packages import importr
rpy2.robjects.numpy2ri.activate()

# Set up our R namespaces
R = rpy2.robjects.r
DTW = importr('dtw')

# Generate our data
idx = np.linspace(0, 2*np.pi, 100)
template = np.cos(idx)
query = np.sin(idx) + np.array(R.runif(100))/10

for i,j in ([0,1], [0,2], [1,2]):
    print "From", i, "to", j, R.dtw(np.transpose(np.array(timeSeries[i])), np.transpose(np.array(timeSeries[j])), keep=True).rx('distance')[0][0]

which output:

From 0 to 1 186.623310713
From 0 to 2 119.769089068
From 1 to 2 272.849560995

So these numbers are about double the distances I have... so I'm not far off, but there's something wrong with what I'm doing. Halp.

Answer 1

Could the difference be caused by the step pattern you are using? Dtw in R defaults to the symmetric2 step pattern. Explanation of the step pattern from Toni Giorgino's paper:

symmetric2 is normalizable, symmetric, with no local slope constraints. Since one diagonal step costs as much as the two equivalent steps along the sides, it can be normalized dividing by $N+M$ (query+reference lengths).

I thought this might be the problem after reading Comparing Dynamic Time Warping in R and Python.

If your problem is the same try adding step=symmetric1 to your R.dtw arguements.

Answer 2

The code shown here is a recursive implementation of dynamic programming used for time series analysis for similiarity, there is though a more optimal implementation named Fast Dynamic Time Warping which significantly reduces time and memory complexity to $O(n)$.

You will see a significant improvement from original DTW although the measured distance will be directly impacted depending on the distance measurement used (usually is the ecludidean distance as shown here). But there are others that would provide better results depending on the time series nature.