動(dòng)態(tài)時(shí)間扭曲算法提供了兩個(gè)速度可能變化的時(shí)間序列之間的距離概念。如果我有 N 個(gè)序列要相互比較，我可以通過成對(duì)應(yīng)用該算法來構(gòu)造一個(gè)具有核對(duì)角線的 NXN 對(duì)稱矩陣。然而，對(duì)于長二維序列來說，這非常慢。因此，我嘗試對(duì)代碼進(jìn)行矢量化以加速該矩陣計(jì)算。重要的是，我還想提取定義最佳對(duì)齊的索引。到目前為止我的成對(duì)比較代碼：import mathimport numpy as npseq1 = np.random.randint(100, size=(100, 2)) #Two dim sequencesseq2 = np.random.randint(100, size=(100, 2))def seqdist(seq1, seq2): # dynamic time warping function ns = len(seq1) nt = len(seq2) D = np.zeros((ns+1, nt+1))+math.inf D[0, 0] = 0 cost = np.zeros((ns,nt)) for i in range(ns): for j in range(nt): cost[i,j] = np.linalg.norm(seq1[i,:]-seq2[j,:]) D[i+1, j+1] = cost[i,j]+min([D[i, j+1], D[i+1, j], D[i, j]]) d = D[ns,nt] # distance matchidx = [[ns-1, nt-1]] # backwards optimal alignment computation i = ns j = nt for k in range(ns+nt+2): idx = np.argmin([D[i-1, j], D[i, j-1], D[i-1, j-1]]) if idx == 0 and i > 1 and j > 0: matchidx.append([i-2, j-1]) i -= 1 elif idx == 1 and i > 0 and j > 1: matchidx.append([i-1, j-2]) j -= 1 elif idx == 2 and i > 1 and j > 1: matchidx.append([i-2, j-2]) i -= 1 j -= 1 else: break matchidx.reverse() return d, matchidx[d,matchidx] = seqdist(seq1,seq2) #try it
1 回答
阿晨1998
TA貢獻(xiàn)2037條經(jīng)驗(yàn) 獲得超6個(gè)贊
這是對(duì)代碼的一種重寫，使其更適合numba.jit. 這并不完全是矢量化解決方案，但我發(fā)現(xiàn)該基準(zhǔn)測(cè)試速度提高了 230 倍。
from numba import jit
from scipy import spatial
@jit
def D_from_cost(cost, D):
# operates on D inplace
ns, nt = cost.shape
for i in range(ns):
for j in range(nt):
D[i+1, j+1] = cost[i,j]+min(D[i, j+1], D[i+1, j], D[i, j])
# avoiding the list creation inside mean enables better jit performance
# D[i+1, j+1] = cost[i,j]+min([D[i, j+1], D[i+1, j], D[i, j]])
@jit
def get_d(D, matchidx):
ns = D.shape[0] - 1
nt = D.shape[1] - 1
d = D[ns,nt]
matchidx[0,0] = ns - 1
matchidx[0,1] = nt - 1
i = ns
j = nt
for k in range(1, ns+nt+3):
idx = 0
if not (D[i-1,j] <= D[i,j-1] and D[i-1,j] <= D[i-1,j-1]):
if D[i,j-1] <= D[i-1,j-1]:
idx = 1
else:
idx = 2
if idx == 0 and i > 1 and j > 0:
# matchidx.append([i-2, j-1])
matchidx[k,0] = i - 2
matchidx[k,1] = j - 1
i -= 1
elif idx == 1 and i > 0 and j > 1:
# matchidx.append([i-1, j-2])
matchidx[k,0] = i-1
matchidx[k,1] = j-2
j -= 1
elif idx == 2 and i > 1 and j > 1:
# matchidx.append([i-2, j-2])
matchidx[k,0] = i-2
matchidx[k,1] = j-2
i -= 1
j -= 1
else:
break
return d, matchidx[:k]
def seqdist2(seq1, seq2):
ns = len(seq1)
nt = len(seq2)
cost = spatial.distance_matrix(seq1, seq2)
# initialize and update D
D = np.full((ns+1, nt+1), np.inf)
D[0, 0] = 0
D_from_cost(cost, D)
matchidx = np.zeros((ns+nt+2,2), dtype=np.int)
d, matchidx = get_d(D, matchidx)
return d, matchidx[::-1].tolist()
assert seqdist2(seq1, seq2) == seqdist(seq1, seq2)
%timeit seqdist2(seq1, seq2) # 1000 loops, best of 3: 365 μs per loop
%timeit seqdist(seq1, seq2) # 10 loops, best of 3: 86.1 ms per loop
以下是一些變化：
cost是使用 計(jì)算的spatial.distance_matrix。的定義
idx被一堆丑陋的 if 語句取代，這使得編譯的代碼更快。min([D[i, j+1], D[i+1, j], D[i, j]])替換為min(D[i, j+1], D[i+1, j], D[i, j])，即我們不取列表的最小值，而是取三個(gè)值的最小值。這導(dǎo)致了令人驚訝的加速jit。matchidx被預(yù)先分配為 numpy 數(shù)組，并在輸出之前截?cái)酁檎_的大小。
添加回答
舉報(bào)
0/150
提交
取消