R/diss.R

Defines functions distanciaCorCruTotal distanciaCorCruLagK diss.PDC diss.FRECHET diss.DTWARP diss.EUCL loo1nn.cv diss noindicesdiss pairwise.diss cluster.evaluation Sim diss.CID diss.NCD diss.CDM compression.lengths diss.COR diss.DWT wavelet.feature.extraction testIgualdadMaharajHCLUST pvalues.clust diss.CORT corrtemporder1 diss.AR.LPC.CEPS cepstral calc.cepstral.coef arma2ar seasontoplain diss.SPEC.ISD diss.SPEC.GLK multidiss.SPEC.ISD multidiss.SPEC.GLK multidiss.SPEC.LLR diss.SPEC.LLR distance.W.DLS distance.W.LK multidiss.interp.SPEC integrate.GLK interp.SPEC.GLK integrate.ISD interp.SPEC.LS leastsquares.spec integrate.divergenceW interp.W.LK interp.SPEC.LOGLIKELIHOOD trapez plotsmoothspec simetricDivergenceW divergenceW likelihood.optim Spectral.AB funcionKh diss.INT.PER diss.PER diss.AR.MAH distance.MAH.SIMP distance.MAH.EXT maharajahextended diss.AR.PIC find_ar_model_force diss.PACF diss.ACF nternal.autocorr.dist check.equal.length.ts ts.sanity.check ts.freq.check common.ts.sanity.check range.normalize z.normalize

Documented in cluster.evaluation diss diss.ACF diss.AR.LPC.CEPS diss.AR.MAH diss.AR.PIC diss.CDM diss.CID diss.COR diss.CORT diss.DTWARP diss.DWT diss.EUCL diss.FRECHET diss.INT.PER diss.NCD diss.PACF diss.PDC diss.PER diss.SPEC.GLK diss.SPEC.ISD diss.SPEC.LLR loo1nn.cv pvalues.clust

#######################################################################################
#################################  AUXILIARY FUNCTIONS ################################
#######################################################################################


z.normalize = function(x) {
  (x - mean(x)) / sd(x)
}

range.normalize = function(x) {
  minim <- min(x)
  maxim <- max(x)
  (x -minim) / (maxim - minim)    
}

.common.ts.sanity.check <- function(x) {
  if (missing(x)) {
    stop("At least one series is missing!")
  }
  if (any(is.na(x))) {
    stop("NA in the series")
  }
  if (!is.numeric(x)) {
    stop("Series must be numeric")
  }
  #check length
  if (length(x) < 2) {
    stop("Incorrect length of the series")
  }
  if (!is.null(dim(x))) {
    stop("Incorrect dimension of the series, please input univarate series")
  }
}
.ts.freq.check <- function(x, y) {
  if (is.ts(x) && is.ts(y)) { #check their frequencies
    cbind(x,y)
  }
}
.ts.sanity.check <- function(x,y) {
  .common.ts.sanity.check(x)
  .common.ts.sanity.check(y)  
  .ts.freq.check(x,y)
}

#check if series have equal length, a requisite of some functions
.check.equal.length.ts <- function(x,y) {
  if (length(x) != length(y)) {
    stop("Time series must have the same length")
  }  
  
}

################################################################################
#######################  AUTOCORRELATION AND PARTIAL ###########################
################################################################################


#weighted distance of acf and pacf coefficients
.internal.autocorr.dist <- function(rhox, rhoy, p=NULL, omega=NULL) {
    if ( length(rhox) != length(rhoy) ) {#check compatible coefficient vectors
        stop("The amount of autocorrelation coefficients must be the same, maybe lag.max greater than the length of one of the series")
    }
    if (is.null(omega)) { #if there is no weighting matrix
        if (!is.null(p)) { #check if there is gemoetrical decay parameter
            omega <- diag(p*(1-p)**(1:length(rhox)))
        }
        else { #no weightinh matrix and no geomtrical decay parameter, use identity matrix
            omega <- diag(length(rhox))
        }
    }
    sqrt(t(rhox - rhoy) %*% omega %*% (rhox - rhoy)) #weighted euclidean distance
}


diss.ACF <- function(x, y ,  p=NULL,  omega=NULL, lag.max=50) {
    .ts.sanity.check(x, y)
    rhox <- acf(x, lag.max=lag.max, plot=FALSE)$acf[-1]
    rhoy <- acf(y, lag.max=lag.max, plot=FALSE)$acf[-1]
    .internal.autocorr.dist( rhox, rhoy, p, omega)
}



diss.PACF <- function(x, y, p = NULL, omega=NULL, lag.max=50) {
    .ts.sanity.check(x, y)
    rhox <- as.vector(pacf(x, lag.max=lag.max, plot=FALSE)$acf)
    rhoy <- as.vector(pacf(y, lag.max=lag.max, plot=FALSE)$acf)
    .internal.autocorr.dist( rhox, rhoy, p, omega)
}





#######################################################
##########  distance Piccolo  #########################
#######################################################


#try to find the ar coefficients of a AR series, if no order found, try forcing
find_ar_model_force = function(x, permissive) {
    arx <- ar(x) #first, try to find automatically
    if (arx$order < 1) { #if no order found, try forcing order 1
        if (permissive) {
            arx <- ar(x, aic=FALSE, order.max = 1)
        }
        if (arx$order < 1) stop("Could not find a valid AR order for the series")
    }
    arx
}


#######################################

diss.AR.PIC <- function(x, y, order.x=NULL, order.y=NULL, permissive=TRUE) {
    .ts.sanity.check(x, y)
    #if order NULL use ar AIC, else use arima fitting of the given order
    PIx <- NULL
    if (is.null(order.x)) { #no ARIMA order, use AR
        PIx <- find_ar_model_force(x, permissive)$ar
    } else { 
        if ((order.x[1]) < 1) stop("The arima order must have AR coefficients, they are used for the distance")
        arim <- arima(x, order.x)
        PIx <- arim$coef[1:order.x[1]] #get the AR coeff off the arima model
    }
    
    PIy <- NULL
    if (is.null(order.y)) { #no ARIMA order, use AR
        PIy <- find_ar_model_force(y, permissive)$ar
    } else {
        if ((order.y[1]) < 1) stop("The arima order must have AR coefficients, they are used for the distance")
        arim <- arima(y, order.y)
        PIy <- arim$coef[1:order.y[1]] #get the AR coeff off the arima model
    }
    
    k <- max(c(length(PIx), length(PIy))) #get the maximun order
    if (k < 1) {
       stop("Could not find any AR coefficients")
    }
    
    PRIMAx <- rep(0,k) #fill with zeroes to the greatest AR order between series x and y (k)
    if (length(PIx) > 0) {
        PRIMAx[1:length(PIx)] <- PIx
    }
    
    PRIMAy <- rep(0,k)
    if (length(PIy) > 0) {
        PRIMAy[1:length(PIy)] <- PIy
    }
    
    as.numeric( dist(rbind(PRIMAx,PRIMAy)) ) #compute the euclidean distance between the zero padded AR coefficients
}


