library(fields)



# set up the data set to be easier to type
  data(COmonthlyMet)
  good<- !is.na(CO.tmin.MAM.climate)
  xHW1<- as.matrix(CO.loc[good,])
  yHW1<- CO.tmin.MAM.climate[good]
  elevHW1<- CO.elev[good]

x.star<- rbind( c(-105.2708,40.2247) )
data(RMelevation)
elev.star<- interp.surface(RMelevation, x.star)
# fiting exponential  (nu =.5)
theta.grid<- exp(seq(log(.15),log(10),,50))
MLE.fit0<-mKrig.MLE( xHW1, yHW1,Z= elevHW1,Covariance="Matern",smoothness=.5,
                   par.grid=list( theta=theta.grid))
fields.style()
plot( theta.grid, MLE.fit0$summary[,2], xlab="theta", log="x", ylab="ln Profile Like")
yline( max(MLE.fit0$summary[,2]) - qchisq(.95,1)/2 , lwd=2, col="grey20", lty=2)
xline( MLE.fit0$cov.args.MLE$theta, lwd=3, col="grey")
fit0<- mKrig(  xHW1, yHW1, Z= elevHW1,Covariance="Matern",smoothness=.5,
                            cov.args= MLE.fit0$cov.args.MLE, lambda= MLE.fit0$lambda.MLE)


##   dev.copy2pdf(file="/Users/nychka/Home/Current_talks/SpatialStatsCU/Lecture7/pix/like0.pdf")

# some other plots
plot( MLE.fit0$summary[,1], MLE.fit0$summary[,2], xlab="eff df", ylab="ln Profile Like")
yline( max(MLE.fit0$summary[,2]) - qchisq(.95,1)/2 , lwd=2, col="grey20", lty=2)

quilt.plot( cbind((MLE.fit0$summary[,4]), log(MLE.fit0$summary[,5])), MLE.fit0$summary[,2]  )


### search over  both theta and smoothness -- this takes a while to run.

theta.grid<- exp(seq(log(.15),log(1.5),,30))
sm.grid<- seq( .5,2.0,,15)
par.grid<- make.surface.grid(
                        list( smoothness=sm.grid, theta=theta.grid))
 lambda.start<- rep( c( .5,rep( NA, 30-1)), 15)                            
MLE.fit1<-mKrig.MLE( xHW1, yHW1,Z= elevHW1,
                         lambda=lambda.start,Covariance="Matern",
                         par.grid=par.grid, relative.tolerance=1e-6)
fields.style()
par( mar=c(4,4,1,1))
out.p<-  as.surface( par.grid, MLE.fit1$summary[,2])
surface( out.p, log="y")
contour( out.p, levels=c( max(out.p$z)- qchisq(.95,2)/2), col="grey80", lwd=3, add=TRUE)

theta.MLE.index<- rep( NA, 15)
for( k in 1:15){
theta.MLE.index[k]<- which.max(out.p$z[k,])
}
lines( sm.grid, theta.grid[theta.MLE.index], col="magenta", lwd=3)
xlines( sm.grid[c(1,8,15)], col="magenta", lty=3)
##   dev.copy2pdf(file="/Users/nychka/Home/Current_talks/SpatialStatsCU/Lecture7/pix/like1.pdf")

# comparing prediction surfaces
# refine estimate of lambda for these parameter values
theta1<- theta.grid[theta.MLE.index[1]]
sm1<- sm.grid[1]
llambda1<- mKrig.MLE(xHW1, yHW1, Z= elevHW1,lambda=.2,
                     theta= theta1, smoothness=sm1,Covariance="Matern",
                     relative.tolerance= 1e-6)$summary[,6]                      
theta2<- theta.grid[theta.MLE.index[8]]
sm2<- sm.grid[8]
llambda2<- mKrig.MLE(xHW1, yHW1, Z= elevHW1,lambda=.2,
                     theta= theta2, smoothness=sm2,Covariance="Matern",
                     relative.tolerance= 1e-6)$summary[,6]                      
theta3<- theta.grid[theta.MLE.index[15]]
sm3<- sm.grid[15]
llambda3<- mKrig.MLE(xHW1, yHW1, Z= elevHW1,lambda=.2,
                     theta= theta3, smoothness=sm3,Covariance="Matern",
                     relative.tolerance= 1e-6)$summary[,6]
# Kriging and OLS
fit1<- mKrig( xHW1, yHW1, Z= elevHW1,Covariance="Matern",
              smoothness=sm1, theta=theta1,
              lambda=exp(llambda1))
fit2<- mKrig( xHW1, yHW1, Z= elevHW1,Covariance="Matern",
              smoothness=sm2, theta=theta2,
              lambda=exp(llambda2))
fit3<- mKrig( xHW1, yHW1, Z= elevHW1,Covariance="Matern",
              smoothness=sm3, theta=theta3,
              lambda=exp(llambda3))
OLS.coef<- coefficients(lm( yHW1 ~ xHW1 + elevHW1))
temp<-cbind( OLS.coef,fit1$d, fit2$d, fit3$d)
dimnames(temp)<- list( c("intercept", "lon", "lat", "elev"), c("OLS", ".5", "1.25", "2.0"))
print(round(temp,5))

