# fields, Tools for spatial data
# Copyright 2004-2007, 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)
predict.derivative( 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)
  predict.derivative( 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")
### test of surface evaluation function
  outd2<-predict.surface.derivative( obj, list(x=seq(-.5,.5,,24),y= seq( -.5,.5,,24)))
  look<- as.surface( xp, outd[,1])
  test.for.zero( outd2$z[,,1],look$z, tag="Tps derivative x1 surface function")
  look<- as.surface( xp, outd[,2])
  test.for.zero( outd2$z[,,2],look$z, tag="Tps derivative x2 surface function")

cat("done with dervative tests", fill=TRUE)
options( echo=TRUE)