##################################################
########### distance Maharaj  ####################
##################################################

#regression model
maharajahextended <- function( x, k ) {
    
    X <- x[-(1:k)]
    
    TT <- length(x)
    
    Wx <- matrix(ncol=k, nrow=(TT-k))
    
    for (i in 1:(TT -k) ) {
        Wx[i,] <- x[(k +i - 1):(i)]
    }
    result <- list()
    result$X <- X
    result$Wx <- Wx
    result
}

#extended distance see reference article in the documentation
distance.MAH.EXT <- function( x, y, k) {
    MX <- maharajahextended(x, k)
    MY <- maharajahextended(y, k)
    w <- dim(MX$Wx)[2]
    h <- dim( MX$Wx)[1]
    bigW <- matrix(0, nrow=2*h, ncol=2*w)
    for ( j in 1:w) {
        for (i in 1:h) {
            bigW[i,j] <- MX$Wx[i,j]  
        }
    }
    for ( j in 1:w) {
        for (i in 1:h) {
            bigW[h+i,w+j] <- MY$Wx[i,j]  
        }
    }
    Epsil <- cov(cbind(x,y))
    bigW
    Epsil
    Iden <- diag(1, length(x) - k)
    Iden
    V <- kronecker(Epsil, Iden  )
    IV <- solve(V)
    tryCatch( { #these operations can go wrong, it they fail, try a smaller order
        PI <- solve(t(bigW) %*% IV %*% bigW) %*% t(bigW) %*% IV %*% c(MX$X, MY$X)
        R <- cbind( diag(1, k), diag(-1,k) )
        result <- list()
        result$statistic <- t(R %*% PI) %*% solve( R %*% (t(bigW) %*% IV %*% bigW ) %*% t(R)) %*% ( R%*%PI)
        result$p_value <- pchisq(result$statistic, k, lower.tail=F)
        result
    }, error = function(e) { if (k>1) {
        distance.MAH.EXT(x,y,k-1)
    }
      else { stop("Could not find valid AR order")}
    })
}



distance.MAH.SIMP = function( x, y, PIx, PIy, var.pred.x, var.pred.y, k, permissive=TRUE ) {
    PRIMAx <- rep(0,k)   #fill with zeroes
    if (length(PIx) > 0) {
        PRIMAx[1:length(PIx)] <- PIx
    }
    PRIMAy <- rep(0,k)
    if (length(PIy) > 0) {
        PRIMAy[1:length(PIy)] <- PIy
    }
    
    covx <- acf(x, lag.max=k-1, type="covariance", plot=FALSE)$acf
    covy <- acf(y, lag.max=k-1,type="covariance", plot=FALSE)$acf
    Rx <- matrix(nrow=length(covx),ncol=length(covx))
    Ry <- Rx
    for (i in 1:length(covx) ) {
        indices <- (c(i:1, 2:(length(covx)-i+1)))
        Rx[i,] <- covx[ indices[1:length(covx)]]
        Ry[i,] <- covy[ indices[1:length(covx)]]
    }
    V <- (solve(solve(Rx)*(var.pred.x) + solve(Ry)*(var.pred.y)))
    dif <- (PRIMAx - PRIMAy)
    
    D <- sqrt(length(x)) * ( dif %*% V %*% dif)
    
    list(statistic=D, p_value=pchisq(D, k, lower.tail=F))
}


diss.AR.MAH = function( x, y, dependence = FALSE, permissive = TRUE) {
    .ts.sanity.check(x, y)
    arx <- find_ar_model_force(x, permissive)
    ary <- find_ar_model_force(y, permissive)
    k <- max(c(length(arx$ar), length(ary$ar)))
    if (k < 1) {
        if (k < 1) stop("Could not find a valid AR order for the series")
    }

    if (dependence) {
        .check.equal.length.ts(x,y)
        distance.MAH.EXT(x, y, k)
    }
    else {
        distance.MAH.SIMP(x, y, arx$ar, ary$ar, arx$var.pred, ary$var.pred, k, permissive)
    }
}

######################################################
########  PERIODOGRAM BASED DISTANCES  ###############
######################################################

diss.PER <- function(x,y, logarithm=FALSE, normalize=FALSE) {
    .ts.sanity.check(x, y)
    .check.equal.length.ts(x,y)
    Ix <- spec.pgram(x, plot=F)$spec
    Iy <- spec.pgram(y, plot=F)$spec
    if (normalize) {
        Ix <- Ix/var(x)
        Iy <- Iy/var(y)
    }
    if (logarithm) {
        Ix <- log(Ix)
        Iy <- log(Iy)
    }
    dist(rbind(Ix,Iy))/(length(Ix))
}


diss.INT.PER <- function(x,y, normalize=TRUE) {
    .ts.sanity.check(x, y)
    .check.equal.length.ts(x,y)
    Ix <- spec.pgram(x, plot=F)
    Iy <- spec.pgram(y, plot=F)
    Cx <- 1
    Cy <- 1
    if (normalize) {
        Cx <- sum(Ix$spec)
        Cy <- sum(Iy$spec)
    }
    sum ( abs(cumsum(Ix$spec)/Cx - cumsum(Iy$spec)/Cy) )
}


################################################
###  SPECTRAL DENSITY APPROXMATION DISTANCES ###
################################################


### maximum likelihood functions ###

#kernel function
funcionKh <- function( value, h) {
    value <- value/h
    dnorm( -(value**2) ) / h
}

#function to be optimized for maximum likelihood, see referenced papers
Spectral.AB <- function (  ABvec, lambda, Yks, lambdas, h) {
    acum <- 0
    a <- ABvec[1]
    b <- ABvec[2]
    
    -sum ( ( -exp(Yks - a -b*(lambdas - lambda) ) + Yks - a -b*(lambdas - lambda) ) * funcionKh( lambdas - lambda, h) )
    
}

likelihood.optim <- function(  lambda, Yks, lambdas, h) {
    
    startA = lambdas[1]
    startB = 0
    
    optim(c(startA,startB), Spectral.AB, lambda=lambda, Yks = Yks, lambdas=lambdas, method="L-BFGS-B",lower = c(min(Yks) , -101), upper= c(max(Yks),101), h =h)$par[1]
}
### end of maximul likelihood


#diveregence function
divergenceW <- function( x, alpha) {
    if ((alpha > 0) & (alpha< 1)) {
        log(alpha*x + 1 - alpha) - alpha*log(x)
    }
    else {
        stop("condition 0 < alpha < 1 not fulfilled")
    }
}

simetricDivergenceW <- function(x,alpha) {
    divergenceW(x,alpha) + divergenceW(1/x,alpha)  
}

