fast-DTW算法 python实现

之前已经介绍过了DTW算法，现在根据文章 toward accurate dynamic time warping in linear time and space,以及别人实现的fastdtw代码分析fast-DTW算法。参考博客：http://www.cnblogs.com/kemaswill/archive/2013/04/18/3029078.html 
简单讲讲fast-DTW，该算法最主要有两个部分，第一个是约束 
 
将搜索空间约束在阴影位置，减少了搜索的次数。 
第二个是抽象（abstraction）,将图像的像素合并，1/1–>1/2–>1/4–>1/8…知道可以确定路径。如下图，当像素合并为1/8时，已经可以确定路径，即从左下角到右上角。接着再像素粒度细化，从1/8回到1/4确定该像素下的路径，接着1/2，最后1/1。 
 
接下来我们看看Python实现的fastdtw算法。

def fastdtw(x, y, radius=1, dist=lambda a, b: abs(a - b)):
    min_time_size = radius + 2

    if len(x) < min_time_size or len(y) < min_time_size:
        return dtw(x, y, dist=dist)

    x_shrinked = __reduce_by_half(x)
    y_shrinked = __reduce_by_half(y)
    distance, path = fastdtw(x_shrinked, y_shrinked, radius=radius, dist=dist)
    window = __expand_window(path, len(x), len(y), radius)
    return dtw(x, y, window, dist=dist)

 
 
1
2
3
4
5
6
7
8
9
10
11
12
 
 
1
2
3
4
5
6
7
8
9
10
11
12

从代码可以看出来这个过程实现是递归的，并且调用了三个函数，__expand_window，dtw，__reduce_by_half。我们先去弄懂这三个函数，再回头来把整个fastdtw算法理顺。 
首先我们来看这个dtw算法。仔细理一理思路，就能很清晰的分析这个代码了，window就是限制了的搜索范围，只要在这个范围内按照这个条件搜索就行。 
D[i, j] = min((D[i-1, j][0]+dt, i-1, j), (D[i, j-1][0]+dt, i, j-1),(D[i-1, j-1][0]+dt, i-1, j-1), key=lambda a: a[0])
这其中用到了python的collections的defaultdict模块，只要你传入一个默认的工厂方法，那么请求一个不存在的key时， 便会调用这个工厂方法使用其结果来作为这个key的默认值。寻找路径的时候是从终点寻到起点即可。

def dtw(x, y, window=None, dist=lambda a, b: abs(a - b)):
    len_x, len_y = len(x), len(y)
    if window is None:
        window = [(i, j) for i in xrange(len_x) for j in xrange(len_y)]
    window = ((i + 1, j + 1) for i, j in window)
    D = defaultdict(lambda: (float('inf'),))
    D[0, 0] = (0, 0, 0)
    for i, j in window:
        dt = dist(x[i-1], y[j-1])
        D[i, j] = min((D[i-1, j][0]+dt, i-1, j), (D[i, j-1][0]+dt, i, j-1),(D[i-1, j-1][0]+dt, i-1, j-1), key=lambda a: a[0])
    path = []
    i, j = len_x, len_y
    while not (i == j == 0):
        path.append((i-1, j-1))
        i, j = D[i, j][1], D[i, j][2]
    path.reverse()
    return (D[len_x, len_y][0], path)
 
 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
 
 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

接着，看看constraint是怎么实现的，即dtw中的window。这一部分没太弄明白，我运行了一下这个程序，当2x2–> 4x4扩大的时候，path是整个大小，而4x4–>8x8扩大的时候path大小为56。我表示很惊讶，要是我写这个代码，写出来的path大小不会这么多。 

paper里面写时间复杂度的时候提到 

假设我们是在8*8的的大小下寻找搜索空间，按照paper的公式来说应该是56个，而按照图示和描述应该是44个，不信自己可以去数一数。当然我也知道当N趋于无穷时，搜索空间的差异就显得很小了。

def __expand_window(path, len_x, len_y, radius):
    path_ = set(path)
    for i, j in path:
        for a, b in ((i + a, j + b)
                     for a in xrange(-radius, radius+1)
                     for b in xrange(-radius, radius+1)):
            path_.add((a, b))

    window_ = set()
    for i, j in path_:
        for a, b in ((i * 2, j * 2), (i * 2, j * 2 + 1),
                     (i * 2 + 1, j * 2), (i * 2 + 1, j * 2 + 1)):
            window_.add((a, b))

    window = []
    start_j = 0
    for i in xrange(0, len_x):
        new_start_j = None
        for j in xrange(start_j, len_y):
            if (i, j) in window_:
                window.append((i, j))
                if new_start_j is None:
                    new_start_j = j
            elif new_start_j is not None:
                break
        start_j = new_start_j

    return window
 
 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
 
 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

最后返回去看整个代码，fastdtw就是一直递推调用自己，同时将两个序列二分压缩，直到达到可以直接计算出path的大小时，返回每一次调用。然后将粒度细化，就是图二的整个过程。