Distance Time Warp (DTW) function implementation in mathematica

7

3

The DTW (Distance time warp) is given here http://en.wikipedia.org/wiki/Dynamic_time_warping

There is an implementation of DTW in R as well as MATLAB. However i am unable to find something for DTW in Mathematica.

Can someone help me how to implement dtw in Mathematica ??

The DTW function in R is given here:

function (x, y = NULL, dist.method = "Euclidean", step.pattern = symmetric2, 
window.type = "none", keep.internals = FALSE, distance.only = FALSE, 
open.end = FALSE, open.begin = FALSE, ...) 
lm <- NULL
if (is.null(y)) {
    if (!is.matrix(x)) 
        stop("Single argument requires a global cost matrix")
    lm <- x
}
else if (is.character(dist.method)) {
    x <- as.matrix(x)
    y <- as.matrix(y)
    lm <- proxy::dist(x, y, method = dist.method)
}
else if (is.function(dist.method)) {
    stop("Unimplemented")
}
else {
    stop("dist.method should be a character method supported by proxy::dist()")
}
wfun <- .canonicalizeWindowFunction(window.type)
dir <- step.pattern
norm <- attr(dir, "norm")
if (!is.null(list(...)$partial)) {
        warning("Argument `partial' is obsolete. Use `open.end' instead")
        open.end <- TRUE
    }
    n <- nrow(lm)
    m <- ncol(lm)
    if (open.begin) {
        if (is.na(norm) || norm != "N") {
            stop("Open-begin requires step patterns with 'N' normalization (e.g. asymmetric, or R-J types (c)). See papers in citation().")
        }
        lm <- rbind(0, lm)
        np <- n + 1
        precm <- matrix(NA, nrow = np, ncol = m)
        precm[1, ] <- 0
    }
    else {
        precm <- NULL
        np <- n
    }
    gcm <- globalCostMatrix(lm, step.matrix = dir, window.function = wfun, 
        seed = precm, ...)
    gcm$N <- n
gcm$M <- m
    gcm$call <- match.call()
gcm$openEnd <- open.end
    gcm$openBegin <- open.begin
gcm$windowFunction <- wfun
    lastcol <- gcm$costMatrix[np, ]
if (is.na(norm)) {
}
else if (norm == "N+M") {
    lastcol <- lastcol/(n + (1:m))
}
else if (norm == "N") {
    lastcol <- lastcol/n
}
else if (norm == "M") {
    lastcol <- lastcol/(1:m)
}
gcm$jmin <- m
    if (open.end) {
        if (is.na(norm)) {
            stop("Open-end alignments require normalizable step patterns")
        }
        gcm$jmin <- which.min(lastcol)
    }
    gcm$distance <- gcm$costMatrix[np, gcm$jmin]
if (is.na(gcm$distance)) {
        stop("No warping path exists that is allowed by costraints")
    }
    if (!is.na(norm)) {
        gcm$normalizedDistance <- lastcol[gcm$jmin]
    }
    else {
        gcm$normalizedDistance <- NA
    }
    if (!distance.only) {
        mapping <- backtrack(gcm)
        gcm <- c(gcm, mapping)
    }
    if (open.begin) {
        gcm$index1 <- gcm$index1[-1] - 1
        gcm$index2 <- gcm$index2[-1]
        lm <- lm[-1, ]
        gcm$costMatrix <- gcm$costMatrix[-1, ]
        gcm$directionMatrix <- gcm$directionMatrix[-1, ]
    }
    if (!keep.internals) {
        gcm$costMatrix <- NULL
        gcm$directionMatrix <- NULL
    }
    else {
        gcm$localCostMatrix <- lm
        if (!is.null(y)) {
            gcm$query <- x
            gcm$reference <- y
    }
}
class(gcm) <- "dtw"
return(gcm)

}

curiousman

Posted 2014-10-25T05:29:17.850

Reputation: 79

There is some here http://forums.wolfram.com/mathgroup/archive/2003/Jul/msg00544.html

– Nasser – 2014-10-25T05:39:44.837

Already saw it. It is not helping me that much. I am dealing with numbers itself and not the string – curiousman – 2014-10-25T05:57:35.720

3I believe WarpingDistance is the function you're looking for which is new in v11. – rfrasier – 2017-04-02T11:53:39.007

Answers

6

Dynamic time warping (DTW) is a built-in function by the name of WarpingDistance as of Mathematica 11.

C. E.

Posted 2014-10-25T05:29:17.850

Reputation: 67 448

6

Here is a direct (recursive) implementation of the DTW. You need three pieces: a distance function, construction of the DTW matrix, and a function to find the best path through the DTW matrix.

(*distance function*)
dist[s_, t_] := Abs[s - t];
(*boundary conditions*)
dtw[1, 1] = dist[x[[1]], y[[1]]];
dtw[1, j_] := dtw[1, j] = dist[x[[1]], y[[j]]] + dtw[1, j - 1];
dtw[i_, 1] := dtw[i, 1] = dist[x[[i]], y[[1]]] + dtw[i - 1, 1];
(*main recursion*)
dtw[i_, j_] := dtw[i, j] = dist[x[[i]], y[[j]]] + 
    Min[dtw[i - 1, j - 1], dtw[i - 1, j], dtw[i, j - 1]];
(*finding best path through dtwMatrix*)    
pathFind[{i_, j_}] := Module[{nbhd},
   nbhd = {{i, Min[j + 1, nY]}, {Min[i + 1, nX], j}, {Min[i + 1, nX], 
      Min[j + 1, nY]}};
   nbhd[[First[Ordering[Map[dtwMat[[#[[1]], #[[2]]]] &, nbhd]]]]]];

To use this, set up inputs x and y (here we demonstrate with sinusoids)

x = Table[Sin[2 Pi 6 t], {t, 0, 3, 0.01}];
y1 = Table[Sin[2 Pi 7 t], {t, 0, 1.5, 0.005}];
y2 = Table[Sin[2 Pi 3 t], {t, 1.5, 3, 0.005}]; 
y = Flatten[{y1, y2}]; nX = Length[x]; nY = Length[y];

Distance matrix distMat & cumulative distance matrix dtwMat:

distMat = Outer[dist, x, y];
dtwMat = dtwPath = Outer[dtw, Range[nX], Range[nY]];
bestPath = NestWhileList[pathFind, {1, 1}, (#[[1]] < nX) && (#[[2]] < nY) &];
Map[(dtwPath[[#[[1]], #[[2]]]] = None) &, bestPath];

Display distance matrix and highlight best path:

GraphicsRow[{ArrayPlot[distMat], ArrayPlot[dtwPath, Background -> Blue]}]

enter image description here

Find the length of the shortest path:

dtw[nX, nY]
66.1387

Of course, it's probably easier to use the built in version, as suggested by C. E.!

bill s

Posted 2014-10-25T05:29:17.850

Reputation: 62 963