#plot the soothed spectral density
plotsmoothspec <- function ( lambdasX, YksX, lambdasY, YksY, myf, hX, hY, n, method="Maximum Likelihood") {
    if ( n < 1) {
        n <- 500
    }
    baseX <- seq(min(lambdasX) + 0.001, max(lambdasX) - 0.001, length.out=min(500,n) )
    specX <- NULL
    if (pmatch(method , c("Maximum Likelihood", "Least Squares")) == 2) {
        specX <- (myf(YksX, lambdasX,  hX, baseX))
    } else {
        specX <- exp(myf(baseX, YksX, lambdasX, hX))
    }
    baseY <- seq(min(lambdasY)+ 0.001, max(lambdasY)- 0.001, length.out=min(500, n) )
    specY <- NULL
    if (pmatch(method , c("Maximum Likelihood", "Least Squares")) == 2) {
        specY <- (myf(YksY,lambdasY, hY, baseY))
    } else {
        specY <- exp(myf(baseY, YksY, lambdasY, hY))
    }
    plot.default(baseX, specX,type="l", col="red", ylim=c( min(specX,specY), max(specX, specY) ) ,
                 main=paste("Approx. spectral density by ", method), xlab="frequency",ylab="spectrum")
    lines(baseY, specY,type="l", col="blue")
    legend("topright", pch=16, col=c("red", "blue"), legend= c("x", "y") )
    
    
}

.vectorized.lk.optim <- Vectorize(likelihood.optim,"lambda")  #needed for likelihood.optim to accept a vector, required for integrate


#trapezoind integration, taken from another package
trapez <- function(x,y) {
    idx = 2:length(x)
    return (as.double( (x[idx] - x[idx-1]) %*% (y[idx] + y[idx-1])) / 2)
}

interp.SPEC.LOGLIKELIHOOD <- function(x, n) {
    pgram <- spec.pgram(x, plot=FALSE)
    Yks <- log(pgram$spec)
    lambdas <- pgram$freq
    interplambdas <- seq(min(lambdas), max(lambdas), length.out=n)
    hX <- 0.93*dpill(lambdas, Yks)
    approx( interplambdas, .vectorized.lk.optim( interplambdas, Yks, lambdas, hX ) )
}

interp.W.LK <- function(x, n) {
    interps <- interp.SPEC.LOGLIKELIHOOD(x, n)
    interps$y <- exp(interps$y)
    interps
}

integrate.divergenceW <- function(base, x, y, alpha ) {
    val <- simetricDivergenceW( x / y, alpha)
    trapez(base, val)
}

leastsquares.spec <- function( Yk, lambdas, h, lambdaeval) {
    d <- data.frame(Yk)
    d$lambdas <- lambdas
    lp <- locpol(Yk~lambdas, d, bw=h,kernel=gaussK, xeval=lambdaeval )
    ord <- order(lambdaeval) #trick to get the original order of the lambas, locpol sorts the input vector
    ord2 <- order(ord)   #second part of the trick
    lp$lpFit$Yk[ord2]
}

interp.SPEC.LS <- function(x, n)  {
    pgram <- spec.pgram(x, plot=FALSE)
    Yk <- pgram$spec
    lambdas <- pgram$freq
    h <- dpill(lambdas, Yk)
    interplambdas <- seq(min(lambdas), max(lambdas), length.out=n)
    ys <- leastsquares.spec(Yk, lambdas, h, interplambdas)
    ys[ys<0.0001] <- 0.0001 #no zeroes allowed
    list( x = interplambdas, y = ys)
}

integrate.ISD <- function(base, x, y) {
    trapez(base, (x - y)^2)
}

interp.SPEC.GLK <- function(x, n) {
    pgram <- spec.pgram(x, plot=FALSE)
    Yk <- log(pgram$spec)
    lambdas <- pgram$freq
    h <- 0.93*dpill(lambdas, Yk)
    ys <- .vectorized.lk.optim( lambdas, Yk, lambdas, h )
    list( x = lambdas, y = list(mu = ys, Z = Yk) )
}

integrate.GLK <- function( base, x, y) {
    Z <- x$Z - y$Z
    mu <- x$mu - y$mu
    sum(Z - mu - 2*log(1 + exp(Z - mu))) - sum( Z - 2*log(1 + exp(Z)))
}
#generic linear interpolation approximation of the spectral dissimilarities
#interpfun is a function to calculate the spectrum aproximation
#n the size of the grid for interpolation
#intergrationfun is calculates the sum/integration of the diferences
multidiss.interp.SPEC <- function( series, n, interpfun, integrationfun, ...) {
    l <- length(series)
    dists <- matrix(0, l, l)
    
    if ( sum(diff(unlist(lapply(series, function(x) length(x))))) != 0 ) {
      stop("SPECtral dissimilarity requires equal length time series")
    }
    
    #get the interpolated values
    interps <- lapply(series, interpfun, n)
    minbase <- unlist(lapply(interps, function(x) {min(x$x)}))
    maxbase <- unlist(lapply(interps, function(x) {max(x$x)}))
    base <- seq( max(minbase), min(maxbase), length.out=n)
    ##calc the function with the interpolated values
    for (i in 1:(l-1)) {
        for (j in (i+1):l) {
            d <- integrationfun( base, interps[[i]]$y , interps[[j]]$y, ...)
            dists[i,j] <- d
            dists[j,i] <- d
        }
    }
    as.dist(dists)
}

distance.W.LK <- function(x,y, alpha, plot=FALSE, n=length(x)) {
    
    
    
    pgx <- spec.pgram(x,plot=FALSE)
    YksX <- log(pgx$spec)
    lambdasX <- pgx$freq
    pgy <- spec.pgram(y,plot=FALSE)
    YksY <- log(pgy$spec)
    lambdasY <- (pgy$freq)
    
    hX <- 0.93*dpill(lambdasX, YksX)
    hY <- 0.93*dpill(lambdasY, YksY)
    
    integrateaux <- function( lambda ) {
        xx <- exp(.vectorized.lk.optim(lambda, YksX, lambdasX, hX))
        yy <- exp(.vectorized.lk.optim(lambda, YksY, lambdasY, hY))
        simetricDivergenceW(  xx / yy, alpha)
    }
    lambdas <- c(max(min(c(lambdasX, lambdasY))), min(max(c(lambdasX, lambdasY))))
    a <- 0
    if (n > 0) {
        a <- multidiss.interp.SPEC(list(x,y), n, interp.W.LK, integrate.divergenceW, alpha)
    } else {
        tryCatch( {
            a <- integrate(integrateaux, min(lambdas), max(lambdas))$value
        }, error = function (e) {
            warning("Failed approximation with window from plug.in method, increasing window...")
            hX <- 2*hX
            hY <- 2*hY
            a <- integrate(integrateaux, min(lambdas), max(lambdas))$value
        })
    }
    if (plot) {
        plotsmoothspec(lambdasX, YksX, lambdasY, YksY, .vectorized.lk.optim, hX, hY, n)
    }
    a
}

distance.W.DLS <- function(x, y, alpha, plot=FALSE, n=length(x)) {
    .ts.sanity.check(x, y)
    .check.equal.length.ts(x,y) 
    YksX <- (spec.pgram(x,plot=FALSE)$spec)
    lambdasX <- spec.pgram(x,plot=FALSE)$freq
    YksY <- (spec.pgram(y,plot=FALSE)$spec)
    lambdasY <- (spec.pgram(y,plot=FALSE)$freq)
    hX <- dpill(lambdasX, YksX)
    hY <- dpill(lambdasY, YksY)
    
    integrateaux <- function( lambda ) {
        xx <- leastsquares.spec(YksX, lambdasX, hX, lambda )
        yy <- leastsquares.spec(YksY, lambdasY, hY, lambda )
        xx[xx<0.0001] <- 0.0001
        yy[yy<0.0001] <- 0.0001
        simetricDivergenceW(  xx / yy, alpha)
    }
    lambdas <- c(max(min(c(lambdasX, lambdasY))), min(max(c(lambdasX, lambdasY))))
    a <- 0
    if (n > 0) {
        a <- multidiss.interp.SPEC(list(x,y), n, interp.SPEC.LS, integrate.divergenceW, alpha)
    } else {
        a <- integrate(integrateaux, min(lambdas), max(lambdas), subdivisions=100)$value
    }
    if (plot) {
        plotsmoothspec(lambdasX, YksX, lambdasY, YksY, leastsquares.spec, hX, hY, n, "Least Squares")
    }
    
    a
}


diss.SPEC.LLR <- function(x,y, alpha=0.5, method="DLS", plot=FALSE, n=length(x)) {
    .ts.sanity.check(x, y)
    .check.equal.length.ts(x,y)
    typedist <- 0
    type <-  (pmatch(method, c("DLS", "LK" )))
    if (is.na(type)) {
        stop(paste("Unknown method", method))
    } else if (type == 1) {
        typedist <- distance.W.DLS(x,y, alpha, plot, n)
    }
    else if (type == 2) {
        typedist <- distance.W.LK(x,y, alpha, plot, n)
    }
    typedist
}

multidiss.SPEC.LLR <- function(series, method="DLS", alpha=0.5, plot=FALSE, n=length(series[[1]])) {
    if ( n > 0) {
        interpfun <- NULL
        type <-  (pmatch(method, c("DLS", "LK" )))
        if (is.na(type)) {
            stop(paste("Unknown method", method))
        } else if (type == 1) {
            interpfun <- interp.SPEC.LS
        }
        else if (type == 2) {
            interpfun <- interp.W.LK
        }
        multidiss.interp.SPEC(series, n, interpfun, integrate.divergenceW, alpha)
    } else {
        pairwise.diss( series, noindicesdiss(diss.SPEC.LLR), alpha, method, plot, n)
    }
}
multidiss.SPEC.GLK <- function(series, plot=FALSE) {
        multidiss.interp.SPEC(series, floor(length(series[[1]])/2), interp.SPEC.GLK, integrate.GLK)
}
multidiss.SPEC.ISD<- function(series, plot=FALSE,  n=length(series[[1]])) {
    if (n > 0) {
        multidiss.interp.SPEC(series, n, interp.SPEC.LOGLIKELIHOOD , integrate.ISD)
    } else {
        pairwise.diss( series, noindicesdiss(diss.SPEC.ISD), plot, n)
    }
}

    
diss.SPEC.GLK <- function(x,y, plot=FALSE ) {
    .ts.sanity.check(x, y)
    .check.equal.length.ts(x,y)
        
    interpx <- interp.SPEC.GLK(x, length(x)/2) #the value n is ignored only used for compatibility with
    interpy <- interp.SPEC.GLK(y, length(y)/2) #multidiss.SPEC
    
    if (plot) {
        YksX <- log(spec.pgram(x,plot=FALSE)$spec)
        lambdasX <- spec.pgram(x,plot=FALSE)$freq
        YksY <- log(spec.pgram(y,plot=FALSE)$spec)
        lambdasY <- (spec.pgram(y,plot=FALSE)$freq)
        hX <- 0.93*dpill(lambdasX, YksX)
        hY <- 0.93*dpill(lambdasY, YksY)
        plotsmoothspec(lambdasX, YksX, lambdasY, YksY, .vectorized.lk.optim, hX, hY, length(x)/2)
    }
    #distance GLK
    integrate.GLK( interpx$x, interpx$y, interpy$y)
}



#distancia ISD
diss.SPEC.ISD <- function(x,y, plot=FALSE, n=length(x)) {
    .ts.sanity.check(x, y)
    .check.equal.length.ts(x,y)
 
    
    YksX <- log(spec.pgram(x,plot=FALSE)$spec)
    lambdasX <- spec.pgram(x,plot=FALSE)$freq
    YksY <- log(spec.pgram(y,plot=FALSE)$spec)
    lambdasY <- (spec.pgram(y,plot=FALSE)$freq)
    
    hX <- 0.93*dpill(lambdasX, YksX)
    hY <- 0.93*dpill(lambdasY, YksY)
    
    integraISDaux <- function(lambda) {
        (.vectorized.lk.optim(lambda, YksX, lambdasX  , hX) - .vectorized.lk.optim(lambda,YksY, lambdasY,  hY))**2
    }
    a <- 0
    if (n > 0) {
        a <- multidiss.interp.SPEC(list(x,y), n, interp.SPEC.LOGLIKELIHOOD , integrate.ISD)
    } else {
        tryCatch( {
            a <- integrate(integraISDaux, min(lambdasX), max(lambdasX) )$value
        }, error = function(e) {
            hX <- 2*hX
            hY <- 2*hY
            a <- integrate(integraISDaux, min(lambdasX), max(lambdasX) )$value
        })
    }
    if (plot) {
        plotsmoothspec(lambdasX, YksX, lambdasY, YksY, .vectorized.lk.optim, hX, hY, n)
    }  
    a
}


#########################################################
##############  distance CEPSTRAL  ######################
#########################################################

#combines non seasonal and seasonal autregressive coefficients in a single non-seasonal coef vector
#coef+period to plain coef, and the mutliplication of nonseasonal and seasonal parts
.seasontoplain = function(ar,sar,period) {
  #initialization to common length vectors
  maxlen = max(length(ar), length(sar)*period)
  arcoef = rep.int(0, maxlen)
  sarcoef = rep.int(0, maxlen)
  
  arcoef[1:length(ar)] = ar
  #transform coefficients + period to plain vector
  sarcoef[(0:length(sar)*period)] = sar
  
  #multiply the coefficients (sort of polynomial multiplication)
  endcoef = rep(0,2*maxlen)
  for (i in 1:maxlen) {
    for (j in 1:maxlen) {
      endcoef[i + j] = endcoef[i + j] - arcoef[i]*sarcoef[j]
    }
  }
  endcoef[1:length(arcoef)] = endcoef[1:length(arcoef)] + arcoef
  endcoef[1:length(sarcoef)] = endcoef[1:length(sarcoef)] + sarcoef
  endcoef
}

