#Necessary libraries
library(INLA)
library(geoR)
set.seed(022017)

#set.seed(231171)


#Size of grid (note you might want to reduce this if the computation is slow
nrow=20
ncol=30
n = nrow*ncol

#Indices in both directions
i.m = matrix(rep(1:ncol,each=nrow),nrow=nrow)
j.m = matrix(rep(1:nrow,ncol),ncol=ncol)
#Covariate
x.m = sin(j.m+0*i.m)
image(x.m,zlim=c(-2,2),col=rainbow(256))

#Converting matrices to vectors, inla requires input to be vectors
x = inla.matrix2vector(x.m)
i = inla.matrix2vector(i.m)
j = inla.matrix2vector(j.m)

#Parameter values
beta = c(2,0.3);sigma.y=1
#Making covariance matrix through distance matrix
d = as.matrix(dist(cbind(i,j),diag=TRUE,upper=TRUE))
C = sigma.y^2*matern(d,2,3)
par(mfrow=c(1,1))
image(beta[1]+beta[2]*x.m,zlim=c(-1,5),col=rainbow(256))
image(C,zlim=c(-2,2),col=rainbow(256),ylim=rev(c(0 ,1)))


#Simulating data using Cholesky decomposition
C.chol = t(chol(C))
y = beta[1]+x*beta[2] + C.chol%*%matrix(rnorm(n),ncol=1)
#Additive Gaussian observation noise
z = y+0.3*rnorm(n)

y.m = inla.vector2matrix(y,nrow,ncol)
z.m = inla.vector2matrix(z,nrow,ncol)

#Plotting covariate and latent process and observations
par(mfrow=c(2,3))
image(beta[1]+beta[2]*x.m,zlim=c(-1,5),col=rainbow(256))

image(y.m,zlim=c(-1,5),col=rainbow(256))
image(z.m,zlim=c(-1,5),col=rainbow(256))

#Making dataframe for inla call, using delta for spatially correlated
#random effect
GridIndex = 1:n
data=data.frame(z=z,x=x,GridIndex=GridIndex)

#Defining model through a formula statement
#Note that the matern2d model assumes GridIndex is defined on a grid with
#dimensions nrow and ncol

formula= z ~ 1 + x +  f(GridIndex, model="matern2d", nu=3, nrow=nrow, ncol=ncol,hyper=c(1,1))

#inla call using empirical Bayes approach
result.eb=inla(formula, family="gaussian", data=data,
  control.predictor = list(compute = TRUE),
  control.inla=list(int.strategy="eb"))
summary(result.eb)

#inla call using Bayes approach
result.b=inla(formula, family="gaussian", data=data, control.predictor = list(compute = TRUE))
summary(result.b)

#Calculate predicted values and plotting it
pred.y.eb = result.eb$summary.fixed[1,1]+result.eb$summary.fixed[2,1]*x+
            result.eb$summary.random$GridIndex$mean
pred.y.m.eb = inla.vector2matrix(pred.y.eb,nrow,ncol)
image(pred.y.m.eb,zlim=c(-1,5),col=rainbow(256))

pred.y.b = result.b$summary.fixed[1,1]+result.b$summary.fixed[2,1]*x+
           result.b$summary.random$GridIndex$mean
pred.y.m.b = inla.vector2matrix(pred.y.b,nrow,ncol)
image(pred.y.m.b,zlim=c(-1,5),col=rainbow(256))

## poisson data
intensity=y-1;  # changed just to reduce the intensity slightly
z = rpois(n,exp(intensity))
z.m = inla.vector2matrix(z,nrow,ncol)

data=data.frame(z=z.m,x=x,GridIndex=GridIndex)
formula= z ~ 1 + x +  f(GridIndex, model="matern2d", nu=3, nrow=nrow, ncol=ncol
result.eb.pois=inla(formula, family="poisson", data=data,  control.inla=list(int.strategy="eb"))

summary(result.eb.pois)

pred.y.eb.pois = result.eb.pois$summary.fixed[1,1]+result.eb.pois$summary.fixed[2,1]*x+
            result.eb.pois$summary.random$GridIndex$mean
pred.y.m.eb.pois = inla.vector2matrix(pred.y.eb.pois,nrow,ncol)
par(mfrow=c(2,2))
image(intensity,zlim=c(-2,5),col=rainbow(256),ylim=rev(c(0 ,1)))
image(z.m,zlim=c(-1,15),col=rainbow(256),ylim=rev(c(0 ,1)))
image(pred.y.m.eb.pois,zlim=c(-2,5),col=rainbow(256),ylim=rev(c(0 ,1)))