https://github.com/cran/dtw
Tip revision: 904b9c1b7d11a1ee46753dc815fb71c7ca61d0ac authored by Toni Giorgino on 30 November 2007, 00:00:00 UTC
version 0.1-2
version 0.1-2
Tip revision: 904b9c1
globalCostMatrix.R
########################################
## Compute the cost matrix from a local distance matrix
## $Id:$
`globalCostMatrix` <-
function(lm,
step.pattern="s",
window.type="none",
window.size=10) {
## TODO. Coerce lm to matrix?
wfun<-window.type;
if(!is.function(wfun)) {
wt<-pmatch(window.type,c("sakoechiba","itakura","none"));
if(is.na(wt)) {
wfun<-noWindow;
warning("Unsupported argument for window.type: using no windowing");
} else if(wt==1) { #SAKOE band
wfun<-sakoeChibaWindow;
} else if(wt==2) { #ITAKURA pg.
wfun<-itakuraWindow;
} else if(wt==3) { #no windowing
wfun<-noWindow;
} else {
wfun<-noWindow;
warning("Should not happen on window.type: using no windowing");
}
}
# Step patterns: \delta i, \delta j, cost
# ...
# all deltas MUST be positive (otherwise we violate monotonicity)
# #TODO pmatch
dir<-step.pattern; # custom step pattern
if(dir == "s") { # symmetric DTW
dir<-matrix(c(0,1,1,
1,0,1,
1,1,1), nrow=3, ncol=3,byrow=TRUE);
} else if(dir == "a") { # asymmetric DTW
dir<-matrix(c(1,0,1,
1,1,1,
1,2,1), nrow=3, ncol=3,byrow=TRUE);
}
# i = 1 .. n in query sequence, on first index, ie rows
# j = 1 .. m on template sequence, on second index, ie columns
# Note: template is usually drawn vertically, up-wise
query.size <- n <- dim(lm)[1];
template.size <- m <- dim(lm)[2];
# number of step patterns allowed
nsteps<-dim(dir)[1];
# clear the cost and step matrix
cm <- matrix(NA,nrow=n,ncol=m);
sm <- matrix(NA,nrow=n,ncol=m);
# initializer
cm[1,1] <- lm[1,1];
# now walk through the matrix, column-wise and row-wise,
# and recursively compute the accumulated distance. Unreachable
# elements are handled via NAs (removed)
for (j in 1:m) {
for (i in 1:n) {
d<-lm[i,j];
clist<-c();
slist<-c();
for (s in 1:nsteps) {
## ii,jj is the cell from where potentially we could
## come from. TODO add windowing here.
ii<-i-dir[s,1]; # previous step in inp
jj<-j-dir[s,2]; # previous step in tpl
if(ii>=1 && jj>=1) { # element exists?
cc<- dir[s,3]; # step penalty
clist<-c(clist,cm[ii,jj]+cc*d); #build a cost list
slist<-c(slist,s); #remember whence we came
}
}
## there could no NAs at this point
# store in cost matrix
if(length(which.min(clist))>0) { # odd way to exclude NANs
minc<-which.min(clist); # pick the least cost
cm[i,j]<-clist[minc];
sm[i,j]<-slist[minc];
}
}
}
out<-list();
out$costMatrix<-cm; # to get distance
out$directionMatrix<-sm; # to backtrack
out$stepPatterns<-dir; # to backtrack
return(out);
}