R/metric.DTW.R In fda.usc: Functional Data Analysis and Utilities for Statistical Computing

Documented in metric.DTWmetric.TWEDmetric.WDTW

#' @title DTW: Dynamic time warping
#'
#' @description Computes distances time warping for functional data
#'
#' @param fdata1 Functional data 1 or curve 1. If \code{fdata} class, the dimension of \code{fdata1$data} object is (\code{n1} x \code{m}), where \code{n1} is the number of curves and \code{m} are the points observed in each curve. #' @param fdata2 Functional data 2 or curve 2. If \code{fdata} class, the dimension of \code{fdata2$data} object is (\code{n2} x \code{m}), where \code{n2} is the number of curves and \code{m} are the points observed in each curve.
#' @param p Lp norm, by default it uses \code{p = 2}
#' @param w Vector of weights with length \code{m}, If \code{w = 1} approximates the metric Lp by Simpson's rule. By default it uses \code{w = 1}
#' @param lambda \code{numeric} lambda value (0 by default)
#' @param wmax \code{numeric} maximum value of weight, (1 by default)
#' @param g  \code{numeric} \code{g=0} (constant), \code{0.05} (linear) by default, 0.25 \code{sigmoid}, 3 two weight values
#' @param nu  \code{numeric} constant value, (0 by default)
#'
#' @aliases metric.DTW metric.WDTW metric.TWED
#' @details Three optins:
#' \itemize{
#'   \item DTW: Dynamic time warping
#'   \item WDTW: Weight Dynamic time warping
#'   \item TWED: twed
#' }
# Futher information in \url{https://timeseriesclassification.com/}
#'
#' @return DTW matrix
#'
#' @examples
#' \dontrun{
#' data(tecator)
#' metric.DTW(tecator$absorp.fdata[1:4,]) #' ab=tecator[[1]] #' D1=fda.usc:::DTW(ab$data[1,],ab$data[2,],p=2) #' aa1=fda.usc:::findPath(D1$D)
#' D2=fda.usc:::DTW(ab$data[1,],ab$data[2,],p=2,w=5)
#' aa2=fda.usc:::findPath(D2$D) #' D3=fda.usc:::WDTW(ab$data[1,],ab$data[2,],p=2,g=0.05) #' aa3=fda.usc:::findPath(D3$D)
#' D4=fda.usc:::TWED(ab$data[1,],ab$data[2,],p=2,lambda=0,nu=0)
#' aa4=fda.usc:::findPath(D4$D) #' par(mfrow=c(2,2)) #' plot(c(ab[1:2])) #' segments(ab$argvals[aa1[,1]],ab[1]$data[aa1[,1]],ab$argvals[aa1[,2]],ab[2]$data[aa1[,2]]) #' plot(c(ab[1:2])) #' segments(ab$argvals[aa2[,1]],ab[1]$data[aa2[,1]],ab$argvals[aa2[,2]],ab[2]$data[aa2[,2]],col=2) #' plot(c(ab[1:2])) #' segments(ab$argvals[aa3[,1]],ab[1]$data[aa3[,1]],ab$argvals[aa3[,2]],ab[2]$data[aa3[,2]],col=3) #' plot(c(ab[1:2])) #' segments(ab$argvals[aa4[,1]],ab[1]$data[aa4[,1]],ab$argvals[aa4[,2]],ab[2]$data[aa4[,2]],col=4) #' } #' #' @references #' Jeong, Y. S., Jeong, M. K., & Omitaomu, O. A. (2011). Weighted dynamic time warping for #' time series classification. \emph{Pattern Recognition}, 44(9), 2231-2240 #' #' @seealso See also \code{\link{semimetric.basis}} and \code{\link{semimetric.NPFDA}} #' #' @author #' Manuel Febrero-Bande, Manuel Oviedo de la Fuente \email{manuel.oviedo@@udc.es} #' #'@keywords cluster #' @rdname metric.DTW #' @export metric.DTW=function(fdata1, fdata2=NULL, p=2, w=min(ncol(fdata1),ncol(fdata2))){ C1=match.call() same=FALSE if (is.null(fdata2)) { fdata2 <- fdata1 same = TRUE } if (is.fdata(fdata1)){ nas1=is.na.fdata(fdata1) if (any(nas1)) warning("fdata1 contains ",sum(nas1), " curves with some NA value \n") if (is.null(fdata2)) {fdata2=fdata1;same=TRUE} else if (!is.fdata(fdata2)) {fdata2=fdata(fdata2,fdata1$argvals,fdata1$rangeval,fdata1$names)}
nas2=is.na.fdata(fdata2)
if (any(nas2)) warning("fdata2 contains ",sum(nas1), " curves with some NA value \n")
data1=fdata1[["data"]]
data2=fdata2[["data"]]
} else {
data1=data.matrix(fdata1)
if (is.null(fdata2)) {data2=data1;same=TRUE} else {data2=data.matrix(fdata2)}
}
#if (all(data1==data2)) same=TRUE
D=matrix(0,nrow=nrow(data1),ncol=nrow(data2))
if (same){
for (i in 1:(nrow(data1)-1)){
for (j in (i+1):nrow(data2)){
D[i,j]=DTW(data1[i,],data2[j,],p=p,w=w)$dist }} D=(D+t(D)) } else { for (i in 1:nrow(data1)){ for (j in 1:nrow(data2)){ D[i,j]=DTW(data1[i,],data2[j,],p=p,w=w)$dist
}}
}
attr(D,"call")="metric.DTW"
attr(D,"par.metric")=list(p=p,w=w)
return(D)
}

