Skip to content
Snippets Groups Projects
Commit 7972c536 authored by juigmend's avatar juigmend
Browse files

Upload New File

parent 05c2da8b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 387 additions and 0 deletions
dtw_minwarplag.ipynb 0 → 100644
+ 387
− 0
View file @ 7972c536
%% Cell type:markdown id: tags:
<center>
*******************************************************************************************
### DYNAMIC TIME WARPING,
### MINIMUM-WARP OPTIMAL PATH, AND POINTWISE LAG
<br>
##### 30 JULY 2023
##### Juan Ignacio Mendoza Garay
##### doctoral student
##### Department of Music, Art and Culture Studies
##### University of Jyväskylä
*******************************************************************************************
</center>
#### INFORMATION:
* Description:
Demonstrates the most simple (as in "easy to understand", not necessarily faster)
algorithm for Dynamic Time Warping (also called "classical" version), an algorithm
to traceback the optimal path (also "classical"), and novel algorithms (as far as I am aware)
to traceback the optimal path with minimum time-warping and its pointwise lag.
* Instructions:
Edit the values indicated with an arrow like this: <---
Comment/uncomment or change values as suggested by the comments.
Run the program, close your eyes and hope for the best.
*******************************************************************************************
%% Cell type:code id: tags:
``` python
import numpy as np
from tslearn.metrics import dtw_path
import matplotlib.pyplot as plt
from scipy.spatial.distance import cdist
from time import time
```
%% Cell type:markdown id: tags:
*******************************************************************************************
#### TEST SIGNALS:
%% Cell type:code id: tags:
``` python
# EXAMPLE SET 1:
if True: # <---
x_length = 40 # <---
x = np.arange(0,x_length)
y = np.zeros(x.size)
y_1 = y.copy()
y_2 = y.copy()
y_1[10:20] = np.arange(0,10,1) # <---
#y_2 = y_1 # <--- same shape, no lag
y_2[15:25] = np.arange(0,10,1) # <--- same shape, lag forwards (delay)
#y_2[5:15] = np.arange(0,10,1) # <--- same shape, lag backwards (anticipation)
#y_2[12:27] = np.arange(0,15,1) # <--- different shape, lag forwards (delay)
#y_2[3:18] = np.arange(0,15,1) # <--- different shape, lag backwards (anticipation)
# EXAMPLE SET 2:
# Uses x, y_1, and y_2 of example set 1.
if False: # <---
y_1[10:20] = [1,2,3,4,5,5,4,3,2,1] # <---
y_2[15:25] = [1,2,3,4,5,5,4,3,2,1] # <---
# EXAMPLE SET 3:
# Uses x of example set 1.
if False: # <---
y_1_period = 3 # <---
y_2_period = 2 # <---
y_2_phase_shift = 0.5 # <---
y_1 = np.sin( x / y_1_period )
y_2 = np.sin( ( (x / y_2_period ) + (np.pi * y_2_phase_shift ) ))
# EXAMPLE SET 4:
if False: # <---
y_1 = np.array([8,7,6,5,5,6,7,8]) # <---
y_2 = np.array([5,4,4,5,6,7,8,9]) # <---
x = np.arange(0,y_1.size)
y_1_rs = y_1.reshape(-1,1)
y_2_rs = y_2.reshape(-1,1)
plt.plot(y_1)
plt.plot(y_2); # the semicolon prevents printing the output of the last command
```
%% Cell type:markdown id: tags:
*******************************************************************************************
#### METHOD 1:
Using the 'tslearn' library. Simple, fast, easy. Gets the job done with no hassle.
%% Cell type:code id: tags:
``` python
dm_1 = cdist( y_1_rs, y_2_rs , 'cityblock' )
tic = time()
optimal_path_1, dtw_score_1 = dtw_path(y_1, y_2)
print('computation time of DTW 1 = '+str(time()-tic))
op_1_x = [col[0] for col in optimal_path_1]
op_1_y = [col[1] for col in optimal_path_1]
plt.imshow(dm_1.T, origin='lower')
plt.plot(op_1_x, op_1_y, ':w',linewidth=2)
plt.title('DTW 1 = '+str(round(dtw_score_1,3)));
```
%% Cell type:markdown id: tags:
The figure above shows a heat-map of the distance matrix and the classical DTW optimal path.
%% Cell type:markdown id: tags:
*******************************************************************************************
#### METHOD 2:
Using hand-made, home-brewed, simple, raw and no-BS code written with love by your humble servant.
%% Cell type:code id: tags:
``` python
# REFERENCES:
# https://tslearn.readthedocs.io/en/stable/user_guide/dtw.html
# https://en.wikipedia.org/wiki/Dynamic_time_warping
# Müller, M. (2007). Dynamic time warping. Information retrieval for music and motion, 69-84.
tic = time()
dm_2 = cdist( y_1_rs, y_2_rs , 'cityblock' )
# initialisation:
C = np.empty((dm_2.shape[0]+1,dm_2.shape[1]+1))
C[:] = np.inf
C[0,0] = 0
# main loop:
for i in range(1,C.shape[0]):
for j in range(1,C.shape[1]):
C[i,j] = dm_2[i-1,j-1] + min(C[i-1, j], C[i, j-1], C[i-1, j-1]) # goes barebones
# C[i,j] = dm_2[i-1,j-1]**2 + min(C[i-1, j], C[i, j-1], C[i-1, j-1]) # gets cocky
dtw_score_2 = C[i,j] # barebones
#dtw_score_2 = np.sqrt(C[i,j]) # cocky
# traceback optimal path:
tic_1 = time()
i, j = np.array(C.shape)-2
optimal_path_2 = [[i,j]]
while (i > 0) or (j > 0):
tb = np.argmin((C[i, j], C[i, j + 1], C[i + 1, j]))
if tb == 0:
i -= 1
j -= 1
elif tb == 1:
i -= 1
else: # (tb == 2):
j -= 1
optimal_path_2.append([i,j])
print('computation time of DTW 2 = '+str(time()-tic))
print('computation time of classical optimal path = '+str(time()-tic_1))
op_2_x = [col[0] for col in optimal_path_2]
op_2_y = [col[1] for col in optimal_path_2]
plt.imshow(dm_2.T, origin='lower')
plt.plot(op_2_x, op_2_y, ':w',linewidth=2)
plt.title('DTW 2 = '+str(round(dtw_score_2,3)));
```
%% Cell type:markdown id: tags:
Notice that the 'barebones' (default) version of method 2 returns a DTW based on absolute differences, not on the Euclidean distance as in the 'cocky' version used by the tslearn library.
%% Cell type:code id: tags:
``` python
# traceback minimum-warp optimal path:
tic = time()
i, j = np.array(C.shape)-2
optimal_path_3 = [[i,j]]
while (i > 0) or (j > 0):
this_query = (C[i, j], C[i, j + 1], C[i + 1, j])
this_min = np.min(this_query)
these_argmins = np.where(this_query==this_min)
if these_argmins[0][0] == 0:
if (np.any(these_argmins[0][:] == 1)) and (i > j):
i -= 1
elif (np.any(these_argmins[0][:] == 2)) and (j > i):
j -= 1
else:
i -= 1
j -= 1
elif these_argmins[0][0] == 1:
i -= 1
elif these_argmins[0][0] == 2:
j -= 1
optimal_path_3.append([i,j])
print('computation time of minimum-warp optimal path = '+str(time()-tic))
op_3_x = [col[0] for col in optimal_path_3]
op_3_y = [col[1] for col in optimal_path_3]
plt.imshow(dm_2.T, origin='lower')
plt.plot(op_2_x, op_2_y, ':w',linewidth=2)
plt.plot(op_3_x, op_3_y, ':r', linewidth=2)
plt.title('DTW 2 = '+str(round(dtw_score_2,3)));
```
%% Cell type:markdown id: tags:
The algorithm for the minimum-warp optimal path (in red) tries to get closer to the diagonal.
%% Cell type:code id: tags:
``` python
pm = np.zeros( dm_2.shape )
for i_op in range(0,len(optimal_path_3)):
pm[op_3_x[i_op],op_3_y[i_op]] = 1
ax = plt.gca()
ax.imshow(pm.T, origin='lower')
ax.plot( range(0,dm_2.shape[0]) , ':w',linewidth=1);
xy_grid = np.arange(-0.5,dm_2.shape[0]+0.5)
ax.grid(1,linestyle=':')
ax.set_xticks(xy_grid)
ax.set_yticks(xy_grid)
ax.set_xticklabels('')
ax.set_yticklabels('')
ax.set_title('minimum-warp optimal path');
```
%% Cell type:code id: tags:
``` python
# pointwise lag of minimum-warp optimal path:
pm_r = np.fliplr(pm).T # rotate to get anti-diagonals (lag domain)
lags = np.zeros(pm_r.shape[0])
i_lags = 0
i_adiag = -pm_r.shape[0] + 1
while i_adiag <= pm_r.shape[0]: # iterate through anti-adiagonals
this_adiag = pm_r.diagonal(i_adiag)
#print('this_adiag = %s'%this_adiag)
if any(this_adiag):
i_center = int( np.floor(this_adiag.size / 2) )
if not this_adiag[int(i_center)] : # value where the anti-adiagonal intersects the main adiagonal
i_match = int(np.where(this_adiag)[0]) # where the path intersects this anti-diagonal
lags[i_lags] = i_center - i_match
i_adiag += 2
i_lags += 1
else: # see if there's anything in the next diagonal
i_adiag += 1
plt.plot(lags)
plt.title('pointwise lag');
```
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment