R/alignment.R

In heike/x3prplus: Algorithms for Matching Bullet Lands

#' @importFrom stats cor
bulletAlign_nist <- function (data, value = "l30", mincor = .8)  {
    bullet <- NULL
    b12 <- unique(data$bullet)
    if (length(b12) != 2) 
        stop("Two surfaces should be compared\n\n")
    data$val <- data.frame(data)[, value]
    
    subLOFx1 <- subset(data, bullet == b12[1])
    subLOFx2 <- subset(data, bullet == b12[2])
    subLOFx1$y <- subLOFx1$y - min(subLOFx1$y)
    subLOFx2$y <- subLOFx2$y - min(subLOFx2$y)
    
    whichmin <- which.min(c(length(subLOFx1$y), length(subLOFx2$y)))
    shorter <- list(subLOFx1$val, subLOFx2$val)[[whichmin]]
    longer <- list(subLOFx1$val, subLOFx2$val)[[3 - whichmin]]
    
    longer_na <- c(rep(0, length(shorter)), longer, rep(0, (length(shorter))))
    
    mycors <- NULL
    for (i in 1:(length(longer_na) - length(shorter))) {
        longersub <- longer_na[i:(i + length(shorter) - 1)]
        
        corval <- cor(shorter, longersub, use = "pairwise.complete.obs")
        
        mycors <- c(mycors, corval)
    }
    
    lag <- which.max(mycors) - length(shorter)
    if (max(mycors, na.rm = TRUE) < mincor) lag <- 0
    
    incr <- min(diff(sort(unique(subLOFx1$y))))
    
    mydat <- if (whichmin == 1) subLOFx1 else subLOFx2
    mydat2 <- if (whichmin == 1) subLOFx2 else subLOFx1
    
    if (lag < 0) {
        mydat2$y <- mydat2$y + lag * incr
    } else {
        mydat$y <- mydat$y + lag * incr
    }
    
    bullets <- rbind(data.frame(mydat), data.frame(mydat2))
    list(ccf = max(mycors, na.rm = TRUE), lag = lag * incr, bullets = bullets)
}

#' Align two surface cross cuts according to maximal correlation
#' 
#' The bullet with the first name serves as a reference, the second bullet is shifted.
#' @param data data frame consisting of at least two surface crosscuts as given by function \code{bulletSmooth}.
#' @param value string of the variable to match. Defaults to l30, the variable returned from function \code{bulletSmooth}.
#' @return list consisting of a) the maximal cross correlation, b) the lag resulting in the highest cross correlation, and c) same data frame as input, but y vectors are aligned for maximal correlation between the 
#' @export
#' @importFrom stats na.omit
bulletAlign <- function(data, value = "l30") {
    bullet <- NULL
    b12 <- unique(data$bullet)
    
    if (length(b12) != 2) stop("Two surfaces should be compared\n\n")
    
    data$val <- data.frame(data)[, value]
    miny <- min(data$y, na.rm=T)
    
    subLOFx1 <- subset(data, bullet==b12[1])
    subLOFx2 <- subset(data, bullet==b12[2]) 
    
    # ccf assumes that both time series start at the same y
    # shift series into the same origin
    subLOFx1$y <- subLOFx1$y - min(subLOFx1$y)
    subLOFx2$y <- subLOFx2$y - min(subLOFx2$y)
    
    ccf <- ccf(subLOFx1$val, subLOFx2$val, plot = FALSE, lag.max=1500, 
               na.action = na.omit)
    lag <- ccf$lag[which.max(ccf$acf)]
    incr <- min(diff(sort(unique(subLOFx1$y))))
    
    subLOFx2$y <- subLOFx2$y +  lag * incr # amount of shifting 
    bullets <- rbind(data.frame(subLOFx1), data.frame(subLOFx2))
    #  bullets$y <- bullets$y + miny # we can, but we don't have to shift the series back. This is rather cosmetic.
    list(ccf=max(ccf$acf), lag = lag * incr, bullets=bullets)
}