# Load Packages
library(maps)          #for state borders
library(ggplot2)      # for plotting spatial map
library(extRemes)     #for fit the GEV
library(fExtremes)    # for getting the GEV quantiles
library(fields)       # for surface to plot surface plot
library(geoR)         # for use as.geodata
library(spBayes)     # for bayesian kriging
library(latex2exp)    #  LaTeX math formulas to R's
library(RColorBrewer) # Color schemes for maps
library(reshape2)    #to use the melt function
library(data.table)   #Fast aggregation of large data
library(gstat)        # to Create gstat objects and variogram()
library(sp)           # for using coordinates

# #########multiplot function ---------------------------------------------
multiplot = function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  
  numPlots = length(plots)
  
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  
  if (numPlots==1) {
    print(plots[[1]])
    
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}


get.samples_1 <- function(Y,X,X1,name,n.samples,burn,thin1,displ=FALSE){
  #Fit variogram and perform MCMC
  # obtain the variogram and the starting values of beta
  bestmod= lm(Y ~ ., data = X)#fit a linear regression and get the residuals
  geod = as.geodata(cbind(bestmod$residuals,X$Long,X$Lat),coords.col=2:3,data.col=1)
  break1=seq(0,50,0.5)#the range of distance for computing the variogram
  vg = variog(geod,breaks=break1,messages=displ)
  tmp=summary(bestmod)
  
  #starting values
  beta_st=tmp$coefficients[,1]
  sigma.sq = median(vg$v)-min(vg$v)#initial sill value
  tau.sq=min(vg$v)#initial nugget value
  phi.val = median(vg$u)#initial expoential decay value
  
  #weekly informative priors
  sigma.sq_prior=c(10,0.1)# shape and scale parameter
  tau.sq_prior=c(10,0.1)# shape and scale parameter
  phi_prior=c(0.1,15)# unif(min dist,max dist)
  beta_prior=list(rep(0,length(beta_st)), diag(1000,length(beta_st)))# meann 0 and variance=1000 for each coefficients
  
  # prior distributions defined (limited distributions available... see help(spLM)
  priors = list("beta.norm"=beta_prior, "phi.unif"=phi_prior,
                "sigma.sq.ig"=sigma.sq_prior,"tau.sq.IG"=tau.sq_prior)
  # starting values for MCMC resampling
  starting = list("phi"=phi.val, "sigma.sq"=sigma.sq,"beta"=beta_st, "tau.sq"=tau.sq)
  # adjustment factor for MCMC sampling routine
  tuning = list("phi"=0.1, "sigma.sq"=0.1,"tau.sq"=0.1) 
  
  # Perform Monte Carlo Markov Chain Analysis
  bat.fm = spLM(Y~., data=X, coords=as.matrix(X[,1:2]), priors=priors, tuning=tuning, starting=starting, cov.model = "exponential", n.samples = n.samples, verbose = displ, n.report = n.samples/2)
  #add some small amount to duplicated coordinates
  dup = duplicated(bat.fm$coords); bat.fm$coords[dup] <- bat.fm$coords[dup] + 1e-3
  # burn-in samples
  bat.fm = spRecover(bat.fm, start=((1/2)*n.samples)+1, thin=thin1, verbose=displ)
  y1=spPredict(bat.fm,pred.coords = as.matrix(X1[,1:2]),as.matrix(cbind(rep(1,nrow(X1)),X1)), start=((1/2)*n.samples)+1, thin=thin1, verbose=displ)
  y=y1$p.y.predictive.samples
  
  median.y = apply(y, 1, FUN = median) #median
  median.yse = apply(y, 1, FUN = sd)/sqrt(nrow(y)) #standard error
  #plot of the posterior distribution of the coefficients
  p=ggplot() +
    geom_histogram(aes(as.vector(y),..density..),fill="gray50",col='black') +
    labs(title = TeX(paste('Histogram of Posterior $\\',name,'$',sep = "")),
         x = "Predictions")+
    #theme_light()+
    theme(plot.title = element_text(size = 13, face = "bold", hjust = 0.5))
  #median and standard error
  Pred=as.data.frame(cbind(median.y,median.yse))
  names(Pred)=c("fit","se")
  #saving all the importance information in a list
  result=list(y=y,Pred=Pred,p=p,acep_r=bat.fm$acceptance)
  # rm(fig.base)
  result
}


# plots of the median coefficients over a grid

plot_sp_grid=function(Pred,Y,X1,xgrid,ygrid,nglobe,index1,name,xlim1,ylim1,namex,namey,labx,laby,r_pal,outl=TRUE){
  
  par( oma=c(0,0,0,3)) # save some room for the legend
  par(mfrow = c(1, 3))
  par(mar = c(4.5, 6, 3, 1.5))
  if (outl){
    zlim1=range(Pred$fit,Y)
  }else{
    zlim1=range(Pred$fit)}
  # ,zlim=zlim1
  quilt.plot(X1[,1], X1[,2],Y,axes=F,xlab="Longitude",ylab="Latitude",main=TeX(paste('mle $\\',name,'$',sep="")),zlim=zlim1,ylim=ylim1,xlim=xlim1,col=r_pal,cex.main=1.4,cex.lab=1.3)
  # grid(col="gray70",lty=2)
  US(add=TRUE, col="gray30", lwd=2)
  box()
  axis(1, at=namex, labels=labx)
  axis(2, at=namey, labels=laby)
  
  zfull = rep(NaN,nglobe)
  zfull[index1]=Pred$fit
  zmat = matrix(zfull,nrow=length(xgrid),ncol=length(ygrid))
  image.plot(xgrid,ygrid,zmat,ylim=ylim1,xlim=xlim1,main =TeX(paste('Median Posterior $\\',name,'$ on DEM Grid',sep="")),axes=F,xlab="Longitude",ylab="Latitude",cex.main=1.4,cex.lab=1.3,col=r_pal)
  # contour(xgrid,ygrid,zmat,ylim=ylim1,xlim=xlim1,add=TRUE,nlev=4,lwd=1)
  US(add=TRUE, col="gray30", lwd=2)
  # grid(col="gray30",lty=2)
  box()
  axis(1, at=namex, labels=labx,cex.axis=1)
  axis(2, at=namey, labels=laby,cex.axis=1)
  
  zfull = rep(NaN,nglobe)
  zfull[index1]=Pred$se
  zmat = matrix(zfull,nrow=length(xgrid),ncol=length(ygrid))
  image.plot(xgrid,ygrid,zmat,ylim=ylim1,xlim=xlim1,main =TeX(paste('Median Posterior $\\',name,'$ on DEM Grid',sep="")),axes=F,xlab="Longitude",ylab="Latitude",cex.main=1.4,cex.lab=1.3,col=r_pal)
  # contour(xgrid,ygrid,zmat,ylim=ylim1,xlim=xlim1,add=TRUE,nlev=4,lwd=1)
  US(add=TRUE, col="gray30", lwd=2)
  # grid(col="gray30",lty=2)
  box()
  axis(1, at=namex, labels=labx,cex.axis=1)
  axis(2, at=namey, labels=laby,cex.axis=1)
  
  
}