#Tutorial-10
n<-20
N<-n^2
x<-expand.grid((1:n)/n,(1:n)/n)
plot(x[,1],x[,2])
dismat<-as.matrix(dist(x))
cov<-exp(-10*dismat^2)
mean<-rep(0,N)
library(MASS)
y<-mvrnorm(1,mean,cov)
library(lattice)
wireframe(matrix(y,20,20),drape=TRUE)
wireframe(matrix(exp(y),20,20),drape=TRUE)

#generate 100 realizations from a log-Gaussian process and look at the mean and 
#sd surface
M<-1000
y<-mvrnorm(M,mean,cov)
z<-exp(y)
zmean<-c()
ymean<-c()
zsd<-c()
ysd<-c()
for (i in 1:N) {
  zmean[i]<-mean(z[,i]) 
  ymean[i]<-mean(y[,i])
  zsd[i]<-sqrt(var(z[,i]))
  ysd[i]<-sqrt(var(y[,i]))
  }
wireframe(matrix(zmean,n,n),drape=TRUE)
wireframe(matrix(zsd,n,n),drape=TRUE)
wireframe(matrix(ymean,n,n),drape=TRUE)
wireframe(matrix(ysd,n,n),drape=TRUE)


#we will select a training set to predict a Gaussian surface based on
y<-mvrnorm(1,mean,cov)
xtrain<-cbind(1:20,sample(1:20,20))/20
#xtrain<-expand.grid(seq(1,20,2),seq(1,20,2))/20
xtrain<-as.matrix(xtrain)
nt<-dim(xtrain)[1]
ind<-c()
for (j in 1:nt) {
  for (i in 1:400) {
    
    if (x[i,1]==xtrain[j,1] & x[i,2]==xtrain[j,2]) ind[j]<-i
    
  }
}

ind<-sort(ind)
ord<-order(ind)
xtrain<-xtrain[ord,]
ytrain<-y[ind]
plot(xtrain[,1],xtrain[,2],type="n")
text(xtrain[,1],xtrain[,2],round(ytrain,digits=2))
library(scatterplot3d)
scatterplot3d(xtrain[,1],xtrain[,2],ytrain,highlight=TRUE,type="h",pch=16)

#our test set, i.e., where we are going to make prediction is the grid we started with
x<-as.matrix(x)
dist1<-as.matrix(dist(rbind(xtrain,x)))
bigcov<-exp(-10*dist1)+diag(rep(.001,N+nt))
R<-bigcov[1:nt,1:nt]#+diag(rep(.001,nt))
Rinv<-solve(R)
c<-bigcov[1:nt,(nt+1):(N+nt)]
M<-bigcov[(nt+1):(N+nt),(nt+1):(N+nt)]
yhat<-t(c)%*%Rinv%*%ytrain
wireframe(matrix(y,n,n),drape=TRUE)
wireframe(matrix(yhat,n,n),drape=TRUE)
MSE<-mean((yhat-y)^2)
ersurf<-matrix((yhat-y)^2,n,n)
wireframe(ersurf,drape=TRUE)

yhat[ind]