.arma2ar <- function(ar, ma, lag) {
  coef<- ARMAacf(ar, ma, lag)
  acf2AR(coef)[lag,]
}

.calc.cepstral.coef <- function( ARx, h ) {
  CEPSTRALx <- 1:h
  CEPSTRALx[1] <- ARx[1]   
  if (length(ARx) >= 2) {
    for (i in 2:length(ARx)) {
      acum <- 0
      for (m in 1:(i-1) ) {
        acum <- acum + ( 1- m/i)*ARx[m]*CEPSTRALx[i-m]
      }
      CEPSTRALx[i] <- ARx[i] + acum   
    }
  }
  
  
  if (h > length(ARx)) {
    for (i in (length(ARx)+1):h) {
      acum <- 0
      for (m in 1:length(ARx)) {
        acum <- acum + (1 - m/i)*ARx[m]*CEPSTRALx[i-m]
      }
      CEPSTRALx[i] <- acum        
    }
  }
  CEPSTRALx
}


cepstral <- function(x, h, order=NULL, seasonal, permissive) {
    ARx <- NULL
    
    SAR <- 0
    period <- 0
    if (is.null(order)) { #if order null, fit automatically
        ARx <- ar(x,order.max=min(length(x)-1,h))
    } else {
        if ((order[1]) < 1) stop("The arima order must have AR coefficients, they are used for the distance")
        arim <- arima(x, order, seasonal) #calc the ARIMA coefs from the imposed model
        #convert ARMA to AR
        arcoef = arim$coef[1:arim$arma[1]][0:arim$arma[1]]  #get the coefs
        macoef = arim$coef[(1+arim$arma[1]):(1+arim$arma[1]+arim$arma[2] -1)][0:arim$arma[2]]
        ARx$ar <- .arma2ar(ar=arcoef, ma = macoef, 50) #convert

        #the seasonal part ARMA to AR
        Scoef <- arim$coef[ (1 + arim$arma[1] + arim$arma[2]):(1 + arim$arma[1] + arim$arma[2]+ arim$arma[3] + arim$arma[4]) ]
        SAR <- Scoef[0:arim$arma[3]]
        SMA <- Scoef[ - (0:arim$arma[3])]
	if (length(SAR) == 0) {
	  SAR = 0
	}
	if (length(SMA) == 0) {
	  SMA = 0
	}

        SAR <- .arma2ar(ar=SAR, ma = SMA, 50)
        period <- arim$arma[5]
    }
    
    if (length(ARx$ar) < 1) {
        if (permissive) {
            warning("Cepstral distance, error on the selection of the AR order, 0 by AIC, forcing 1")
            ARx  <- ar(x, aic=FALSE, order.max=1)
        }
        if (length(ARx)<1) {
            stop("Could not find any AR coefficient")      
        }
    }
    
    coef <- .seasontoplain(ARx$ar, SAR, period)
    
    .calc.cepstral.coef(coef, h)
    
}


diss.AR.LPC.CEPS <- function(x, y, k=50, order.x=NULL, order.y= NULL,
                    seasonal.x = list(order = c(0, 0, 0), period = NA),
                    seasonal.y = list(order = c(0, 0, 0), period = NA),
                    permissive=TRUE) {
    .ts.sanity.check(x, y)
    if (!is.list(seasonal.x) || !is.list(seasonal.y)) {
            stop("Invalid seasonal part")
    }
    cpx <- cepstral(x, k, order.x, seasonal.x, permissive)
    cpy <- cepstral(y, k, order.y, seasonal.y, permissive)
    as.numeric(dist(rbind(cpx,cpy)))
}


#############################################################################
#################   Temporal Correlation Distance   #########################
#############################################################################

##CHOUAKRIA-DOUZAL

corrtemporder1 <- function (x, y) {
    p <- length(x)
    sum((x[2:p] - x[1:(p-1)]) * (y[2:p] - y[1:(p-1)])) / ( sqrt( sum((x[2:p] - x[1:(p-1)])^2) ) * sqrt( sum((y[2:p] - y[1:(p-1)])^2) ))
}

diss.CORT <- function( x, y, k=2, deltamethod="Euclid") {
    .ts.sanity.check(x, y)
    .check.equal.length.ts(x,y)
    corrt <- corrtemporder1(x,y)
    type <-  (pmatch(deltamethod, c("Euclid", "Frechet", "DTW")))
    typedist <- 0
    if (is.na(type)) {
        stop(paste("Unknown method", deltamethod))
    } else if (type == 1) {
        typedist <- as.numeric( dist(rbind(x,y)) )
    }
    else if (type == 2) {
        typedist <- diss.FRECHET(x,y)
    }
    else if (type == 3) {
        typedist <- dtw(x,y, dist.method="Manhattan", distance.only=T)$distance
    }
    
    (2/( 1+ exp(k*corrt)))*typedist
    
}



