R/interpolate.R

In e1071: Misc Functions of the Department of Statistics, Probability Theory Group (Formerly: E1071), TU Wien

interpolate <- function(x, a, adims=lapply(dimnames(a), as.numeric),
                        method="linear"){

    if(is.vector(x)) x<- matrix(x, ncol=length(x))

    if(!is.array(a))
        stop("a is not an array")
    
    ad <- length(dim(a))
    
    method <- pmatch(method, c("linear", "constant"))
    if (is.na(method)) 
        stop("invalid interpolation method")
    
    if(any(unlist(lapply(adims, diff))<0))
        stop("dimensions of a not ordered")
    
    retval <- rep(0, nrow(x))
    bincombi <- bincombinations(ad)

    convexcoeff <- function(x, y) {
        ok <- y>0
        x[ok] <- y[ok]-x[ok]
        x
    }

    for(n in 1:nrow(x)){
        
        ## the "leftmost" corner of the enclosing hypercube 
        leftidx <- rep(0, ad)
        xabstand <- rep(0, ad)
        aabstand <- rep(0, ad)
        
        for(k in 1:ad){
            if(x[n,k] < min(adims[[k]]) || x[n,k] > max(adims[[k]]))
                stop("No extrapolation allowed")
            else{
                leftidx[k] <- max(seq(adims[[k]])[adims[[k]] <= x[n,k]])
                ## if at the right border, go one step to the left
                if(leftidx[k] == length(adims[[k]]))
                    leftidx[k] <- leftidx[k] - 1
                
                xabstand[k] <- x[n,k] - adims[[k]][leftidx[k]]
                aabstand[k] <- adims[[k]][leftidx[k]+1] - 
                    adims[[k]][leftidx[k]]
            }
        }

        coefs <- list()
        if(method==1){
            for(k in 1:(2^ad)){
                retval[n] <- retval[n] +
                    element(a, leftidx+bincombi[k,]) *
                        prod((aabstand-
                              convexcoeff(xabstand,
                                          aabstand*bincombi[k,]))/aabstand)
            }
        }
        else if(method==2){
            retval[n] <- element(a, leftidx)
        }
    }
    
    names(retval) <- rownames(x)
    retval
}