R/interp_Grid.R

In IRSN/smint: Smooth Multivariate Interpolation for Gridded and Scattered Data

#' Grid linear interpolation in arbitrary dimension through repeated
#' one-dimensional interpolations.
#'
#' The one-dimensional interpolation is performed using an
#' interpolation function given in \code{interpFun1d} or a Cardinal
#' Basis function given in \code{cardinalBasis1d}.  In the second
#' case, it is required here that the one-dimensional interpolation
#' method is \emph{linear} w.r.t. the vector of interpolated values.
#' For each dimension, a one-dimensional interpolation is carried
#' out, leading to a collection of interpolation problems each with a
#' dimension reduced by one.
#'
#' When a cardinal basis is used, the \emph{same basis} is used for
#' all interpolations w.r.t. the same variable, so the number of call
#' to the \code{cardinalBasis1d} function is equal to the interpolation
#' dimension, making this method \emph{very fast} compared to the
#' general grid interpolation, see the \bold{Examples} section.
#'
#' By default, i.e. when both \code{interpFun1d} and \code{cardinalBasis1d}
#' are \code{NULL}, a Lagrange cardinal basis interpolation is used. So
#' the 1d interpolation is the broken line interpolation.
#' 
#' @title Grid interpolation in arbitrary dimension through Cardinal
#' Basis
#'
#' @aliases interp_Grid
#' 
#' @usage
#' interp_Grid(X, Y, Xout,
#'             interpFun1d = NULL,
#'             cardinalBasis1d = NULL,
#'             intOrder = NULL,
#'             useC = TRUE, trace = 1L,
#'             ...)
#' 
#' @param X An object that can be coerced into \code{Grid}. This
#' can be a data.frame or a matrix in \emph{Scattered Data style} in
#' which case the column number is equal to the spatial dimension
#' \eqn{d}, and the row number is then equal to the number of nodes
#' \eqn{n}. But it can also be a \code{Grid} object previously
#' created.  A data frame or matrix \code{X} will first be coerced
#' into \code{Grid} by using the the S3 method
#' \code{\link{as.Grid}}.
#' 
#' @param Y Response to be interpolated. It must be a vector of
#' length \eqn{n} equal to the number of nodes. When \code{X} has
#' class \code{Grid}, the order of the elements in \code{Y} must
#' conform to the order of the nodes as given in \code{X}.
#'
#' @param Xout Interpolation locations. Can be a vector or a
#' matrix. In the first case, the length of the vector must be equal
#' to the spatial dimension \eqn{d} as given by \code{xLev} or
#' \code{X}.  In the second case, each row will be considered as a
#' response to be interpolated, and the number of columns of
#' \code{Xout} must be equal to the spatial dimension.
#'
#' @param interpFun1d The function to interpolate in 1d. This
#' function must have as its first 3 formals 'x', 'y', and 'xout', as
#' does \code{\link{approx}}.  It must also return the vector of
#' interpolated values as DOES NOT \code{approx}. In most cases, a
#' simple wrapper will be enough to use an existing method of
#' interpolation, see \bold{Examples}.
#' 
#' @param cardinalBasis1d Function evaluating an interpolation
#' Cardinal Basis. This function must have as its first 2 formals
#' 'x', and 'xout'.  It must return a matrix with \code{length(x)}
#' columns and \code{length(xout)} rows. The \eqn{j}-th column is the
#' vector of the values of the \eqn{j}-th cardinal basis function on
#' the vector \code{xout}. \emph{This formal argument can be used
#' only when} \code{interpFun1d} \emph{is not used}.
#'
#' @param trace Level of verbosity.
#'
#' @param intOrder Order of the one-dimensional interpolations.  Used
#' only when \code{interpFun1d} is given, and not when
#' \code{cardinalBasis1d} is given. Must be a permutation of
#' \code{1:d} where \code{d} is the spatial dimension. This argument
#' is similar to the argument of the \code{\link{aperm}} method.  By
#' default, the interpolation order is \eqn{d}, \eqn{d-1}, \dots,
#' \eqn{1}, corresponding to \code{intOrder = d:1L}. Note that for
#' the first element of \code{intOrder}, a vector of length
#' \code{nrow(Xout)} is passed as \code{xout} formal to
#' \code{interpFun1d}, while the subsequent calls to
#' \code{interpFun1} use a \code{xout} of length \code{1}. So the
#' choice of the first element can have an impact on the computation
#' time.
#' 
#' @param useC Logical. Used only when \code{interpFun1d} is given,
#' and not when \code{cardinalBasis1d} is given. If \code{TRUE}, the
#' computation is done using ' a C program via the
#' \code{.Call}. Otherwise, the computation is ' done entirely in R
#' by repeated use of the \code{apply} function.  ' Normally
#' \code{useC = TRUE} should be faster, but it is not always ' the
#' case.
#'
#' @param ... Further arguments to be passed to \code{interpCB}. NOT
#' IMPLEMENTED YET.
#' 
#' @export
#'
#' @importFrom stats approx
#' 
#' @return A single interpolated value if \code{Xout} is a vector or
#' a row matrix.  If \code{Xout} is a matrix with several rows, the
#' result is a vector of interpolated values, in the order of the
#' rows of \code{Xout}.
#' 
#' @author Yves Deville
#'
#' @examples
#' set.seed(12345)
#' ##========================================================================
#' ## Select Interpolation Function. This function must have its first 3
#' ## formals 'x', 'y', and 'xout', as does 'approx'. It must also return
#' ## the vector of interpolated values as DOES NOT 'approx'. So a wrapper
#' ## must be xritten.
#' ##=======================================================================
#' myInterpFun <- function(x, y, xout) approx(x = x, y = y, xout = xout)$y
#'
#' ##=======================================================================
#' ## Example 1
#' ## ONE interpolation, d = 2. 'Xout' is a vector.
#' ##=======================================================================
#' myFun1 <- function(x) exp(-x[1]^2 - 3 * x[2]^2)
#' myGD1 <- Grid(nlevels = c("X" = 8, "Y" = 12))
#' Y1 <- apply_Grid(myGD1, myFun1)
#' Xout1 <- runif(2)
#' GI1 <- interp_Grid(X = myGD1,  Y = Y1, Xout = Xout1,
#'                    interpFun1d = myInterpFun)
#' c(true = myFun1(Xout1), interp = GI1)
#'
#' ##=======================================================================
#' ## Example 2
#' ## ONE interpolation, d = 7. 'Xout' is a vector.
#' ##=======================================================================
#' d <- 7; a <- runif(d); myFun2 <- function(x) exp(-crossprod(a, x^2))
#' myGD2 <- Grid(nlevels = rep(4L, time = d))
#' Y2 <- apply_Grid(myGD2, myFun2)
#' Xout2 <- runif(d)
#' GI2 <- interp_Grid(X = myGD2,  Y = Y2, Xout = Xout2,
#'                    interpFun1d = myInterpFun)
#' c(true = myFun2(Xout2), interp = GI2)
#'
#' ##=======================================================================
#' ## Example 3
#' ## 'n' interpolations, d = 7. 'Xout' is a matrix. Same grid data and
#' ## response as before
#' ##=======================================================================
#' n <- 30
#' Xout3 <- matrix(runif(n * d), ncol = d)
#' GI3 <- interp_Grid(X = myGD2,  Y = Y2, Xout = Xout3,
#'                    interpFun1d = myInterpFun)
#' cbind(true = apply(Xout3, 1, myFun2), interp = GI3)
#' 
#' ##======================================================================
#' ## Example 4
#' ## 'n' interpolation, d = 5. 'Xout' is a matrix. Test the effect of the
#' ## order of interpolation.
#' ##=======================================================================
#' d <- 5; a <- runif(d); myFun4 <- function(x) exp(-crossprod(a, x^2))
#' myGD4 <- Grid(nlevels = c(3, 4, 5, 2, 6))
#' Y4 <- apply_Grid(myGD4, myFun4)
#' n <- 100
#' Xout4 <- matrix(runif(n * d), ncol = d)
#' t4a <- system.time(GI4a <- interp_Grid(X = myGD4,  Y = Y4, Xout = Xout4,
#'                                        interpFun1d = myInterpFun))
#' t4b <- system.time(GI4b <- interp_Grid(X = myGD4,  Y = Y4, Xout = Xout4,
#'                                        interpFun1d = myInterpFun,
#'                                        intOrder = 1L:5L))
#' cbind(true = apply(Xout4, 1, myFun4), inta = GI4a, intb = GI4b)
#'
#' ##======================================================================
#' ## Example 5 
#' ## 'n' interpolation, d = 5 using a GIVEN CARDINAL BASIS.
#' ## 'Xout' is a matrix. 
#' ##=======================================================================
#'
#' ## Natural spline for use through Cardinal Basis in 'gridIntCB'
#' myInterpCB5 <-  function(x, xout) cardinalBasis_natSpline(x = x, xout = xout)$CB
#' 
#' ## Natural spline for use through an interpolation function
#' myInterpFun5 <- function(x, y, xout) {
#'      spline(x = x, y = y, n = 3 * length(x), method = "natural", xout = xout)$y
#' }
#' 
#' n <- 10
#' ## Since we use a natural spline, there must be at least 3 levels for
#' ## each dimension.
#' myGD5 <- Grid(nlevels = c(3, 4, 5, 4, 6))
#' Y5 <- apply_Grid(myGD5, myFun4)
#' Xout5 <- matrix(runif(n * d), ncol = d)
#' 
#' t5a <- system.time(GI5a <- interp_Grid(X = myGD5, Y = Y5, Xout = Xout5,
#'                                        interpFun1d = myInterpFun5))
#' t5b <- system.time(GI5b <- interp_Grid(X = myGD5, Y = Y5, Xout = Xout5,
#'                                        interpFun1d = myInterpFun5,
#'                                        useC = FALSE))
#' t5c <- system.time(GI5c <- interp_Grid(X = myGD5, Y = Y5, Xout = Xout5,
#'                                        cardinalBasis1d = myInterpCB5))
#' df <- data.frame(true = apply(Xout5, 1, myFun4),
#'                  gridInt_C = GI5a, gridInt_R = GI5b, gridIntCB = GI5c)
#' head(df)
#' rbind(gridInt_C = t5a, gridInt_R = t5b, gridIntCB = t5c)
#' 
interp_Grid <- function(X, Y, Xout,
                        interpFun1d = NULL,
                        cardinalBasis1d = NULL,
                        intOrder = NULL,
                        useC = TRUE,
                        trace = 1L,
                        ...) {

    if (is.null(cardinalBasis1d)) {
        if (is.null(interpFun1d)) {
            type <- "CB"
            cardinalBasis1d <- function(x, xout){
                cardinalBasis_lagrange(x = x, xout = xout)$CB
            }
        } else {
            type <- "fun"
        }
    } else {
        if (!is.null(interpFun1d)) {
            stop("only one of the formals 'interpFun1d' and",
                 " 'cardinalBasis1s' can be given")
        }
        type <- "CB"
    }

    if (type == "CB") {
        res <- gridIntCB(X, Y, Xout,
                         interpCB = cardinalBasis1d,
                         intOrder = intOrder,
                         trace = trace,
                         ...) 
    } else {
        res <- gridInt(X, Y, Xout,
                       interpFun = interpFun1d,
                       intOrder = intOrder,
                       useC = useC,
                       trace = trace,
                       ...) 
    }
    
    res
    
}