##################################################
######  maharaj clustering algorithm #############
##################################################
#input, distance matrix
pvalues.clust <- function(pvalues,significance) {
    distancias <- pvalues
    significacion <- significance
    distancias <- as.matrix(distancias)
    tam <- dim(distancias)[1]
    distancias
    combinaciones <- combn(tam,2)
    plandist <- 1:ncol(combinaciones)
    #create a vector with the distances
    for (i in 1:length(plandist)) {
        plandist[i] <- distancias[combinaciones[1,i], combinaciones[2,i]]
    }
    
    ord <- order(plandist, decreasing=TRUE)
    plandist <- plandist[ord]
    combinaciones <- combinaciones[,ord]
    
    is_in_setlist <- function( element, setlist) {
        for (i in 1:length(setlist)) {
            if (element %in% setlist[[i]]) {
                return(TRUE)
            }
        }
        return(FALSE)
    }
    
    find_which_set <- function ( element, setlist) {
        for (i in 1:length(setlist)) {
            if (element %in% setlist[[i]]) {
                return(i)
            }
        }
        return(0)
    }
    
    
    is_pvalue_of_element_less_than_significance_with_any_in_set <- function( element, setid, setlist, significance,distances) {
        sum(distances[element, setlist[[setid]]] < significance) > 0
    }
    
    add_to_set <- function ( element, setindex, setlist) {
        setlist[[setindex]] <- union(setlist[[setindex]], element)
        setlist
    }
    
    create_new_set <- function( element, setlist) {
        setlist[[length(setlist) + 1]] <- element
        setlist
    }
    
    is_all_series_already_in_a_cluster <- function(series, setlist) {
        already = TRUE
        for (i in series) {
            already <- already & (find_which_set(i,setlist) != 0)
        }
        already
    } 
    
    are_pvalues_of_all_pairs_across_clusters_greater_than_significance <- function( clusterone, clustertwo, setlist, significance, distances) {
        seriesone <- setlist[[clusterone]]
        seriestwo <- setlist[[clustertwo]]
        aregreater <- TRUE
        for (i in seriesone) {
            aregreater <- aregreater & (sum(( distances[i, seriestwo] > significance)) == length(seriestwo) )
        }
        for (i in seriestwo) {
            aregreater <- aregreater & (sum(( distances[i, seriesone] > significance)) == length(seriesone) )
        }
        aregreater
    }
    
    merge_sets <- function (setone, settwo, setlist) {
        setlist[[settwo]] <- c( setlist[[settwo]] ,setlist[[setone]])
        setlist[-setone]
    }
    
    if (plandist[1] < significacion) {
        grupos <- as.list(1:tam)
    } else {
        grupos <- list()
        grupos[[1]] <- combinaciones[,1]
        for (i in 2:length(combinaciones[1,])) {
            if (plandist[i] < significacion) { #p-value < significance = YES
                conj <- find_which_set( combinaciones[1,i], grupos )
                if (conj == 0) {
                    grupos <- create_new_set(combinaciones[1,i], grupos)
                }
                conj <- find_which_set( combinaciones[2,i], grupos )
                if (conj == 0) {
                    grupos <- create_new_set(combinaciones[2,i], grupos)
                } 
                #each remaining serie to its own cluster
                for (j in i:length(combinaciones[1,])) {
                    conj <- find_which_set( combinaciones[1,j], grupos )
                    if (conj == 0) {
                        grupos <- create_new_set(combinaciones[1,j], grupos)
                    }
                    conj <- find_which_set( combinaciones[2,j], grupos )
                    if (conj == 0) {
                        grupos <- create_new_set(combinaciones[2,j], grupos)
                    } 
                }
                break;
            } else { #p-value < significance = NO
                conj <- is_in_setlist( combinaciones[1,i], grupos )
                conj <- conj + is_in_setlist( combinaciones[2,i], grupos )
                if (conj < 2) {  #is each(ALL) series already in a cluster = NO
                    conj <- find_which_set( combinaciones[1,i], grupos )
                    if (conj > 0) { #one of the series in the pair already in a cluster = YES (x)
                        if (is_pvalue_of_element_less_than_significance_with_any_in_set( combinaciones[2,i], conj, grupos, significacion, distancias  )    ) {
                            grupos <- create_new_set(combinaciones[2,i], grupos)
                        } else {
                            grupos <- add_to_set( combinaciones[2,i], conj, grupos )
                        }
                    } else {
                        conj <- find_which_set( combinaciones[2,i], grupos )
                        if (conj > 0) {#one of the series in the pair already in a cluster = YES (y)
                            if (is_pvalue_of_element_less_than_significance_with_any_in_set( combinaciones[1,i], conj, grupos, significacion, distancias  )    ) {
                                grupos <- create_new_set(combinaciones[1,i], grupos)
                            } else {
                               grupos <- add_to_set( combinaciones[1,i], conj, grupos )
                            }
                        } else {#one of the series in the pair already in a cluster = NO
                            grupos <- create_new_set( combinaciones[1,i], grupos ) # create a new cluster with the two series
                            setid <- find_which_set( combinaciones[1,i], grupos )
                            grupos <- add_to_set( combinaciones[2,i], setid, grupos )
                        }
                    } 
                } else {  #is each series already in a cluster = YES
                    conj1 <- find_which_set( combinaciones[1,i], grupos )
                    conj2 <- find_which_set( combinaciones[2,i], grupos )
                    if (are_pvalues_of_all_pairs_across_clusters_greater_than_significance(conj1, conj2, grupos, significacion, distancias)) {
                        #merge
                        grupos <- merge_sets( conj1, conj2, grupos)
                    } else {
                        #do nothing
                    }
                }
            }
        }  
    }
    
    result <- 1:tam
    for ( i in 1:length(grupos) ) {
        for (j in grupos[[i]]) {
            result[j] <- i
        }
    }
    result
}


#check if the maharaj algorithm is equivalent to hierarchical clustering with 
# a cut level
testIgualdadMaharajHCLUST <- function( distancias, pvalor) {
    clusterp <- pvalues.clust(as.dist(distancias), pvalor)
    clusterh <- (hclust(as.dist(-distancias), method="single"))
    print(clusterp)
    plot(clusterh)
}


####################################################################
############# Feature Extraction Based on Wavelet Transform ########
####################################################################

wavelet.feature.extraction <- function(series) {
    
    calcEnergies <- function(wavdecomp) {
        level <- length(wavdecomp$data) -1
        energyD <- rep.int(0,level)
        energyD[1] <- sum(wavdecomp$data[[1]]**2)
        for ( i in 1:level) {
            energyD[i] <- sum((wavdecomp$data[[i]])**2)
        }
        return (energyD)
    }
    
    #fill with zeroes
    max.level <- as.integer(ceiling(logb(length(series[1,]),base=2)))
    true.level <- as.integer(floor(logb(length(series[1,]),base=2)))
    if (max.level != true.level) {
        npad <- 2**max.level - length(series[1,])
        for (i in 1:npad) {
            series <- cbind(series, rep(0, nrow(series)))
        }
    }
    
    energies <- matrix(0, nrow=nrow(series), ncol = max.level)
    for ( i in 1:nrow(series) ) {
        wavdecomp <- wmtsa::wavDWT(series[i,], n.levels=max.level, wavelet="haar")
        energies[i,] <- calcEnergies(wavdecomp)
    }
    
    sumEnergies <- colSums(energies)
    
    final_level <- max.level
    for (i in 1:(max.level-1)) {
        if (sumEnergies[i] < sumEnergies[i+1]) {
            final_level <- i
            break
        }
    }
    wavdecomp <- wavDWT(series[1,], n.levels=final_level, wavelet="haar")
    out.series <- wavdecomp$data[[final_level+1]]
    for (i in 2:nrow(series)) {
        wavdecomp <- wavDWT(series[i,], n.levels=final_level, wavelet="haar")
        out.series <- rbind(out.series, wavdecomp$data[[final_level+1]])
    }
    out.series
}

diss.DWT <- function(series) {
    if ( length(dim(series)) == 2 ) {
        if ( dim(series)[1] < 2 ) {
            stop( "diss.DWT needs a minimum of 2 series to compute the distance, incorrect amount provided" )
        }
    } else {
        stop( "diss.DWT series matrix with incorrect dimensions" ) 
    }
    wt <- wavelet.feature.extraction( series )
    dist(wt)
}

