rm(list=ls())

library(spBayes)
library(MASS)
library(GenKern)

set.seed(1234)

########################
#make some data
########################
n <- 50
x <- runif(n, 0, 100)
y <- runif(n, 0, 100)

D <- as.matrix(dist(cbind(x, y)))

phi <- 3/100
sigmasq <- 50
tausq <- 20
mu <- 150

s <- (sigmasq*exp(-phi*D))
w <-  mvrnorm(1, rep(0, n), s)
Y <- mvrnorm(1, rep(mu, n) + w, tausq*diag(n))
X <- as.matrix(rep(1, length(Y)))


lmObject <- lm(Y ~ X)
lmObject.cov <- vcov(lmObject)

write(t(Y), "Ytest.txt", ncolumns=1)
write(t(X), "Xtest.txt", ncolumns=1)
write(t(cbind(x,y)), "Sitestest.txt", ncolumns=2)

##Create Dataframe for all data

simData <- as.data.frame(cbind(x, y, X, Y))
names(simData) <- c("x", "y", "X", "Y")

########################
#some priors for sigma^2
#tau^2 and phi
########################
#IG sigma^2 and tau^2
a.sig <- 2 
b.sig <- 100
a.tau <- 2
b.tau <- 100

#U phi
b.phi <- 3/3
a.phi <- 3/150


##Specify the priors, hyperparameters, and variance parameter starting values.
sigmasq.prior <- prior(dist="IG", shape=a.sig, scale=b.sig)
tausq.prior <- prior(dist="IG", shape=a.tau, scale=b.tau)
##the prior on phi corresponds to a prior of 3-150 meters
##for the effective range (i.e., -log(0.05)/3, -log(0.05)/150, when
##using the exponential covariance function).
phi.prior <- prior(dist="UNIF", a=a.phi, b=b.phi)

var.update.control <- list("sigmasq"=list(sample.order=0, starting=30, tuning=0.25, prior=sigmasq.prior),
                           "tausq"=list(sample.order=0, starting=30, tuning=0.25, prior=tausq.prior),
                           "phi"=list(sample.order=0, starting=0.006, tuning=0.25, prior=phi.prior))

##specify the number of samples
run.control <- list("n.samples" = 50000)

##specify some model, assign the prior and starting values for the regressors
model <- Y~1

##note, precision of 0 is same as flat prior.
beta.prior <- prior(dist="NORMAL", mu=rep(0, 1), precision=diag(0, 1))

beta.control <- list(beta.update.method="gibbs", beta.starting=rep(0, 1), beta.prior=beta.prior)

model.coords <- cbind(x,y)

ggt.sp.out <-ggt.sp(formula=model,
                    data=simData,
                    coords=model.coords,
                    var.update.control=var.update.control,
                    beta.control=beta.control,
                    cov.model="exponential",
                    run.control=run.control, DIC=TRUE, DIC.start=30001, verbose=TRUE)

##print(ggt.sp.out$accept)
print(ggt.sp.out$DIC)


##Get the effective range of spatial dependence.
eff.range <- -log(0.05)/ggt.sp.out$p.samples[,"phi"]

##Plot the chains with the mcmc function in coda.
mcmc.obj <- mcmc(cbind(ggt.sp.out$p.samples, eff.range))
##plot(mcmc.obj)

##Summarize the post burn-in posterior distributions
##We use the "start" option in the mcmc function in coda
mcmc.obj.post.burnin <- mcmc(cbind(ggt.sp.out$p.samples, eff.range), start=30001)
print(summary(mcmc.obj.post.burnin))


##Customized summary

myqnt <- c(0.50, 0.025, 0.975)
start <- 30001
beta.hat <- quantile(ggt.sp.out$p.samples[,"(Intercept)"], myqnt)
sigmasq.hat <- quantile(ggt.sp.out$p.samples[, "sigmasq"], myqnt)
tausq.hat <- quantile(ggt.sp.out$p.samples[,"tausq"], myqnt)
phi.hat <- quantile(ggt.sp.out$p.samples[,"phi"], myqnt)

beta.hat
sigmasq.hat
tausq.hat
phi.hat

##pdf("ConvergencePlots.pdf")
par(mfrow=c(2,2))
plot(ggt.sp.out$p.samples[start:nrow(ggt.sp.out$p.samples),"(Intercept)"],type="l",main="mu",ylab="",xlab="samples")
abline(h=mu, col="red")
plot(ggt.sp.out$p.samples[start:nrow(ggt.sp.out$p.samples),"sigmasq"],type="l",main="sigma^2",ylab="",xlab="samples")
abline(h=sigmasq, col="red")
plot(ggt.sp.out$p.samples[start:nrow(ggt.sp.out$p.samples),"tausq"],type="l",main="tau^2",ylab="",xlab="samples")
abline(h=tausq, col="red")
plot(3/ggt.sp.out$p.samples[start:nrow(ggt.sp.out$p.samples),"phi"],,type="l",main="3/phi",ylab="",xlab="samples")
abline(h=3/phi, col="red")
##dev.off()


##############################################
##recover random spatial effect
##Y = X^tB + w + e, where w is the random
##spatial effect and e is pure error
##############################################

##get the spatial random effect w
sp.effect.out <- sp.effect(ggt.sp.out, start=start, thin=1)
w <- rowMeans(sp.effect.out$pp.samples)

par(mfrow=c(1,2))
int.obs <- interp(x, y, Y)
image(int.obs, xlab="x", ylab="y", main="Observed Y")
contour(int.obs, add=TRUE)

int.w <- interp(x, y, w)
image(int.w, xlab="x", ylab="y", main="Spatial random effects (w)")
contour(int.w, add=TRUE)

###############################################
##prediction
###############################################

##resample
n.pred <- 35

simData.pred <- simData[sample(1:nrow(model.coords), n.pred),]
pred.coords <- cbind(simData.pred$x, simData.pred$y)

pred <- sp.predict(ggt.sp.out, pred.joint=TRUE, pred.coords=pred.coords, pred.covars=simData.pred, start=start)

par(mfrow=c(1,2))
int.obs <- interp(x,y, Y)
image(int.obs, xlab="x", ylab="y", main="Observed Y")
contour(int.obs, add=TRUE)

int.pred <- interp(simData.pred$x, simData.pred$y, rowMeans(pred$pp.samples))
image(int.pred, xlab="x", ylab="y", main="Predicted Y")
contour(int.pred, add=TRUE)