# R/ipol.R In chebpol: Multivariate Interpolation

#### Documented in ipol

```#' Create interpolating function.
#'
#' Create an interpolating function from given values. Several interpolation methods are
#' supported.
#'
#' \code{ipol} is a wrapper around various interpolation methods in package \pkg{chebpol}.
#' Which arguments to specify depends on the method. The interpolation methods are described
#' in \code{vignette("chebpol",package="chebpol")}.
#'
#' The method \code{"chebyshev"} needs only the number of Chebyshev knots in
#' each dimension. This is either the \code{"dim"} attribute of the array
#' \code{val}, or the \code{dims} argument if \code{val} is a function.  Also
#' the intervals can be specified if different from [-1, 1].
#'
#' The method \code{"uniform"} is similar to the \code{"chebyshev"}, but
#' uniformly spaced knots are created.  The argument
#' \code{intervals} generally goes with \code{dims} when something else than
#' standard intervals \code{[-1, 1]} are used.
#'
#' The methods \code{"multilinear"}, \code{"fh"} (Floater-Hormann), \code{"stalker"},
#' \code{"hstalker"}, and
#' \code{"general"} needs the argument \code{grid}.  These are the methods
#' which can use arbitrary Cartesian grids.
#' The stalker spline
#' is described in \code{vignette("stalker",package="chebpol")}. The Floater-Hormann
#' method (\code{"fh"}) also needs the \code{k} argument, the degree of the blending
#' polynomials. It defaults to 4.
#'
#' The method \code{"polyharmonic"} needs the arguments \code{knots} and
#' \code{k}. In addition it can take the logical argument \code{normalize}
#' for normalizing the knots to the unit hypercube. The default is \code{NA}, which
#' uses normalization if any of the knots are outside the unit hypercube. Also, a logical
#' \code{nowarn} is accepted, it is used to suppress a warning in case the system can't
#' be solved exactly and a least squares fallback method is used.
#'
#' The method \code{"simplexlinear"} needs the argument \code{knots}. It creates a
#' Delaunay triangulation from the knots, and does linear interpolation in each simplex
#' by weighting the vertex values with the barycentric coordinates.
#'
#' If knots are required, but the grid argument is given, knots are constructed as
#' \code{t(expand.grid(grid))}
#'
#' The \code{"crbf"} is the multilayer compact radial basis function
#' interpolation in ALGLIB (\url{http://www.alglib.net/interpolation/fastrbf.php}).
#' It is only available if ALGLIB was available at
#' compile time. It takes the extra arguments \code{"rbase"}, \code{"layers"}, and
#' \code{"lambda"}. These are discussed in the ALGLIB documentation.
#'
#' There are also some usage examples and more in \code{vignette("chebpol")} and \code{vignette('chebusage')}.
#'
#' @param val array or function. Function values on a grid, or the function
#' itself. If it is the values, the \code{"dim"}-attribute must be
#' appropriately set. If it is a function, it will be evaluated in the grid points.
#' @param dims integer vector. The number of grid points in each dimension. Not
#' needed if \code{val} is an array or \code{grid} is used.
#' @param intervals list of length 2 numeric vectors. The lower and upper bound
#' in each dimension. Not used if \code{grid} is specified.
#' @param grid list. Each element is a vector of ordered grid-points for a
#' dimension.  These need not be Chebyshev-knots, nor evenly spaced.
#' @param knots matrix. Each column is a point in an M-dimensional space.
#' @param k numeric. Additional value, used with some methods.
#' @param method character. The interpolation method to use.
#' @param ... Further arguments to the function, if \code{is.function(val)}. And some
#' extra arguments for interpolant creation described in section Details.
#' for the given function. The argument \code{x} can be a matrix of column
#' vectors which are evaluated in parallel in a number of threads.  The
#' function yields values for arguments outside the hypercube as well, though
#' it will typically be a poor approximation.  \code{threads} is an integer
#' specifying the number of parallel threads which should be used when
#' evaluating a matrix of column vectors.
#' @examples
#'
#' ## evenly spaced grid-points
#' su <- seq(0,1,length.out=10)
#' ## irregularly spaced grid-points
#' s <- su^3
#' ## create approximation on the irregularly spaced grid
#' ml1 <- ipol(exp, grid=list(s), method='multilin')
#' fh1 <- ipol(exp, grid=list(s), method='fh')
#' ## test it, since exp is convex, the linear approximation lies above
#' ## the exp between the grid points
#' ml1(su) - exp(su)
#' fh1(su) - exp(su)
#'
#' ## multi dimensional approximation
#' f <- function(x) 10/(1+25*mean(x^2))
#' # a 3-dimensional 10x10x10 grid, first and third coordinate are non-uniform
#' grid <- list(s, su, sort(1-s))
#'
#' # make multilinear, Floater-Hormann, Chebyshev and polyharmonic spline.
#' ml2 <- ipol(f, grid=grid, method='multilin')
#' fh2 <- ipol(f, grid=grid, method='fh')
#' hst <- ipol(f, grid=grid, method='hstalker')
#' ch2 <- ipol(f, dims=c(10,10,10), intervals=list(0:1,0:1,0:1), method='cheb')
#' knots <- matrix(runif(3*1000),3)
#' ph2 <- ipol(f, knots=knots, k=2, method='poly')
#' sl2 <- ipol(f, knots=knots, method='simplexlinear')
#' # my alglib is a bit slow, so stick to 100 knots
#' if(havealglib()) crb <- ipol(f, knots=knots[,1:100], method='crbf',
#'   rbase=2, layers=5, lambda=0)
#' # make 7 points in R3 to test them on
#' m <- matrix(runif(3*7),3)
#' rbind(true=apply(m,2,f), ml=ml2(m), fh=fh2(m), cheb=ch2(m), poly=ph2(m), sl=sl2(m),hst=hst(m),
#' crbf=if(havealglib()) crb(m) else NULL )
#'
#' @export
ipol <- function(val,dims=NULL,intervals=NULL,grid=NULL,knots=NULL,k=NULL,
method=c('chebyshev','multilinear','fh','uniform','general','polyharmonic',
'simplexlinear', 'hstalker', 'stalker','crbf', 'oldstalker'),
...) {
method <- match.arg(method)
args <- list(...)
switch(method,
chebyshev={
if(is.function(val) && is.null(dims)) stop('dims must be specified')
if(is.function(val)) return(chebappxf.real(val,dims,intervals,...))
if(is.null(dim(val)) && !is.null(dims)) dim(val) <- dims
return(chebappx.real(val,intervals))
},
multilinear={
if(is.null(grid)) stop('grid must be specified for multi linear interpolation')
if(!is.list(grid)) grid <- list(grid)
grid <- lapply(grid,as.numeric)
if(unsortedgrid(grid)) stop("grid must be distinct ordered values")
blend <- args[['blend']]
if(is.null(blend)) return(mlappx.real(val,grid,...))
newarg <- args[-match(c('blend'),names(args), nomatch=0)]
return(blenddef(do.call(mlappx.real,c(list(val,grid), newarg)),blend))
},
simplexlinear={
if(is.null(knots)) {
if(!is.null(grid))
knots <- t(expand.grid(grid))
else
stop('knots must be specified for simplex linear interpolation')
}
blend <- args[['blend']]
if(is.null(blend)) return(slappx.real(val,knots,...))
newarg <- args[-match(c('blend'),names(args), nomatch=0)]
return(blenddef(do.call(slappx.real,c(list(val,knots), newarg)),blend))
},
fh={
if(is.null(grid)) stop('grid must be specified for Floater-Hormann interpolation')
if(!is.list(grid)) grid <- list(grid)
if(unsortedgrid(grid)) stop("grid must be distinct ordered values")
grid <- lapply(grid,as.numeric)
if(is.null(k)) k <- pmin(4,sapply(grid,length)-1)
return(fhappx.real(val,grid,d=k))
},
uniform={
if(is.function(val) && is.null(dims)) stop('Must specify dims for uniform intervals')
if(is.null(intervals) && !is.null(grid)) {
intervals <- lapply(grid,range)
warning("intervals constructed from ranges of grid-argument")
}
if(is.function(val)) return(ucappxf.real(val,dims,intervals,...))
if(is.null(dim(val)) && !is.null(dims)) dim(val) <- dims
return(ucappx.real(val,intervals))
},
general={
if(is.null(grid)) stop('grid must be specified for general interpolation')
if(!is.list(grid)) grid <- list(grid)
if(unsortedgrid(grid)) stop("grid must be distinct ordered values")
grid <- lapply(grid,as.numeric)
if(is.function(val)) return(chebappxgf.real(val,grid,...))
return(chebappxg.real(val,grid))
},
polyharmonic={
if(is.null(knots)) {
if(!is.null(grid))
knots <- t(expand.grid(grid))
else
stop('Must specify knots for polyharmonic splines.')
}
if(is.null(k)) k <- 3
# Now, ... may contain normalize and nowarn in addition to
# arguments to val if it is a function
# peel off normalize and nowarn
normalize <- args[['normalize']]
if(is.null(normalize)) normalize <- NA
nowarn <- args[['nowarn']]
if(is.null(nowarn)) nowarn <- FALSE
newarg <- args[-match(c('normalize','nowarn'),names(args), nomatch=0)]
return(do.call(polyh.real,c(list(val,knots,k,normalize,nowarn),newarg)))
},
oldstalker={
if(is.null(grid)) stop('grid must be specified for stalker interpolation')
if(!is.list(grid)) grid <- list(grid)
if(unsortedgrid(grid)) stop('grid must be distinct ordered values')
grid <- lapply(grid,as.numeric)
if(is.null(k)) k = 2
blend <- args[['blend']]
if(is.null(blend)) return(oldstalkerappx(val,grid,r=k,...))
newarg <- args[-match(c('blend'),names(args), nomatch=0)]
return(blenddef(do.call(oldstalkerappx,c(list(val,grid,r=k), newarg)),blend))
},
hstalker={
if(is.null(grid)) stop('grid must be specified for stalker interpolation')
if(!is.list(grid)) grid <- list(grid)
if(unsortedgrid(grid)) stop('grid must be distinct ordered values')
grid <- lapply(grid,as.numeric)
blend <- args[['blend']]
if(is.null(blend)) return(hstalkerappx(val,grid,...))
newarg <- args[-match(c('blend'),names(args), nomatch=0)]
return(blenddef(do.call(hstalkerappx,c(list(val,grid), newarg)),blend))

},
stalker={
if(is.null(grid)) stop('grid must be specified for stalker interpolation')
if(!is.list(grid)) grid <- list(grid)
if(unsortedgrid(grid)) stop('grid must be distinct ordered values')
grid <- lapply(grid,as.numeric)
blend <- args[['blend']]
if(is.null(blend)) return(stalkerappx(val,grid,...))
newarg <- args[-match(c('blend'),names(args), nomatch=0)]
return(blenddef(do.call(stalkerappx,c(list(val,grid), newarg)),blend))

},
crbf={
if(is.null(knots)) stop('Must specify knots for radial basis functions.')
if(length(args) == 0) {
# old interface
if(is.null(k)) {
rbase <- 2; layers <- 5; lambda = 0
} else {
rbase <- k[1]; layers <- k[2]; lambda <- k[3]
}
return(rbf.alglib(val,knots,rbase,layers,lambda))
}
rbase <- args[['rbase']]; if(is.null(rbase)) rbase <- 2
layers <- args[['layers']]; if(is.null(layers)) layers <- 5L
lambda <- args[['lambda']]; if(is.null(lambda)) lambda <- 0
newarg <- args[-match(c('rbase','layers','lambda'),names(args), nomatch=0)]
return(do.call(rbf.alglib,c(list(val,knots,rbase,layers,lambda),newarg)))
},
stop('Unknown interpolation method: ', method)
)
}

unsortedgrid <- function(g) {
any(sapply(g,function(s) is.unsorted(s,strictly=TRUE) && is.unsorted(-s,strictly=TRUE)))
}
```

## Try the chebpol package in your browser

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

chebpol documentation built on Dec. 9, 2019, 5:08 p.m.