##########################################################
############  CORRELATION BASED DISTANCES ################
##########################################################



diss.COR <- function(x, y, beta = NULL) {
    .ts.sanity.check(x, y)
    correl <- cor(x,y)
    if (is.null(beta)) {
        sqrt(2*(1- correl))
    } else {
        if (beta<0) {
            stop("beta must be greater than 0")
        }
        sqrt( ((1-correl)/(1+correl ))**beta )    
    }
}


####################################################################################
########################## CDM KEOGH 2004 ##########################################
################## Compression based data mining of sequential data ################
####################################################################################

#common part of compression methods,
#calculate the sizes of the compressed series and of their concatenation
.compression.lengths <- function(x, y, type) {      
    methods <- type
    type = match.arg(type, c("gzip", "bzip2", "xz", "min"))
    if (type == "min") { #choose the best compression method of the three 
        methods <- c("gzip", "bzip2", "xz")
    }
    xy <- as.character(c(x,y))
    x <- as.character(x)
    y <- as.character(y)
    cxym <- sapply( methods, function(m) { length( memCompress(xy, type=m) )})
    cxy <- min(cxym)
    cx <- min(sapply( methods, function(m) { length( memCompress(x, type=m) )}))
    cy <- min(sapply( methods, function(m) { length( memCompress(y, type=m) )}))
    list(cx=cx, cy=cy, cxy=cxy)
}

#length of the compressed concatenated series / sum lengths of the compressed series
diss.CDM <- function(x, y, type="min") {
    .ts.sanity.check(x, y)
    comp <- .compression.lengths(x,y, type)    
    comp$cxy / (comp$cx + comp$cy)
}


###################################################################################
####### Clustering by Compression (2005), Cilibrasi, R., Vitanyi, P.M.B.,  ########
######## Normalized Compression Distance ##########################################
###################################################################################
diss.NCD <- function(x,y, type="min") {
    .ts.sanity.check(x,y)
    comp <- .compression.lengths(x,y, type)  
    (comp$cxy - min(comp$cx,comp$cy)) / max(comp$cx, comp$cy)
}

##############################################################
################# COMPLEXITY BASED DISTANCE ##################
#BatistaWangKeogh_2011_SDM_A complexity invariant distance measure for TS#
##############################################################

diss.CID = function(x, y) {
    .ts.sanity.check(x, y)
    .check.equal.length.ts(x,y)
    CED.x <- sqrt( sum( diff(x)^2) ) #complexities of the series
    CED.y <- sqrt( sum( diff(y)^2) )
    denom <- min(CED.x, CED.y)
    if(denom == 0) {
      stop("Cannot divide by zero: A series exists that has complexity zero.")
    }
    CF = max(CED.x, CED.y) / denom #complexity correction factor
    CF * dist(rbind(x,y))
}



#######################################################################
############### CLUSTER SIMILARITY INDEX ##############################
#######################################################################

#gravrilov similarity ratio, used in kalpakis
Sim <- function(Gi, Sj, G, S) {
    2* (sum ( (G==Gi) & (S==Sj) ) ) / ( sum(G==Gi) + sum(S==Sj))
}

cluster.evaluation <- function(G,S) {
    if (length(G) != length(S)) {
        stop("Different amount of elements between cluster solutions")
    }
    if (any(is.na(G)) || any(is.na(S))) {
        stop("NA in the cluster solutions")
    }
    acum <- 0
    gclust <- unique(G)
    sclust <- unique(S)
    for (i in gclust) {
        maxS <- 0
        for (j in sclust) {
            ms <- Sim(i, j, G, S)
            if (ms > maxS) {
                maxS <- ms
            }
        }
        acum <- acum + maxS
    }
    acum/(length(gclust))
}


###################################################################
########################### DISS WRAPPER ##########################
###################################################################

#series is a list with time series
#dissfun is a function that takes the list of series and their indices, and extra paramters
#... parameters for dissfun
pairwise.diss <- function(series, dissfun, ...) {   
    n <- length(series)
    distances <- matrix(0, n, n)
    for (i in 1:(n-1)) {
        for (j in (i+1):n) {
            tryCatch( {
                d <- dissfun( series, i, j, ...)
                distances[i,j] <- d
                distances[j,i] <- d
            }, error = function (e) {
                stop( paste("Applying diss, series (",i,",",j,") produced the following error: ", e) )
            })
        }
    }
    as.dist((distances))
}

noindicesdiss <- function( fun ) {
    function(series, i, j, ...) {
        fun(series[[i]], series[[j]], ...)
    }
}



