[S] RE: interpolation in m dimensions

Mark Bravington FSMG CEFAS (M.V.BRAVINGTON@cefas.co.uk)
Tue, 6 Oct 1998 16:00:43 +0100

> I have a data set allocated on a rectangular grid, say z(x,y) where
> x_seq(from=xb,to=xe,length=xl);y_seq(from=yb,to=ye,length=yl)
> Is there a command in S+ to perform interpolation of z into an arbitrary point
> X0,Y0? (Simon Rosenfeld)


the routine 'mapprox' below does linear interpolation at a set of m-dimensional points onto an m-dimensional grid of values.
If the m-dimensional grid has dim( n1, n2, ..., nm), then the points for interpolation are by default assumed to fall in the range
(1:n1), (1:n2), ..., (1:nm) (per dimension), but you can set the scale of each dimension of the grid using the x.range parameter.

It's probably wise to experiment with this routine-- I haven't used it directly for quite a while. No guarantees!

Hope this helps
Mark Bravington


> # Test interpolation on a fine grid
> yy_ outer( seq( 0, 1, 0.013), seq( 0, 1, 0.013), function( x, y) sin( x*y*pi))
> mapprox( rbind( c( 0.2, 0.7), c( 0.9, 0.4), c( 0.2, 0.3)), yy, x.range=rbind( c( 0, 0), c( 1,1)))
[,1] [,2]
[1,] 0.2 0.7
[2,] 0.9 0.4
[3,] 0.2 0.3

[1] 0.4162331 0.8928420 0.1829592
> sin( pi* c( 0.2*0.7, 0.9*0.4, 0.2*0.3))
[1] 0.4257793 0.9048271 0.1873813
> # not bad

mapprox_ function( x.out, y, x.range=rbind( 1, dim(y))) {
# Linear interpolation of grid values onto points in m dimesions.
# Points outside "x.range" are pushed back to side/edge of grid
# "y" is an m-dimensional array of size n1 x n2 x ... x n.m
# "x.range" (if supplied) is a matrix of dimension "c(2, m)"
# "x.out" is a matrix of dimension "c( n.out, m)"

n.m_ length( dim( y))
n.out_ nrow( x.out)
r_ rep( n.out, n.m)

x_ (x.out-rep( x.range[1,], r)) *
rep( dim( y)-1, r) /
rep( x.range[2,]-x.range[1,], r)
f.out_ x - floor( x)
x_ 1 + floor( x)

below_ x < 1; x[ below]_ 1; f.out[ below]_ 0
above_ x > rep( dim( y), r); x[ above]_ rep( dim( y)-1, r)[above];
f.out[ above]_ 1

dx_ as.matrix( do.call( "expand.grid", rep( list( c( 0, 1)), n.m)))
two.to.n.m_ prod( rep( 2, n.m))
nx_ array( x, c( dim( x), two.to.n.m))
nx_ nx + aperm( array( dx, c( dim( dx), n.out)), c( 3,2,1))
nx_ matrix( aperm( nx, c( 3,1,2)), n.out*two.to.n.m, n.m)
y.out_ matrix( y[nx], two.to.n.m, n.out)

fprod_ matrix( 0, n.m*n.out*two.to.n.m, 3)
fprod[,1]_ rep( 1:n.out, rep( two.to.n.m, n.out))
fprod[,2]_rep( 1:n.m, rep( n.out*two.to.n.m, n.m))
fprod[,3]_1+dx[ rep(1:two.to.n.m,n.out),]

lff_ array( log( c( 1-f.out, f.out)), c( dim( f.out), 2))
fy_ lff[ fprod]
dim( fy)_ c( n.out * two.to.n.m, n.m)
fy_ fy %*% rep.int(1,n.m)
fy_ exp( fy)
dim( fy)_ c( two.to.n.m, n.out)

y.out_ y.out*fy
y.out_ rep.int( 1, two.to.n.m) %*% y.out

return( x=x.out, y=as.vector( y.out))

This message was distributed by s-news@wubios.wustl.edu. To unsubscribe
send e-mail to s-news-request@wubios.wustl.edu with the BODY of the
message: unsubscribe s-news