R/ffcsaps.R

In dplR: Dendrochronology Program Library in R

ffcsaps <- function(y, x=seq_along(y), nyrs=length(y)/2, f=0.5) {
### support functions
    ffppual <- function(breaks, c1, c2, c3, c4, x, left){
        if (left){
            ix <- order(x)
            x2 <- x[ix]
        } else{
            x2 <- x
        }

        n.breaks <- length(breaks)
        if (left) {
            ## index[i] is maximum of a and b:
            ## a) number of elements in 'breaks[-n.breaks]' that are
            ##    less than or equal to x2[i],
            ## b) 1
            index <- pmax(ffsorted(breaks[-n.breaks], x2), 1)
        } else {
            ## index[i] is:
            ## 1 + number of elements in 'breaks[-1]' that are
            ## less than x2[i]
            index <- ffsorted2(breaks[-1], x2)
        }

        x2 <- x2 - breaks[index]
        v <- x2 * (x2 * (x2 * c1[index] + c2[index]) + c3[index]) + c4[index]

        if (left)
            v[ix] <- v
        v
    }

    ffsorted <- function(meshsites, sites) {
        index <- order(c(meshsites, sites))
        which(index > length(meshsites)) - seq_along(sites)
    }

    ffsorted2 <- function(meshsites, sites) {
        index <- order(c(sites, meshsites))
        which(index <= length(sites)) - seq(from=0, to=length(sites)-1)
    }

    ## Creates a sparse matrix A of size n x n.
    ## The columns of B are set to the diagonals of A so that column k
    ## becomes the diagonal in position d[k] relative to the main
    ## diagonal (zero d[k] is the main diagonal, positive d[k] is
    ## above, negative is below the main diagonal).
    ## A value on column j in A comes from row j in B.
    ## This is similar in function to spdiags(B, d, n, n) in MATLAB.
    spdiags <- function(B, d, n) {
        n.d <- length(d)
        A <- matrix(0, n.d * n, 3)
        count <- 0
        for(k in seq_len(n.d)){
            this.diag <- d[k]
            i <- inc(max(1, 1 - this.diag), min(n, n - this.diag)) # row
            n.i <- length(i)
            if(n.i > 0){
                j <- i + this.diag                                 # column
                row.idx <- seq(from=count+1, by=1, length.out=n.i)
                A[row.idx, 1] <- i
                A[row.idx, 2] <- j
                A[row.idx, 3] <- B[j, k]
                count <- count + n.i
            }
        }
        A <- A[A[, 3] != 0, , drop=FALSE]
        A[order(A[, 2], A[, 1]), , drop=FALSE]
    }

### start main function

    y2 <- as.numeric(y)
    ## If as.numeric() does not signal an error, it is unlikely that
    ## the result would not be numeric, but...
    if(!is.numeric(y2)) stop("'y' must be coercible to a numeric vector")
    x2 <- as.numeric(x)
    if(!is.numeric(x2)) stop("'x' must be coercible to a numeric vector")

    n <- length(x2)
    ## quick error check
    if (n < 3) stop("there must be at least 3 data points")
    if(!is.numeric(f) || length(f) != 1 || f < 0 || f > 1)
        stop("'f' must be a number between 0 and 1")
    if(!is.numeric(nyrs) || length(nyrs) != 1 || nyrs <= 1)
        stop("'nyrs' must be a number greater than 1")

    ix <- order(x2)
    zz1 <- n - 1
    xi <- x2[ix]
    zz2 <- n - 2
    diff.xi <- diff(xi)

    ## quick error check
    if (any(diff.xi == 0)) stop("the data abscissae must be distinct")

    yn <- length(y2)

    ## quick error check
    if (n != yn)
        stop("abscissa and ordinate vector must be of the same length")

    arg2 <- -1:1
    odx <- 1 / diff.xi
    R <- spdiags(cbind(c(diff.xi[-c(1, zz1)], 0),
                       2 * (diff.xi[-1] + diff.xi[-zz1]),
                       c(0, diff.xi[-c(1, 2)])),
                 arg2, zz2)
    R2 <- spdiags(cbind(c(odx[-zz1], 0, 0),
                        c(0, -(odx[-1] + odx[-zz1]), 0),
                        c(0, 0, odx[-1])),
                  arg2, n)
    R2[, 1] <- R2[, 1] - 1
    forR <- Matrix(0, zz2, zz2, sparse = TRUE)
    forR2 <- Matrix(0, zz2, n, sparse = TRUE)
    forR[R[, 1:2, drop=FALSE]] <- R[, 3]
    forR2[R2[, 1:2, drop=FALSE]] <- R2[, 3]
    ## The following order of operations was tested to be relatively
    ## accurate across a wide range of f and nyrs
    p.inv <- (1 - f) * (cos(2 * pi / nyrs) + 2) /
        (12 * (cos(2 * pi / nyrs) - 1) ^ 2) / f + 1
    yi <- y2[ix]
    p <- 1 / p.inv
    mplier <- 6 - 6 / p.inv # slightly more accurate than 6*(1-1/p.inv)
    ## forR*p is faster than forR/p.inv, and a quick test didn't
    ## show any difference in the final spline
    u <- as.numeric(solve(mplier * tcrossprod(forR2) + forR * p,
                          diff(diff(yi) / diff.xi)))
    yi <- yi - mplier * diff(c(0, diff(c(0, u, 0)) / diff.xi, 0))
    test0 <- xi[-c(1, n)]
    c3 <- c(0, u / p.inv, 0)
    x3 <- c(test0, seq(from=xi[1], to=xi[n], length = 101))
    cc.1 <- diff(c3) / diff.xi
    cc.2 <- 3 * c3[-n]
    cc.3 <- diff(yi) / diff.xi - diff.xi * (2 * c3[-n] + c3[-1])
    cc.4 <- yi[-n]
    to.sort <- c(test0, x3)
    ix.final <- order(to.sort)
    sorted.final <- to.sort[ix.final]
    tmp <-
        unique(data.frame(sorted.final,
                          c(ffppual(xi, cc.1,cc.2,cc.3,cc.4, test0, FALSE),
                            ffppual(xi, cc.1,cc.2,cc.3,cc.4, x3, TRUE))[ix.final]))
    ## get spline on the right timescale - kludgy
    tmp2 <- tmp
    tmp2[[1]] <- round(tmp2[[1]], 5) # tries to deal with identical() issues
    res <- tmp2[[2]][tmp2[[1]] %in% x2]
    ## deals with identical() issues via linear approx
    if(length(res) != n)
        res <- approx(x=tmp[[1]], y=tmp[[2]], xout=x2)$y
    res
}