# predict on grids (note grids are created automatically but will be 150X150.
sur1<- predict.surface( fit1, drop.Z=TRUE, nx=150, ny=150)
sur2<- predict.surface( fit2, drop.Z=TRUE, nx=150, ny=150)
sur3<- predict.surface( fit3, drop.Z=TRUE, nx=150, ny=150)
# predictions under different models. 
fields.style()
set.panel(1,3)
zr<- range( sur1$z, na.rm=TRUE)
coltab<- tim.colors(256)
par( mar=c(1,1,1,2), oma=c(0,0,0,4))
image( sur1, axes=FALSE, xlab="", ylab="", col=coltab, zlim =zr)
title( ".5")
US( add=TRUE)
image( sur1, axes=FALSE, xlab="", ylab="", col=coltab, zlim =zr)
title( "1.25")
US( add=TRUE)
image( sur1, axes=FALSE, xlab="", ylab="", col=coltab, zlim =zr)
title( "2.0")
US( add=TRUE)
par( oma=c(0,0,0,1))
image.plot( legend.only=TRUE, col=coltab, zlim =zr)

## dev.copy2pdf(file="/Users/nychka/Home/Current_talks/SpatialStatsCU/Lecture7/pix/fits.pdf")
##
# looking at difference surfaces
 zmean<- (sur1$z + sur2$z + sur2$z)/3
fields.style()
set.panel(1,3)
zr<- range( sur1$z-zmean , na.rm=TRUE)
coltab<- tim.colors(256)
par( mar=c(1,1,1,2), oma=c(0,0,0,4))
image( sur1$x, sur1$y, sur1$z-zmean , axes=FALSE, xlab="", ylab="", col=coltab, zlim =zr)
title( ".5")
points( xHW1, pch=".",col="black")
US( add=TRUE)
image( sur1$x, sur1$y, sur2$z-zmean , axes=FALSE, xlab="", ylab="", col=coltab, zlim =zr)
title( "1.25")
points( xHW1, pch=".",col="black")
US( add=TRUE)
image(sur1$x, sur1$y, sur3$z-zmean ,  axes=FALSE, xlab="", ylab="", col=coltab, zlim =zr)
title( "2.0")
US( add=TRUE)
points( xHW1, pch=".",col="black")
par( oma=c(0,0,0,1))
image.plot( legend.only=TRUE, col=coltab, zlim =zr)
## dev.copy2pdf(file="/Users/nychka/Home/Current_talks/SpatialStatsCU/Lecture7/pix/resids.pdf")

# table of key parameters 
temp<- cbind(
       c(fit1$args$smoothness, fit1$args$theta, fit1$rho.MLE, fit1$sigma.MLE, fit1$eff.df),
       c(fit2$args$smoothness, fit2$args$theta, fit2$rho.MLE, fit2$sigma.MLE, fit2$eff.df),
       c(fit3$args$smoothness, fit3$args$theta, fit3$rho.MLE, fit3$sigma.MLE, fit3$eff.df))
dimnames( temp)<- list( c("smooth", "theta", "rho", "sigma", "eff.df"), format( temp[1,]))
print( round(temp[-1,],4))

# finding gamma

fit2<- mKrig( xHW1, yHW1, Z= elevHW1,Covariance="Matern",
              smoothness=sm2, theta=theta2,
              lambda=exp(llambda2))
#
# sanity check using Krig 
#fit.test<- Krig( xHW1, yHW1, Z= elevHW1,Covariance="Matern",
#              smoothness=sm2, theta=theta2,
#              lambda=exp(llambda2))
#
x.star<- rbind( c(-105.2708,40.2247) )
data(RMelevation)
elev.star<- interp.surface(RMelevation, x.star)
nobs<- length( yHW1)
gamma.lyons<- c(predict( fit2, x.star,Z= elev.star, ynew= diag(1,nobs)))

M<- fit2$rho.MLE* Matern( rdist( xHW1, xHW1)/theta2, smoothness=sm2) +
                                    fit2$sigma.MLE^2 * diag( 1, nobs)
kstar<- c(fit2$rho.MLE* Matern( rdist( xHW1, x.star)/theta2, smoothness=sm2))
# the formula for the standard error
SE<- sqrt(
          t(gamma.lyons) %*%M %*%gamma.lyons - 2* t(gamma.lyons) %*%kstar + fit2$rho.MLE
           )
predict.se(  fit2, x.star,Z= elev.star)
#
# plot of gamma weights
dev.off()
par( mar=c(1,1,1,2))
quilt.plot( xHW1, gamma.lyons, axes=FALSE, xlab="", ylab="", col=rainbow(256))
points( x.star, pch=16 , cex=1.5)
US( add=TRUE, lwd=2)
## dev.copy2pdf(file="/Users/nychka/Home/Current_talks/SpatialStatsCU/Lecture7/pix/gamma.pdf")
             
# creating a 150X150 grid on the range of the data
glist<- fields.x.to.grid( xHW1, nx=150, ny=150)

xg<- make.surface.grid( glist)
data( RMelevation)
elevg<- interp.surface( RMelevation, xg)

sur2<- predict( fit2, xg, Z= elevg)
SE2<- predict.se( fit2, xg, Z= elevg)

# plot of prediction surface and the standard errors. 
fields.style()
set.panel( 1,2)
par( mar=c(1,1,3,2))
image.plot( as.surface( glist, sur2), axes=FALSE, xlab="", ylab="", legend.mar=7.1)
title("Predicted field")
image.plot( as.surface( glist,SE2), axes=FALSE, xlab="", ylab="", legend.mar=7.1)
title("Standard errors")
points( xHW1, col="magenta")
##dev.copy2pdf(file="/Users/nychka/Home/Current_talks/SpatialStatsCU/Lecture7/pix/SE.pdf")