#' @rdname metric.DTW
#' @export
metric.WDTW=function(fdata1,fdata2=NULL,p=2,
w=min(ncol(fdata1),ncol(fdata2)),wmax=1,g=0.05){
C1=match.call()
same=FALSE
if (is.null(fdata2)) {
fdata2 <- fdata1
same = TRUE
}
if (is.fdata(fdata1)){
nas1=is.na.fdata(fdata1)
if (any(nas1)) warning("fdata1 contains ",sum(nas1), " curves with some NA value \n")
if (is.null(fdata2)) {fdata2=fdata1;same=TRUE} else if (!is.fdata(fdata2)) {fdata2=fdata(fdata2,fdata1$argvals,fdata1$rangeval,fdata1$names)} nas2=is.na.fdata(fdata2) if (any(nas2)) warning("fdata2 contains ",sum(nas1), " curves with some NA value \n") data1=fdata1[["data"]] data2=fdata2[["data"]] } else { data1=data.matrix(fdata1) if (is.null(fdata2)) {data2=data1;same=TRUE} else {data2=data.matrix(fdata2)} } #if (all(data1==data2)) same=TRUE D=matrix(0,nrow=nrow(data1),ncol=nrow(data2)) if (same){ for (i in 1:(nrow(data1)-1)){ for (j in (i+1):nrow(data2)){ D[i,j]=WDTW(data1[i,],data2[j,],p=p,w=w,wmax=wmax,g=g)$dist
}}
D=(D+t(D))
} else {
for (i in 1:nrow(data1)){
for (j in 1:nrow(data2)){
D[i,j]=WDTW(data1[i,],data2[j,],p=p,w=w,wmax=wmax,g=g)$dist }} } attr(D,"call")="metric.WDTW" attr(D,"par.metric")=list(p=p,w=w,wmax=wmax,g=g) return(D) } #' @rdname metric.DTW #' @export metric.TWED=function(fdata1, fdata2=NULL,p=2,lambda=1,nu=0.05){ C1=match.call() same=FALSE if (is.null(fdata2)) { fdata2 <- fdata1 same = TRUE } if (is.fdata(fdata1)){ nas1=is.na.fdata(fdata1) if (any(nas1)) warning("fdata1 contains ",sum(nas1), " curves with some NA value \n") if (is.null(fdata2)) {fdata2=fdata1;same=TRUE} else if (!is.fdata(fdata2)) {fdata2=fdata(fdata2,fdata1$argvals,fdata1$rangeval,fdata1$names)}
nas2=is.na.fdata(fdata2)
if (any(nas2)) warning("fdata2 contains ",sum(nas1), " curves with some NA value \n")
data1=fdata1[["data"]]
data2=fdata2[["data"]]
} else {
data1=data.matrix(fdata1)
if (is.null(fdata2)) {data2=data1;same=TRUE} else {data2=data.matrix(fdata2)}
}
#if (all(data1==data2)) same=TRUE
D=matrix(0,nrow=nrow(data1),ncol=nrow(data2))
if (same){
for (i in 1:(nrow(data1)-1)){
for (j in (i+1):nrow(data2)){
D[i,j]=TWED(data1[i,],data2[j,],p=p,lambda=lambda,nu=nu)$dist }} D=(D+t(D)) } else { for (i in 1:nrow(data1)){ for (j in 1:nrow(data2)){ D[i,j]=TWED(data1[i,],data2[j,],p=p,lambda=lambda,nu=nu)$dist
}}
}
attr(D,"call")="metric.TWED"
attr(D,"par.metric")=list(p=p,lambda=lambda,nu=nu)
return(D)
}

