derivative.test.R
# fields, Tools for spatial data
# Copyright 2004-2011, Institute for Mathematics Applied Geosciences
# University Corporation for Atmospheric Research
# Licensed under the GPL -- www.gpl.org/licenses/gpl.html
#library( fields, lib.loc="lib.test")
library( fields)
options(echo=FALSE)
test.for.zero.flag<- 1
DD<- cbind( seq(.01,2,,50))
look2<- RadialBasis(DD, dimension=2,M=3,derivative=1)
look1<- ( RadialBasis(DD+1e-5, dimension=2,M=3,derivative=0 )
- RadialBasis(DD-1e-5, dimension=2,M=3,derivative=0))/2e-5
test.for.zero( look1, look2,tol=1e-6, tag="radial basis function exact" )
set.seed( 234)
x<- matrix( runif(10), ncol=2)
ctest<- rep(0,5)
ctest[3]<- 1
stationary.cov( x,x, Covariance="RadialBasis", dimension=2,M=3,derivative=0)-> look0
RadialBasis( rdist(x,x), dimension=2,M=3,derivative=0)-> sanity.look
test.for.zero( look0, sanity.look, tag="sanity test of stationary.cov with RadialBasis")
Rad.cov(x,x,p= (2*3 -2))-> look1
test.for.zero( sanity.look, look1, tag="sanity test of Rad.cov")
sanity.look%*% ctest->look0
stationary.cov( x,x, Covariance="RadialBasis", dimension=2,M=3,
derivative=0, C=ctest)-> look
test.for.zero( look0, look, tag="stat.cov Radbas C multiply")
Rad.cov(x,x,p= (2*3 -2), C=ctest)-> look1
test.for.zero( look0, look1, tag="Rad.cov C multiply")
############################ end of radial basis
DD<- cbind( seq(.01,2,,50))
look2<- Wendland(DD, theta=1.0, dimension=2,k=3,derivative=1)
look1<- (Wendland(DD+1e-5, theta=1.0, dimension=2,k=3)
- Wendland(DD-1e-5, theta=1.0, dimension=2,k=3))/2e-5
test.for.zero( look1, look2,tol=1e-6)
look2<- Wendland(DD, theta=1.5, dimension=2,k=3,derivative=1)
look1<- (Wendland(DD+1e-5, theta=1.5, dimension=2,k=3)
- Wendland(DD-1e-5, theta=1.5, dimension=2,k=3))/2e-5
test.for.zero( look1, look2,tol=1e-6, tag="Wendland exact")
x<- seq( -1,1,,5)
ctest<- rep(0,5)
ctest[3]<- 1
wendland.cov( x,x, k=2, theta=.75)-> look0
Wendland( rdist(x,x)/.75, k=2, dimension=1)-> sanity.look
test.for.zero( look0, sanity.look)
look0%*% ctest->look0
wendland.cov( x,x, k=2, theta=.75, C=ctest, derivative=0)-> look
test.for.zero( look0, look, tag="Wendland C multiply")
wendland.cov( x,x, k=2, theta=1.0, C=ctest, derivative=1)-> look
wendland.cov( x+1e-5, x, k=2, theta=1.0, C=ctest)-
wendland.cov( x-1e-5, x, k=2, theta=1.0, C=ctest)-> look2
look2<- look2/2e-5
test.for.zero( look, look2,tol=1e-7, tag="Wendland.cov theta=1.0")
wendland.cov( x,x, k=2, theta=.75, C=ctest, derivative=1)-> look
wendland.cov( x+1e-5, x, k=2, theta=.75, C=ctest)-
wendland.cov( x-1e-5, x, k=2, theta=.75, C=ctest)-> look2
look2<- look2/2e-5
test.for.zero( look, look2,tol=1e-7, tag="Wendland.cov theta=.75")
stationary.cov( x,x, Covariance="Wendland", dimension=1,
k=2, theta=1.0, C=ctest, derivative=0)-> look
look0<- Wendland( rdist(x,x), k=2, dimension=1)%*%ctest
test.for.zero( look0, look, tag="stationary.cov and exact C multiply for Wendland")
wendland.cov( x,x, k=2,C=ctest, theta=1.0)-> look
look0<- Wendland( rdist(x,x), k=2, dimension=1)%*%ctest
test.for.zero( look0, look, tag=" Wendland C multiply")
####### 2 -d quadratic surface
x<- make.surface.grid( list(x=seq( -1,1,,20), y=seq( -1,1,,20)))
y<- (.123*x[,1] + .234*x[,2])
obj<- mKrig( x,y, lambda=0, cov.function="wendland.cov", k=3, theta=.4)
xp<- make.surface.grid( list(x=seq(-.5,.5,,24),y= seq( -.5,.5,,24)) )
predict( obj, xp, derivative=1)-> outd
test.for.zero( outd[,1],.123, tag="2-d derivs from wend.cov/mKrig")
test.for.zero( outd[,2],.234)
#%%%%%%%% repeat to check derivatives in stationary.cov
x<- make.surface.grid( list(x=seq( -1,1,,20), y=seq( -1,1,,20)))
y<- (.123*x[,1] + .234*x[,2])
obj<- mKrig( x,y, lambda=0, cov.function="stationary.cov",
cov.args=list(k=3, theta=.2, dimension=2, Covariance="Wendland"))
xp<- make.surface.grid( list(x=seq(-.5,.5,,24),y= seq( -.5,.5,,24)) )
predict( obj, xp, derivative=1)-> outd
test.for.zero( outd[,1],.123, tag="2-d derivs from stationary-wend/mKrig")
test.for.zero( outd[,2],.234)
############## quadratic surface
x<- make.surface.grid( list(x=seq( -1,1,,20), y=seq( -1,1,,20)))
y<- (x[,1]**2 - 2* x[,1]*x[,2] + x[,2]**2)/2
############## wendland.cov
obj<- mKrig( x,y, lambda=0, cov.function="wendland.cov", k=3, theta=.8)
xp<- make.surface.grid( list(x=seq(-.5,.5,,24),y= seq( -.5,.5,,24)) )
true<- cbind( xp[,1] - xp[,2], xp[,2]- xp[,1])
############## wendland.cov
predict( obj, xp, derivative=1)-> outd
rmse<-sqrt(mean((true[,1] - outd[,1])**2))/sqrt(mean(true[,1]**2))
test.for.zero( rmse,0, tol=5e-3,relative=FALSE, tag="wendland.cov quad 2-d")
############## stationary cov
x<- make.surface.grid( list(x=seq( -1,1,,20), y=seq( -1,1,,20)))
y<- (x[,1]**3 + x[,2]**3)
obj<- mKrig( x,y, lambda=0, cov.function="stationary.cov",
cov.args=list(k=3, theta=.8, dimension=2, Covariance="Wendland"))
xp<- make.surface.grid( list(x=seq(-.5,.5,,24),y= seq( -.5,.5,,24)) )
true<- cbind( 3*xp[,1]**2 , 3*xp[,2]**2)
predict( obj, xp, derivative=1)-> outd2
rmse<-sqrt(mean((true[,1] - outd2[,1])**2))/sqrt(mean(true[,1]**2))
test.for.zero( rmse,0, tol=1e-2,relative=FALSE,
tag="stationary.cov/Wendland cubic 2-d")
############## stationary cov with radial basis
x<- make.surface.grid( list(x=seq( -1,1,,20), y=seq( -1,1,,20)))
y<- (x[,1]**3 + x[,2]**3)
obj<- Krig( x,y, cov.function="stationary.cov", m=3,
cov.args=list(M=3, dimension=2, Covariance="RadialBasis"))
xp<- make.surface.grid( list(x=seq(-.5,.5,,24),y= seq( -.5,.5,,24)) )
true<- cbind( 3*xp[,1]**2 , 3*xp[,2]**2)
predictDerivative.Krig( obj, xp)-> outd2
look<- as.surface( xp, outd2[,1])
rmse<-sqrt(mean((true[,1] - outd2[,1])**2))/sqrt(mean(true[,1]**2))
test.for.zero( rmse,0, tol=5e-3,relative=FALSE,
tag="stationary.cov/Wendland cubic 2-d")
#########################
x<- make.surface.grid( list(x=seq( -1,1,,20), y=seq( -1,1,,20)))
y<- (x[,1]**3 + x[,2]**3)
obj<- mKrig( x,y, lambda=0, cov.function="wendland.cov", k=3, V=diag(c( 1.1,1.1) ))
xp<- make.surface.grid( list(x=seq(-.5,.5,,24),y= seq( -.5,.5,,24)) )
predict( obj, xp, derivative=1)-> outd
true<- cbind( 3*xp[,1]**2 , 3*xp[,2]**2)
rmse<-sqrt(mean((true[,1] - outd[,1])**2)/mean(true[,1]**2))
test.for.zero( rmse,0, tol=5e-3,relative=FALSE)
obj<- Tps( x,y,lambda=0)
predictDerivative.Krig( obj, xp, derivative=1)-> outd
look<- as.surface( xp, outd[,1])
rmse<-sqrt(mean((true[,1] - outd[,1])**2)/mean(true[,1]**2))
test.for.zero( rmse,0, tol=2e-4,relative=FALSE, tag="Tps derivative x1")
rmse<-sqrt(mean((true[,2] - outd[,2])**2)/mean(true[,2]**2))
test.for.zero( rmse,0, tol=2e-4,relative=FALSE, tag="Tps derivative x2")
cat("done with dervative tests", fill=TRUE)
options( echo=TRUE)