Commit 53406f4e authored by agebhard's avatar agebhard
Browse files

Imported sources

parents
This diff is collapsed.
interp Gridded Bivariate Interpolation for Irregular Data
interpp Pointwise Bivariate Interpolation for Irregular Data
"interp"<-function(x, y, z, xo = seq(min(x), max(x), length = 40),
yo = seq(min(y), max(y), length = 40),
ncp = 0, extrap = FALSE, duplicate = "error", dupfun = NULL)
{
if(!(all(is.finite(x)) && all(is.finite(y)) && all(is.finite(z))))
stop("missing values and Infs not allowed")
if(ncp>25){
ncp <- 25
cat("ncp too large, using ncp=25\n")
}
drx <- diff(range(x))
dry <- diff(range(y))
if(drx == 0 || dry == 0)
stop("all data collinear") # other cases caught in Fortran code
if(drx/dry > 10000 || drx/dry < 0.0001)
stop("scales of x and y are too dissimilar")
n <- length(x)
nx <- length(xo)
ny <- length(yo)
if(length(y) != n || length(z) != n)
stop("Lengths of x, y, and z do not match")
xy <- paste(x, y, sep =",")
i <- match(xy, xy)
if(duplicate=="user" && !is.function(dupfun))
stop("duplicate=\"user\" requires dupfun to be set to a function")
if(duplicate!="error")
{
centre <- function(x) {
switch(duplicate,
mean = mean(x),
median = median(x),
user = dupfun(x))
}
if(duplicate!="strip"){
z <- unlist(lapply(split(z,i), centre))
ord <- !duplicated(xy)
x <- x[ord]
y <- y[ord]
n <- length(x)
}
else{
ord <- (hist(i,plot=F,freq=T,breaks=seq(0.5,max(i)+0.5,1))$counts==1)
x <- x[ord]
y <- y[ord]
z <- z[ord]
n <- length(x)
}
}
else
if(any(duplicated(xy)))
stop("duplicate data points")
zo <- matrix(0, nx, ny)
storage.mode(zo) <- "double"
miss <- !extrap #if not extrapolating use missing values
misso <- matrix(miss, nx, ny)
if(extrap & ncp == 0)
warning("Cannot extrapolate with linear option")
ans <- .Fortran("idsfft",
as.integer(1),
as.integer(ncp),
as.integer(n),
as.double(x),
as.double(y),
as.double(z),
as.integer(nx),
as.integer(ny),
x = as.double(xo),
y = as.double(yo),
z = zo,
integer((31 + ncp) * n + nx * ny),
double(5 * n),
misso = as.logical(misso))
temp <- ans[c("x", "y", "z", "misso")]
temp$z[temp$misso]<-NA
temp[c("x", "y", "z")]
}
"interpp"<-function(x, y, z, xo, yo, ncp = 0, extrap = FALSE,
duplicate = "error", dupfun = NULL)
{
if(!(all(is.finite(x)) && all(is.finite(y)) && all(is.finite(z))))
stop("missing values and Infs not allowed")
if(is.null(xo))
stop("xo missing")
if(is.null(yo))
stop("yo missing")
if(ncp>25){
ncp <- 25
cat("ncp too large, using ncp=25\n")
}
drx <- diff(range(x))
dry <- diff(range(y))
if(drx == 0 || dry == 0)
stop("all data collinear") # other cases caught in Fortran code
if(drx/dry > 10000 || drx/dry < 0.0001)
stop("scales of x and y are too dissimilar")
n <- length(x)
np <- length(xo)
if(length(yo)!=np)
stop("length of xo and yo differ")
if(length(y) != n || length(z) != n)
stop("Lengths of x, y, and z do not match")
xy <- paste(x, y, sep =",")
i <- match(xy, xy)
if(duplicate=="user" && !is.function(dupfun))
stop("duplicate=\"user\" requires dupfun to be set to a function")
if(duplicate!="error")
{
centre <- function(x) {
switch(duplicate,
mean = mean(x),
median = median(x),
user = dupfun(x)
)
}
if(duplicate!="strip"){
z <- unlist(lapply(split(z,i), centre))
ord <- !duplicated(xy)
x <- x[ord]
y <- y[ord]
n <- length(x)
}
else{
ord <- (hist(i,plot=F,freq=T,breaks=seq(0.5,max(i)+0.5,1))$counts==1)
x <- x[ord]
y <- y[ord]
z <- z[ord]
n <- length(x)
}
}
else
if(any(duplicated(xy)))
stop("duplicate data points")
zo <- rep(0, np)
storage.mode(zo) <- "double"
miss <- !extrap #if not extrapolating use missing values
misso <- seq(miss, np)
if(extrap & ncp == 0)
warning("Cannot extrapolate with linear option")
ans <- .Fortran("idbvip",
as.integer(1),
as.integer(ncp),
as.integer(n),
as.double(x),
as.double(y),
as.double(z),
as.integer(np),
x = as.double(xo),
y = as.double(yo),
z = zo,
integer((31 + ncp) * n + np),
double(8 * n),
misso = as.logical(misso))
temp <- ans[c("x", "y", "z", "misso")]
temp$z[temp$misso]<-NA
temp[c("x", "y", "z")]
}
library.dynam("akima")
This library contains an R implementation of the S-Plus function interp().
Version 0.1-1
I used the S-Plus version of interp() (S-Plus 3.3 f. Win 3.11) as starting
point. Then I searched for Akimas idsfft function, which is referenced in
interp().
I found the appropriate Fortran code (1978) in the ACM Collected
Algorithms archive under
http://www.netlib.org/toms/526
(also included here) and splitted it into the Fortran files under src/ .
The test driver ttidbs is also included and can be compiled with "make test".
However, it seems that this code differs a little bit from that used in S-Plus:
It implements "only" bicubic spline interpolation, no linear interpolation as
in S-Plus.
So I modified IDSFFT and added a subroutine IDPTLI which does linear
interpolation within the triangles, generated by IDTANG.
Further changes are
+ REAL -> DOUBLE PRECISION
+ static DIMENSIONs replaced with dynamic
+ option to toggle extrapolation outside of convex hull added in IDSFFT and
IDBVIP.
Because I don't know how to generate NAs in Fortran (I use g77 on Linux),
I added a logical array MISSI, that indicates if a returned value should
be NA. These values will be set to NA after the Fortran call.
+ option to handle duplicate data points added (according to an example in the
S-Plus help page)
+ man pages converted and rewritten
+ data set akima (from S-Plus) added
+ function interpp() added, it evaluates the interpolated function at
arbitraryly choosen points and generates no regular grid as interp() does.
There where some problems with interpp() when using the Fortran version:
- it crashes when compiled with "g77 -O2 -fpic" and called with more than
one output point.
- compilation with "g77 -g -fpic" fails (see src/Makefile for details)
- compilation with "g77 -g" works (and no crashes occur)
These problems do not occur in the C Version (generated by f2c), so it would
be better to use only the src-c tree.
After I finished the above steps I found a more recent version (ACM 761, 1996)
of Akimas interpolation code which uses the tripack package (also available at
ACM as algorithm no. 751) for triangulation and now I'm trying to use it for
the next version of interp().
------------------------------------------------------------------
Albrecht Gebhardt email: albrecht.gebhardt@uni-klu.ac.at
Institut fuer Mathematik Tel. : (++43 463) 2700/837
Universitaet Klagenfurt Fax : (++43 463) 2700/834
Villacher Str. 161
A-9020 Klagenfurt, Austria
------------------------------------------------------------------
akima interpolation of irregularly spaced data
"akima" <-
list(x = c(11.16, 24.2, 12.85, 19.85, 10.35, 24.65, 19.72, 15.91,
0, 20.87, 6.71, 3.45, 19.99, 14.26, 10.28, 4.51, 17.43, 22.8,
0, 7.58, 16.7, 6.08, 1.99, 25, 14.9, 3.22, 0, 9.66, 2.56, 5.22,
11.77, 17.25, 15.1, 25, 12.13, 25, 22.33, 11.52, 14.59, 15.2,
7.54, 5.23, 17.32, 2.14, 0.51, 22.69, 25, 5.47, 21.67, 3.31), y = c(1.24,
16.23, 3.06, 10.72, 4.11, 2.4, 1.39, 7.74, 20, 20, 6.26, 12.78,
4.62, 17.87, 15.16, 20, 3.46, 12.39, 4.48, 1.98, 19.65, 4.58,
5.6, 11.87, 3.12, 16.78, 0, 20, 3.02, 14.66, 10.47, 19.57, 17.19,
3.87, 10.79, 0, 6.21, 8.53, 8.71, 0, 10.69, 10.72, 13.78, 15.03,
8.37, 19.63, 20, 17.13, 14.36, 0.13), z = c(22.15, 2.83, 22.11,
7.97, 22.33, 10.25, 16.83, 15.3, 34.6, 7.54, 30.97, 41.24, 14.72,
10.74, 21.59, 15.61, 18.6, 5.47, 61.77, 29.87, 6.31, 35.74, 51.81,
4.4, 21.7, 39.93, 58.2, 4.73, 50.55, 40.36, 13.62, 6.43, 12.57,
8.74, 13.71, 12, 10.25, 15.74, 14.81, 21.6, 19.31, 26.5, 12.11,
53.1, 49.43, 3.25, 0.6, 28.63, 5.52, 44.08))
akima Waveform Distortion Data for Bivariate Interpolation
\name{akima}
\title{
Waveform Distortion Data for Bivariate Interpolation
}
\arguments{
\item{x,y,z}{
represents a smooth surface of \code{z} values at selected points
irregularly distributed in the \code{x-y} plane.
}}
\section{SOURCE}{
Hiroshi Akima, "A Method of Bivariate Interpolation and
Smooth Surface Fitting for Irregularly Distributed Data
Points",
ACM Transactions on Mathematical Software,
Vol. 4, No. 2, June 1978, pp. 148-159.
Copyright 1978, Association for Computing Machinery, Inc.,
reprinted by permission.
The data was taken from a study of waveform distortion in
electronic circuits, described in:
Hiroshi Akima, "A Method of Bivariate Interpolation and
Smooth Surface Fitting Based on Local Procedures",
CACM,
Vol. 17, No. 1, January 1974, pp. 18-20.
}
\keyword{sysdata}
% Converted by Sd2Rd version 0.2-a3.
\name{interp}
\title{
Gridded Bivariate Interpolation for Irregular Data
}
\usage{
interp(x, y, z, xo=<<see below>>, yo=<<see below>>, ncp=0, extrap=F)
}
\arguments{
\item{x}{
vector of x-coordinates of data points.
Missing values are not accepted.
}
\item{y}{
vector of y-coordinates of data points.
Missing values are not accepted.
}
\item{z}{
vector of z-coordinates of data points.
Missing values are not accepted.
\code{x}, \code{y}, and \code{z} must be the same length and may contain no fewer
than four points. The points of \code{x} and \code{y}
cannot be collinear, i.e, they cannot fall on the same line (two vectors
\code{x} and \code{y} such that \code{y = ax + b} for some \code{a}, \code{b} will not be
accepted). \code{interp} is meant for cases in which you have \code{x}, \code{y}
values scattered over a plane and a \code{z} value for each. If, instead,
you are trying to evaluate a mathematical function, or get a graphical
interpretation of relationships that can be described by a polynomial,
try \code{outer()}.
}
\item{xo}{
vector of x-coordinates of output grid. The default is 40 points
evenly spaced over the range of \code{x}. If extrapolation is not being
used (\code{extrap=F}, the default), \code{xo} should have a range that is
close to or inside of the range of \code{x} for the results to be meaningful.
}
\item{yo}{
vector of y-coordinates of output grid. The default is 40 points
evenly spaced over the range of \code{y}. If extrapolation is not being
used (\code{extrap=F}, the default), \code{yo} should have a range that is
close to or inside of the range of \code{y} for the results to be meaningful.
}
\item{ncp}{
number of additional points to be used in computing partial
derivatives at each data point.
\code{ncp} must be either \code{0} (partial derivatives are not used), or at
least 2 but smaller than the number of data points (and smaller than 25).
}
\item{extrap}{
logical flag: should extrapolation be used outside of the
convex hull determined by the data points?
}
\item{duplicate}{
indicates how to handle duplicate data points. Possible values are
\code{"error"} - produces an error message, \code{"strip"} - remove
duplicate z values, \code{"mean"},\code{"median"},\code{"user"} -
calculate mean , median or user defined function of duplicate z
values.}
\item{dupfun}{this function is applied to duplicate points if \code{duplicate="user"}
}
}
\value{
list with 3 components:
\item{x}{
vector of x-coordinates of output grid, the same as the input
argument \code{xo}, if present. Otherwise, a vector 40 points evenly spaced
over the range of the input \code{x}.
}
\item{y}{
vector of y-coordinates of output grid, the same as the input
argument \code{yo}, if present. Otherwise, a vector 40 points evenly spaced
over the range of the input \code{x}.
}
\item{z}{
matrix of fitted z-values. The value \code{z[i,j]} is computed
at the x,y point \code{x[i], y[j]}. \code{z} has
dimensions \code{length(x)} times \code{length(y)} (\code{length(xo)} times \code{length(yo)}).
}}
\section{NOTE}{
The resulting structure is suitable for input to the
functions \code{contour} and \code{image}. Check the requirements of
these functions when choosing values for \code{xo} and \code{yo}.
}
\description{
If \code{ncp} is zero, linear
interpolation is used in the triangles bounded by data points.
Cubic interpolation is done if partial derivatives are used.
If \code{extrap} is \code{FALSE}, z-values for points outside the convex hull are
returned as \code{NA}.
No extrapolation can be performed if \code{ncp} is zero.
The \code{interp} function handles duplicate \code{(x,y)} points
in different ways. As default it will stop with an error message. But
it can give duplicate points an unique \code{z} value according to the
parameter \code{duplicate} (\code{mean},\code{median} or any other user defined function).
The triangulation scheme used by \code{interp} works well if \code{x} and \code{y} have
similar scales but will appear stretched if they have very different
scales. The spreads of \code{x} and \code{y} must be within four orders of magnitude
of each other for \code{interp} to work.
}
\references{
Akima, H. (1978). A Method of Bivariate Interpolation and
Smooth Surface Fitting for Irregularly Distributed Data Points.
ACM Transactions on Mathematical Software,
\bold{4}, 148-164.
}
\seealso{
\code{\link{contour}}, \code{\link{image}}, \code{\link{approx}}, \code{\link{spline}}, \code{\link{outer}}, \code{\link{expand.grid}}.
}
\examples{
data(akima)
# linear interpolation
akima.li <- interp(akima$x, akima$y, akima$z)
image(akima.li$x,akima.li$y,akima.li$z)
contour(akima.li$x,akima.li$y,akima.li$z,add=T)
points(akima$x,akima$y)
# increase smoothness
akima.smooth <- interp(akima$x, akima$y, akima$z,
xo=seq(0,25, length=100), yo=seq(0,20, length=100))
image(akima.smooth$x,akima.smooth$y,akima.smooth$z)
contour(akima.smooth$x,akima.smooth$y,akima.smooth$z,add=T)
points(akima$x,akima$y)
# use only 15 points (interpolation only within convex hull!)
akima.part <- interp(akima$x[1:15],akima$y[1:15],akima$z[1:15])
image(akima.part$x,akima.part$y,akima.part$z)
contour(akima.part$x,akima.part$y,akima.part$z,add=T)
points(akima$x[1:15],akima$y[1:15])
# spline interpolation, use 5 points to calculate derivatives
akima.spl <- interp(akima$x, akima$y, akima$z,
xo=seq(0,25, length=100), yo=seq(0,20, length=100),ncp=5)
image(akima.spl$x,akima.spl$y,akima.spl$z)
contour(akima.spl$x,akima.spl$y,akima.spl$z,add=T)
points(akima$x,akima$y)
# example with duplicate points
data(airquality)
air <- airquality[(!is.na(airquality$Temp) &
!is.na(airquality$Ozone) &
!is.na(airquality$Solar.R)),]
# gives an error:
air.ip <- interp(air$Temp,air$Solar.R,air$Ozone)
# mean of duplicate points:
air.ip <- interp(air$Temp,air$Solar.R,air$Ozone,duplicate="mean")
}
\keyword{dplot}
% Converted by Sd2Rd version 0.2-a3.
\name{interpp}
\title{
Pointwise Bivariate Interpolation for Irregular Data
}
\usage{
interpp(x, y, z, xo, yo, ncp=0, extrap=F)
}
\arguments{
\item{x}{
vector of x-coordinates of data points.
Missing values are not accepted.
}
\item{y}{
vector of y-coordinates of data points.
Missing values are not accepted.
}
\item{z}{
vector of z-coordinates of data points.
Missing values are not accepted.
\code{x}, \code{y}, and \code{z} must be the same length and may contain no fewer
than four points. The points of \code{x} and \code{y}
cannot be collinear, i.e, they cannot fall on the same line (two vectors
\code{x} and \code{y} such that \code{y = ax + b} for some \code{a}, \code{b} will not be
accepted).
}
\item{xo}{
vector of x-coordinates of points at which to evaluate the interpolating
function.}
\item{yo}{
vector of y-coordinates of points at which to evaluate the interpolating
function.}
\item{ncp}{
number of additional points to be used in computing partial
derivatives at each data point.
\code{ncp} must be either \code{0} (partial derivatives are not used, =
linear interpolation), or at
least 2 but smaller than the number of data points (and smaller than 25).
}
\item{extrap}{
logical flag: should extrapolation be used outside of the
convex hull determined by the data points?}
\item{duplicate}{
indicates how to handle duplicate data points. Possible values are
\code{"error"} - produces an error message, \code{"strip"} - remove
duplicate z values, \code{"mean"},\code{"median"},\code{"user"} -
calculate mean , median or user defined function of duplicate z
values.
}
\item{dupfun}{this function is applied to duplicate points if \code{duplicate="user"}
}
}
\value{
list with 3 components:
\item{x}{
vector of x-coordinates of output points, the same as the input
argument \code{xo}.
}
\item{y}{
vector of y-coordinates of output points, the same as the input
argument \code{yo}.
}
\item{z}{
fitted z-values. The value \code{z[i]} is computed
at the x,y point \code{x[i], y[i]}.
}
}
\section{NOTE}{
Use \code{interp} if interpolation on a regular grid is wanted.
}
\description{
If \code{ncp} is zero, linear
interpolation is used in the triangles bounded by data points.
Cubic interpolation is done if partial derivatives are used.
If \code{extrap} is \code{FALSE}, z-values for points outside the convex hull are
returned as \code{NA}.
No extrapolation can be performed if \code{ncp} is zero.
The \code{interpp} function handles duplicate \code{(x,y)} points in different ways. As default it will stop with an error message. But
it can give duplicate points an unique \code{z} value according to the
parameter \code{duplicate} (\code{mean},\code{median} or any other user defined function).
The triangulation scheme used by \code{interp} works well if \code{x} and \code{y} have
similar scales but will appear stretched if they have very different
scales. The spreads of \code{x} and \code{y} must be within four orders of magnitude
of each other for \code{interpp} to work.
}
\references{
Akima, H. (1978). A Method of Bivariate Interpolation and
Smooth Surface Fitting for Irregularly Distributed Data Points.
ACM Transactions on Mathematical Software,
\bold{4}, 148-164.
}
\seealso{
\code{\link{contour}}, \code{\link{image}}, \code{\link{approx}}, \code{\link{spline}}, \code{\link{outer}}, \code{\link{expand.grid}},\code{\link{interp}}.
}
\examples{
data(akima)
# linear interpolation at points (1,2), (5,6) and (10,12)
akima.lip<-interpp(akima$x, akima$y, akima$z,c(1,5,10),c(2,6,12))
}
\keyword{dplot}
LIBNAME=akima
LD=ld
OBJS= idbvip.o idgrid.o idpdrv.o idsfft.o idxchg.o idcldp.o \
idlctn.o idptip.o idtang.o idptli.o
$(LIBNAME): $(OBJS)
@$(LD) $(SHLIBLDFLAGS) -o $(LIBNAME).so $(OBJS) -lf2c
clean:
@rm -f *.o *.so ttidbs
realclean:
@rm -f Makefile *.o *.so
test: ttidbs
ttidbs: ttidbs.c $(OBJS)
$(CC) ttidbs.c $(OBJS) -o ttidbs -lm -lf2c
# compilation with f77 -g -fpic fails on Linux with g77 0.5.19.1!
# (assembler complains about unknown i386 instructions, so it seems to be
# an g77 internal problem)
# But it works when the PICFLAG is omitted:
#%.o: %.f
# $(F77) -g -c $< -o $@
# compilation with f77 -O2 -fpic works, but now interpp() (it calls idbvip)
# crashes R!
/* ../src/idbvip.f -- translated by f2c (version 19950110).
You must link the resulting object file with the libraries:
-lf2c -lm (in that order)
*/
#include "f2c.h"
/* Common Block Declarations */
struct {
integer nit;
} idlc_;
#define idlc_1 idlc_
struct {
integer itpv;
} idpi_;
#define idpi_1 idpi_
/* Table of constant values */
static integer c__1 = 1;
/* Subroutine */ int idbvip_(md, ncp, ndp, xd, yd, zd, nip, xi, yi, zi, iwk,
wk, missi)
integer *md, *ncp, *ndp;
doublereal *xd, *yd, *zd;
integer *nip;
doublereal *xi, *yi, *zi;
integer *iwk;
doublereal *wk;
logical *missi;
{
/* Initialized data */
static integer lun = 6;
/* Format strings */
static char fmt_2090[] = "(1x/\002 *** IMPROPER INPUT PARAMETER VALUE(\
S).\002/\002 MD =\002,i4,10x,\002NCP =\002,i6,10x,\002NDP =\002,i6,10x,\
\002NIP =\002,i6/\002 ERROR DETECTED IN ROUTINE IDBVIP\002/)";
/* System generated locals */
integer i__1, i__2;
/* Builtin functions */
integer s_wsfe(), do_fio(), e_wsfe();
/* Local variables */