# install.packages("tgp")
library(tgp)
help(package="tgp")
# vignette("tgp")

library(MASS)
plot(mcycle)

moto.bgp = bgp(mcycle[,1],mcycle[,2],XX=seq(from=2,to=58,by=1),
  BTE = c(500, 2000, 1),pred.n=F,verb=2)
plot(moto.bgp)
par(mfrow=c(1,1))
plot(moto.bgp,layout="surf")


X = seq(0, 20, length=100)
XX = seq(0, 20, length=49)
Z = (sin(pi*X/5) + 0.2*cos(4*pi*X/5)) * (X <= 9.6)
lin = (X > 9.6)
Z[lin] = -1 + X[lin]/10
Z = Z + rnorm(length(Z), sd=0.1)
plot(X,Z)

sine.bgp = bgp(X,Z,XX,BTE = c(500, 2000, 1),pred.n=F)
plot(sine.bgp,layout="surf")

sine.krige = predict(sine.bgp, XX=XX, pred.n=F)
plot(sine.krige, center="km", as="ks2")

sine.bgp2 = bgp(X,Z,XX,BTE = c(500, 2000, 1),d.p=c(20,10,20,10),pred.n=F,verb=2)
plot(sine.bgp2,layout="surf")


#install.packages("akima")
library(akima)
scallops=read.table("ch6.dat")
image(interp(scallops[,3],scallops[,4],log(scallops[,5]+1)))
persp(interp(scallops[,3],scallops[,4],log(scallops[,5]+1)),theta=-25,phi=20)
scallops.bgp = bgp(scallops[,c(3,4)],log(scallops[,5]+1),verb=2)
plot(scallops.bgp)

persp(interp(scallops[,3],scallops[,4],log(scallops[,5]+1)),theta=-25,phi=20)
plot(scallops.bgp,layout="surf")