library(fields)


#
# multivariate normal example
 mu<- matrix( c( .1 , -.5), ncol=1, nrow=2)
  Sigma<-  rbind(  c(1,.8), 
                           c(.8,1) )
  v.grid<- make.surface.grid( 
                    list( x= seq( -3,3,,100), y =   seq( -3,3,,100) ))                   
  z<- rep( NA, 100*100)
  for( k in 1:(100*100)){
          v<-  v.grid[k,] # a vector with 2 rows           
          z[k] <-  ( 1/ (2*pi)^( 2/2) ) *
                exp(  - t(v - mu) %*% solve( Sigma) %*% (v - mu)/2 ) *
                                   det( Sigma)^(-2/2)                
   }   
   surface( as.surface( v.grid, z))    

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

# 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]

# two handy functions:
which.max.image<- function (obj) 
{
    ind.z <- which.max.matrix(obj$z)
    return(list(x = obj$x[ind.z[, 1]], y = obj$y[ind.z[, 2]], 
        z = obj$z[ind.z], ind = ind.z))
}
which.max.matrix<- function (z) 
{
    if (!is.matrix(z)) {
        stop("Not a matrix")
    }
    m <- nrow(z)
    n <- ncol(z)
    ind <- which.max(z)
    iy <- trunc((ind - 1)/m) + 1
    ix <- ind - (iy - 1) * m
    return(cbind(ix, iy))
}

  # grid over theta and loglambda using exponential covariance
  llambda.grid<- seq( log(.05), log( 5),,50)
  ltheta.grid<- seq( log(.2),log(4),,50)
  par.grid<- make.surface.grid( list( x= llambda.grid, y=ltheta.grid))
  NG<- nrow( par.grid)
  lgProfileLike<- rep( NA, NG)
  rho.profile<- rep( NA, NG)
  nobs<- nrow( xHW1)
  X<- cbind(rep(1, nobs), xHW1, elevHW1)                             
  for( k in 1:NG){
      cat(k ," ")
      theta.temp<- exp( par.grid[k,2])
      lambda.temp<- exp( par.grid[k,1])
      K<- exp( -rdist( xHW1, xHW1)/theta.temp)
      M<-  K + diag(lambda.temp,nobs)
      Minv<- solve( M)
      d.MLE<- solve( t(X)%*%Minv%*%X)%*%t(X)%*%Minv%*%yHW1
      residual<- yHW1 - X%*%d.MLE
      rho.MLE<- (1/nobs)*t(residual)%*%(Minv%*%residual)
      rho.profile[k] <- rho.MLE
      lgProfileLike[k]<- -nobs/2 *log( 2*pi) - (nobs/2) -
                         (1/2)*determinant(M, logarithm=TRUE)$modulus - (nobs/2)*log(rho.MLE)
### check function from fields
###     out<-  mKrig( xHW1, yHW1, theta=theta.temp,Z=elevHW1, lambda=lambda.temp)
    }
 lgPL1<-  as.surface( par.grid, lgProfileLike)
 surface(lgPL1,ylab=" log theta", xlab="log lambda")
# find maximum onthis grid
 maxLike<- which.max.image(lgPL1)
 ltheta.MLE<-  maxLike$y
 llambda.MLE<- maxLike$x
 lgPLmax<- maxLike$z
 points(  llambda.MLE,ltheta.MLE,col="grey", pch=16)
#  95% CI for log theta and log lambda
 
  contour(lgPL1 , level= c(lgPLmax - qchisq(.95, 2)/2 ), lwd =4 ,
         add= TRUE,labcex=.1, draw.labels=FALSE, col="grey")
  title("log Profile Likelihood with 95% confidence region (grey)")
##   dev.copy2pdf(file="/Users/nychka/Home/Current_talks/SpatialStatsCU/Lecture6/pix/like1.pdf")

  theta.MLE<- exp(ltheta.MLE)
  lambda.MLE<- exp( llambda.MLE)
# an interesting plot
   ind<-    lgProfileLike > lgPLmax - qchisq(.9, 2)/2
   temp<- cbind( exp(par.grid[ind,2]), rho.profile[ind], exp(par.grid[ind,1])*rho.profile[ind])
   dimnames( temp)<- list( NULL, c("theta","rho", "sigma^2"))
   pairs( temp)

  K<- exp( -rdist( xHW1, xHW1)/theta.MLE)
  M<-  K + diag(lambda.MLE,nobs)
  Minv<- solve( M)
  d.MLE<- solve( t(X)%*%Minv%*%X)%*%t(X)%*%Minv%*%yHW1
  residual<- yHW1 - X%*%d.MLE
  rho.MLE<- c((1/nobs)*t(residual)%*%(Minv%*%residual))
  sigma.MLE<- sqrt( lambda.MLE*rho.MLE)
  print( c(theta.MLE, lambda.MLE, rho.MLE))
  print( d.MLE)
  c.MLE<- Minv%*%residual
# out<-  mKrig( xHW1, yHW1, theta=theta.MLE,Z=elevHW1, lambda=lambda.MLE)
# predictions at Lyons, CO
#  40.2247° N, 105.2708° W

x.star<- rbind( c(-105.2708,40.2247) )
data(RMelevation)
elev.star<- interp.surface(RMelevation, x.star)

ghat<-   t(rho.MLE* exp( -rdist( xHW1,x.star)/theta.MLE)) %*% c.MLE
X.star<- matrix(c( 1, x.star, elev.star), nrow=1)

fixed.part<-    X.star%*%d.MLE

Lyons.prediction <- fixed.part + ghat
print( Lyons.prediction)