############################
cumgamma=function(D,w=min(nrow(D),ncol(D))){
n=nrow(D)
m=ncol(D)
DTW=matrix(0,nrow=n+1,ncol=m+1)
DTW[1,]=Inf
DTW[,1]=Inf
DTW[1,1]=0
for (i in 2:(n+1)){for (j in max(2,i-w):min(m+1,i+w)){DTW[i,j]=0}}
for (i in 2:(n+1)){
for (j in max(2,i-w):min(m+1,i+w)){
DTW[i,j]=D[i-1,j-1]+min(DTW[i-1,j],DTW[i,j-1],DTW[i-1,j-1])
}}
return(DTW)
}

############################
findPath=function(D,i=nrow(D),j=ncol(D)){
ca=c(i)
cb=c(j)
while(i>1 & j>1) {
if (i>1 & j>1){
aa=min(D[i-1,j-1],D[i-1,j],D[i,j-1])
if (aa==D[i-1,j-1]) {i=i-1;j=j-1} else if (aa==D[i-1,j]){i=i-1} else {j=j-1}
} else if (j>1){j=j-1} else {i=i-1}
ca=c(ca,i)
cb=c(cb,j)
}
aa=cbind(ca,cb)
aa=aa[-nrow(aa),]-1
return(aa)
}

############################
DTW=function(a,b,w=min(length(a),length(b)),p=2){
D=outer(a,b,"-")
D=abs(D)^p
DTW=cumgamma(D,w=w)
return(list(dist=DTW[length(a)+1,length(b)+1],D=DTW))
}

############################
WDTW=function(a,b,w=length(a),p=2,wmax=1,g=0.05){
#g=0 (constante), 0.05 (linear), 0.25 (sigmoid), 3 (dos pesos)
n=length(b)
wei=outer(1:n,1:n,"-")
wei=wmax/(1+exp(-g*(abs(w)-n/2)))
M=outer(a,b,"-")
D=abs(wei*M)^p
DTW=cumgamma(D,w=w)
return(list(dist=DTW[length(a)+1,length(b)+1],D=DTW))
}

############################
TWED=function(a,b,p=2,lambda=1,nu=1){
n=length(a)
m=length(b)
D=matrix(0,nrow=n+1,ncol=m+1)
D[1,1]=0
D[2,1]=abs(a[1])^p
D[1,2]=abs(b[1])^p
for (i in 3:(n+1)){ D[i,1]=D[i-1,1]+abs(a[i-2]-a[i-1])^p}
for (j in 3:(m+1)){ D[1,j]=D[1,j-1]+abs(b[i-2]-b[i-1])^p}
for (i in 2:(n+1)){
for (j in 2:(m+1)){
if (i>2 & j>2) {d1=D[i-1,j-1]+2*nu*abs(i-j)+abs(a[i-1]-b[j-1])^p+abs(a[i-2]-b[j-2])^p} else {d1=D[i-1,j-1]+nu*abs(i-j)+abs(a[i-1]-b[j-1])^p}
if (i>2) {d2=D[i-1,j]+abs(a[i-1]-a[i-2])^p+lambda+nu} else {d2=D[i-1,j]+abs(a[i-1])^p+lambda}
if (j>2) {d3=D[i,j-1]+abs(b[j-1]-b[j-2])^p+lambda+nu} else {d3=D[i,j-1]+abs(b[j-1])^p+lambda}
D[i,j]=min(d1,d2,d3)
}}
return(list(dist=D[length(a)+1,length(b)+1],D=D))
}


Try the fda.usc package in your browser

Any scripts or data that you put into this service are public.

fda.usc documentation built on Oct. 17, 2022, 9:06 a.m.