diss <- function(SERIES, METHOD, ...) {
  pair.diss.fun = pairwise.diss
    if (!is.matrix(SERIES) && !is.list(SERIES) && !is.mts(SERIES)) {
        stop("list, mts, matrix or data.frame object is required for SERIES ")
    }
    mat.ser <- SERIES
    if (is.mts(SERIES)) {
        SERIES <- t( as.matrix(SERIES))
    }
    
    if (!is.list(SERIES)) {
        tmpser <- SERIES
        SERIES <- list()
        for (i in 1:nrow(tmpser)) {
            SERIES[[i]] <- tmpser[i,]
        }
        names(SERIES) <- rownames(tmpser)
    }
    
    #common check for input parameters
    if (any(is.na(unlist(SERIES)))) {
        stop("NA in the series")
    }
    if (length(SERIES) < 2) {
        stop("Only one series provided")
    }
    
    list.to.matrix <- function(series) {
        n <- length(series)
        k <- length(series[[n]])
        mat.ser <- matrix(0, n, k)
        for (i in 1:n) {
            if ( length( series[[i]]) != k ) {
                stop(paste("diss method",METHOD,"requires same length series"))
            }
            mat.ser[i,] <- series[[i]]
        }
        rownames(mat.ser) <- names(series)
        mat.ser
    }

    out.dist <- NULL
    
    METHODS = c("ACF", "PACF", "AR.MAH", "AR.PIC", "AR.LPC.CEPS", "PER", "INT.PER", "COR", "CORT", "DWT",
                "PDC", "PRED", "MINDIST.SAX", "SPEC.LLR", "SPEC.GLK", "SPEC.ISD", "CDM", "CID", "NCD", "DTWARP", "FRECHET",
                "EUCL")
    diss.method = match.arg(METHOD, METHODS)
    #get the statistic of the MAHARAJ dissimilarity
    diss.AR.MAH.STAT <- function(x,y, ...) {
        diss.AR.MAH(x,y,...)$statistic
    }
    diss.AR.MAH.PVAL <- function(x,y, ...) {
        diss.AR.MAH(x,y,...)$p_value
    }
    diss.fun <- switch(diss.method,
        ACF = diss.ACF,
        PACF = diss.PACF,
        AR.MAH = diss.AR.MAH,
        AR.PIC = diss.AR.PIC,
        AR.LPC.CEPS = diss.AR.LPC.CEPS,
        PER = diss.PER,
        INT.PER = diss.INT.PER,
        COR = diss.COR,
        CORT = diss.CORT,
        DWT = diss.DWT,
        PRED = multidiss.PRED,
        SPEC.LLR = diss.SPEC.LLR,
        SPEC.GLK = diss.SPEC.GLK,
        SPEC.SD = diss.SPEC.ISD,
        MINDIST.SAX = diss.MINDIST.SAX,
        PDC = pdcDist,
        CDM = diss.CDM,
        CID = diss.CID,
        NCD = diss.NCD,
        DTWARP = diss.DTWARP,
        FRECHET = diss.FRECHET,
        EUCL = diss.EUCL)
    
    if (diss.method == "DWT") { #diss dwt is not a pairwise diss, we cannot use proxy::dist
        out.dist <- diss.DWT( list.to.matrix(SERIES) )
    } else if (diss.method == "AR.PIC") {
        multi.PIC <- function(SERIES, i, j, order=NULL, permissive=TRUE, order.x=NULL, order.y=NULL) {
            if (!is.null(order.x) || !is.null(order.y) ) {
                stop("AR.PIC from the diss wrapper function must be called using 'order' argument, not with
                     order.x and order.y, see diss.AR.PIC help page")
            }
             diss.AR.PIC(SERIES[[i]], SERIES[[j]], order[i,], order[j,], permissive)
        }
        out.dist <- pair.diss.fun(SERIES, multi.PIC, ...)
    } else if (diss.method == "AR.LPC.CEPS") {
        multi.CEPS <- function(series, i, j, k=50, order=NULL, seasonal=NULL, permissive=TRUE,
                        order.x=NULL, order.y=NULL, seasonal.x=NULL, seasonal.y=NULL) { #arguments to inform incorrect usage
            if (!is.null(order.x) || !is.null(order.y) || !is.null(seasonal.x) || !is.null(seasonal.y)) {
                stop("AR.LPC.CEPS from the diss wrapper function must be called using 'order' and 'seasonal' arguments, not with
                     order.x, order.y, seasonal.x or seasonal.y arguments, see diss.AR.LPC.CEPS help page")
            }
            if (is.null(seasonal)) {
                seasonal[[i]] <-  list(order=c(0,0,0), period=NA)
                seasonal[[j]] <-  list(order=c(0,0,0), period=NA)
            }
            distance <- diss.AR.LPC.CEPS(series[[i]], series[[j]], k, order[i,], order[j,], seasonal[[i]], seasonal[[j]] )
        }
        out.dist <- pair.diss.fun(SERIES, multi.CEPS, ...)
    } else if (diss.method == "PRED") {
        return( multidiss.PRED(SERIES, ...) ) #TODO proper names
    } else if (diss.method == "AR.MAH") {

        statistic = pair.diss.fun( SERIES, noindicesdiss(diss.AR.MAH.STAT), ...)
        p_value = pair.diss.fun( SERIES, noindicesdiss(diss.AR.MAH.PVAL), ...)
        return( list(statistic=statistic, p_value=p_value) ) #TODO proper naming of the output
    } else if (diss.method == "PDC") {
        out.dist <- pdcDist( t( list.to.matrix(SERIES) ), ...)
    } else if (diss.method == "SPEC.LLR") { #for performance reasons, we must call these in a different way
        out.dist <- multidiss.SPEC.LLR(SERIES, ...)
    } else if (diss.method == "SPEC.GLK") { #for performance reasons, we must call these in a different way
        out.dist <- multidiss.SPEC.GLK(SERIES, ...)
    } else if (diss.method == "SPEC.ISD") { #for performance reasons, we must call these in a different way
        out.dist <- multidiss.SPEC.ISD(SERIES, ...)
    } else {
        out.dist <- pair.diss.fun( SERIES, noindicesdiss(diss.fun), ...)
    }
    out.dist <- as.dist(out.dist)
    names(out.dist) <- names(SERIES)
    out.dist
}

loo1nn.cv <- function(d, G) {
    d <- as.matrix(d)
    diag(d) <- max(d) + 1 #distance with self is not included
    nearest.set <- apply ( d, 1, function(x) { #handle ties
                                    which( x == min(x) )
                                }
                         )
    
    #mode auxiliary function
    #based on http://stackoverflow.com/questions/2547402/standard-library-function-in-r-for-finding-the-mode
    #with random tie breaking
    my.mode <- function(x) {
        ux <- unique(G[x])
        if (length(ux) > 1) {
            ux <- sample(ux)
        }
        prop <- tabulate(match(G[x], ux))
        if (sum(prop == max(prop))>1) {
            warning("There were ties on the voting, selecting one at random")
        }
        ux[which.max(prop)]
    }
    
    #apply the mode, depending on whether apply obtained different amount of ties (list), the same amount(matrix) or no ties(vector) 
    nearest <- switch( class(nearest.set), matrix = apply(nearest.set, 2, my.mode),
                       numeric = apply(as.matrix(nearest.set), 1, my.mode),
                       integer = apply(as.matrix(nearest.set), 1, my.mode),
                       list = sapply(nearest.set, my.mode)
    )
    sum(nearest == G)/length(G)
}

#wrappers for easier discovery of the available functions
diss.EUCL <- function(x, y) {
    dist(rbind(x,y))
}
diss.DTWARP <- function(x,y,...) {
    dtw(x,y, ...)$distance
}
diss.FRECHET <- function(x,y,...) {
    abscissex = 1:length(x)
    abscissey = 1:length(y)
    distFrechet(abscissex,x,abscissey,y, ...)
}
diss.PDC <- function(x,y, ...) {
    pdcDist(cbind(x,y), ...)
}
############################################################################
#######################   OLD STUFF (UNUSED)   #############################
############################################################################



#distancia basada en correlaciones cruzadas
distanciaCorCruLagK = function(x,y,k) {
    muX = mean(x)
    muY = mean(y)
    TT = length(x)
    sum( (x[1:(TT-k)] - muX)*( y[(1+k):TT] - muY)) / (sqrt( sum((x - muX)**2))*sqrt( sum((y - muY)**2)))
}

#distanciaCorCruLagK(x,y,18)

distanciaCorCruTotal = function(x,y) {
    k = length(x)-1
    denom = 0
    for (i in 1:k) {
        denom = denom + distanciaCorCruLagK(x,y,i)
        print(denom)
    }
    
    sqrt( (1 - distanciaCorCruLagK(x,y,0))/denom  )
    denom
}

#distanciaCorCruTotal(x,y)

Try the TSclust package in your browser

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

TSclust documentation built on Nov. 17, 2017, 7:24 a.m.