function D = dtw(x,y,dst,verbose)
% D = dtw(x,y,dst,verbose)
%
% x first sequence
% y second sequence
% dst distance indicator
%     1:    edit distance
%     2:    Euclidean DTW
%     else  Derivative DTW (default)
% verbose plots sequences if enabled
%
% Returns DTW distance between x and y.

if nargin < 2
    error('Too few input parameters.');
end
if nargin < 3
    dst = 0;
end
if nargin < 4
    verbose = 0;
end

switch dst
    case 1
        mydist = @(x,y,i,j) (x(i) ~= y(j))*1.0;
    case 2
        mydist = @(x,y,i,j) (x(i) -  y(j))^2;
    otherwise
        mydist = @(x,y,i,j) (dx(x,i) - dx(y,j))^2;
end

m = length(x);
n = length(y);

gamma = Inf*ones(m+1,n+1);
gamma(1,1) = 0;
for i = 2:m+1
    for j = 2:n+1
        gamma(i,j) = mydist(x,y,i-1,j-1) + min([gamma(i-1,j-1) gamma(i,j-1) gamma(i-1,j)]);
    end
end
D = gamma(end,end);

if verbose
    e = [i-1, j-1];
    while 1
        [~,i_min] = min([gamma(i-1,j-1) gamma(i,j-1) gamma(i-1,j)]);
        if i_min ~= 2
            i = i - 1;
        end
        if i_min ~= 3
            j = j - 1;
        end
        if  i <= 1 || j <= 1
            break
        end
        e = [e; i-1, j-1];
    end

    hold on
    plot(1:m, x, '-r');
    plot(1:m, y, '-b');
    for i = 1:size(e,1)
        plot([e(i,1) e(i,2)+eps], [x(e(i,1)) y(e(i,2))], '-k');
    end
end


function delta_x = dx(x,i)

m = length(x);
if i < 1 || i > m
    delta_x = 0;
    return
end

if i == 1
    % left boundary
    delta_x = x(i + 1) - x(i);
    return
end

if i == m
    % right boundary
    delta_x = x(i) - x(i - 1);
    return
end

delta_x = 1/2*(x(i) - x(i - 1) + (x(i + 1) - x(i - 1))/2);
