# fields is a package for analysis of spatial data written for
  # the R software environment .
  # Copyright (C) 2018
  # University Corporation for Atmospheric Research (UCAR)
  # Contact: Douglas Nychka, nychka@ucar.edu,
  # National Center for Atmospheric Research,
  # PO Box 3000, Boulder, CO 80307-3000
  #
  # This program is free software; you can redistribute it and/or modify
  # it under the terms of the GNU General Public License as published by
  # the Free Software Foundation; either version 2 of the License, or
  # (at your option) any later version.
  # This program is distributed in the hope that it will be useful,
  # but WITHOUT ANY WARRANTY; without even the implied warranty of
  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  # GNU General Public License for more details.

#library( fields, lib.loc="lib.test")

suppressMessages(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)