7
$\begingroup$

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)

}

$\endgroup$
  • $\begingroup$ There is some here forums.wolfram.com/mathgroup/archive/2003/Jul/msg00544.html $\endgroup$ – Nasser Oct 25 '14 at 5:39
  • $\begingroup$ Already saw it. It is not helping me that much. I am dealing with numbers itself and not the string $\endgroup$ – curiousman Oct 25 '14 at 5:57
  • 3
    $\begingroup$ I believe WarpingDistance is the function you're looking for which is new in v11. $\endgroup$ – rfrasier Apr 2 '17 at 11:53
6
$\begingroup$

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

$\endgroup$
6
$\begingroup$

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.!

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.